---
name: stationarity-gages-ii-reference-gauges
description: >
  Tests whether GAGES-II reference (hydrologically unimpaired) gauges exhibit
  flood non-stationarity at a different rate than regulated gauges on the same
  CONUS record window, as a natural experiment isolating climate-driven trend
  signals from flow-regulation artifacts. Downloads USGS NWIS annual peak
  streamflow for ~120 reference and ~40 regulated gauges, computes
  Hamed & Rao (1998) variance-corrected Mann-Kendall z-scores per gauge,
  block-bootstrap 95% confidence intervals for the median Sen slope, and a
  permutation test of the reference-minus-regulated rejection-rate gap under
  the null that GAGES-II reference status is uninformative. Sensitivity
  analyses span record-length thresholds (30/40/50/60 yr) and alternate
  reference-set inclusion rules. Python 3.8 standard library only.
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "hydrology", "flood-frequency", "stationarity", "hamed-rao", "mann-kendall", "block-bootstrap", "gages-ii", "usgs-nwis", "reference-vs-regulated"]
python_version: ">=3.8"
dependencies: []
---

# Do GAGES-II Reference Gauges Show Non-Stationarity When Restricted to Unregulated Basins?

## When to Use This Skill

Use this skill when you need to test whether an apparent non-stationary trend in
a time-series claim (here, CONUS flood magnitudes) is a genuine climate signal
or a reporting/confound artifact, using a negative-control design that contrasts
a putatively unimpaired subset against a regulated comparator. Concretely: use
it when the widely-repeated claim that floods are becoming larger across CONUS
needs to be tested on the GAGES-II reference (hydrologically unimpaired) gauge
subset with Hamed & Rao (1998) autocorrelation-corrected Mann-Kendall and a
label-permutation null. If only the regulated cohort shows non-stationarity,
"climate-driven" trend claims on the full network are confounded by flow
regulation.

## Research Question

Does the rejection rate of the Mann-Kendall stationarity null for annual peak
streamflow differ between GAGES-II reference (unimpaired) gauges and regulated
gauges on a matched CONUS record window, after variance correction for serial
dependence, and is any observed gap explainable by sampling variation in
group composition?

## Prerequisites

- **Python**: 3.8+ standard library ONLY (`math`, `hashlib`, `json`, `random`,
  `urllib.request`). No `pip install`; no `numpy`, `scipy`, or `pandas`.
- **Network**: Internet access to `waterdata.usgs.gov` and
  `nwis.waterdata.usgs.gov` on first run (~183 gauges, each < 100 KB; ~15 MB
  total). Subsequent runs use the SHA256-verified local cache.
- **Disk space**: < 50 MB (cached RDB files + results.json + report.md).
- **Runtime**: first run ~10–25 minutes (network + per-gauge Hamed-Rao +
  2,000-draw gauge bootstrap + 5,000-shuffle permutation, x2 variants).
  Cached reruns: ~2–5 minutes. `--verify` mode: < 2 seconds.
- **Environment variables**: none required.
- **Write access**: needed for the workspace directory
  (`/tmp/claw4s_auto_stationarity-in-gages-ii-reference-gauges/`), its `cache/`
  subdir, and the output files `results.json` and `report.md`.

## Adaptation Guidance

To adapt this skill to a different region, variable, or reference-set
definition, modify only the `DOMAIN CONFIGURATION` block of the script:

- **Change data source**: edit `NWIS_URL_TEMPLATES`, `USGS_REFERENCE_SITES`,
  and `USGS_REGULATED_SITES`. The parser `parse_usgs_rdb()` is specific to
  USGS RDB format — replace it for a different repository, but keep its
  return contract `List[(year:int, value:float)]`.
- **Change analyzed variable**: `load_data()` emits
  `{group_label: {site_id: {year: value}}}`; any annual-value time series
  of the same shape plugs directly into `run_analysis()`.
- **Change record-length sweep**: edit `RECORD_LENGTH_THRESHOLDS`. These
  drive the sensitivity sweep in `run_analysis()` — the per-threshold
  rejection-rate gap between groups is recomputed automatically.
- **Change reference-set definition**: edit `REFERENCE_SET_VARIANTS` to add
  or remove alternate site lists. Each variant re-enters `run_analysis()`
  in its own loop so every downstream rejection rate and CI reflects that
  variant.
- **Keep the statistical core unchanged**: `mann_kendall_S()`,
  `hamed_rao_variance_correction()`, `hamed_rao_pvalue()`, `sen_slope()`,
  `block_bootstrap_median_slope_ci()` (a gauge-level resampling CI, not
  a within-series block bootstrap), `label_permutation_test()`,
  `stratified_label_permutation_test()`, and `benjamini_hochberg()`
  are domain-agnostic and should never change when re-targeting.

## Overview

**Common claim**: CONUS mean annual peak streamflow is rising — a signal
widely cited as evidence that "stationarity is dead" (Milly et al. 2008,
*Science*).

**Confound tested**: Regulated gauges (reservoirs, diversions, impervious
upstream land) dominate naive CONUS aggregates. If the reported signal
collapses on the GAGES-II reference subset (Falcone 2011), the
meta-interpretation of a climate-wide shift is unsupported.

**Methodological hook**: We apply the Hamed & Rao (1998) variance
correction, the mainstream hydrologic remedy for autocorrelation-induced
p-value inflation in Mann-Kendall tests, and compute a *reference minus
regulated* rejection-rate gap. The gap is tested against a permutation
null that reassigns GAGES-II reference status uniformly at random across
the combined gauge pool, preserving record length and region. This null
controls for the possibility that the reference group is merely shorter,
smaller, or in a different climatic belt than the regulated group.

**Effect sizes**: The median Sen slope (cfs/yr and %/decade) per group
is reported with a 95% percentile confidence interval obtained by
resampling gauges with replacement within the group (2,000 draws). This
is a *between-gauge* uncertainty estimate, not a within-gauge
dependent-time-series bootstrap; the autocorrelation-adjusted Mann-Kendall
p-values for each gauge come from the Hamed-Rao variance inflation.

**Null model**: Three layers.

1. *Per-gauge trend test*: Under the null of stationarity, Mann-Kendall S
   has mean zero and a variance that must be inflated for serial
   dependence. The Hamed-Rao inflation factor `n/n*` uses only
   lag-k rank autocorrelations that are individually significant at
   α = 0.05 after detrending, so autocorrelation is not over-corrected.
2. *Unstratified label permutation*: Under the null that
   reference-vs-regulated status is uninformative, permuting the group
   labels across the pooled gauge set should produce rejection-rate gaps
   and median-slope differences at least as extreme as observed with
   probability p. We compute these p-values from 5,000 label shuffles.
3. *HUC2-stratified label permutation (supplementary)*: A tighter null
   that permutes labels only within 2-digit hydrologic unit codes,
   thereby preserving the regional composition of each group. Reports
   a second p-value; where the strata are too thin it falls back to
   the pooled permutation.

**Sensitivity**: Rejection rates recomputed at `MIN_RECORD_YEARS` ∈
{30, 40, 50, 60}, and separately under two reference-set variants
(strict: only gauges with GAGES-II classification = Ref; permissive:
Ref + hydro-climatically unimpaired non-Ref).

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_stationarity-in-gages-ii-reference-gauges/cache
```

**Expected output:** Directory created silently. Exit code 0.

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_stationarity-in-gages-ii-reference-gauges/analyze.py
#!/usr/bin/env python3
"""
GAGES-II reference-vs-regulated flood stationarity test using
Hamed & Rao (1998) variance-corrected Mann-Kendall and block-bootstrap
Sen-slope confidence intervals.

Python 3.8+ standard library only. No external dependencies.
"""

import sys
import os
import json
import math
import hashlib
import random
import time
import urllib.request
import urllib.error

# ═══════════════════════════════════════════════════════════════
# DOMAIN CONFIGURATION — To adapt this analysis to a new domain,
# modify only this section.
# ═══════════════════════════════════════════════════════════════

SEED = 42
WORKSPACE = os.path.dirname(os.path.abspath(__file__))
CACHE_DIR = os.path.join(WORKSPACE, "cache")
RESULTS_FILE = os.path.join(WORKSPACE, "results.json")
REPORT_FILE = os.path.join(WORKSPACE, "report.md")

USER_AGENT = "GagesStationarityHamedRao/1.0 (claw4s-2026; academic)"

# Core analysis parameters
MIN_RECORD_YEARS = 50        # primary inclusion threshold
MAX_RECORD_YEAR = 2023       # ignore partial or future years
MIN_ANNUAL_FOR_PARSE = 20    # gauges with < 20 years are dropped on load
ALPHA = 0.05                 # rejection level for MK and lag-k screening
FDR_Q = 0.10                 # Benjamini-Hochberg false discovery rate
N_BOOTSTRAP_GAUGES = 2000    # gauge-level resamples for group-median CI
N_PERMUTATIONS = 5000        # group-label permutations
STRATIFIED_MIN_STRATUM = 3   # HUC2 strata smaller than this are pooled

RECORD_LENGTH_THRESHOLDS = [30, 40, 50, 60]   # sensitivity sweep
REFERENCE_SET_VARIANTS = ["strict", "permissive"]  # reference-set sensitivity

REQUEST_TIMEOUT = 60
REQUEST_RETRIES = 3
INTER_REQUEST_SLEEP = 0.3

# Plausibility bounds used in --verify assertions.  Loose by design: they
# catch degenerate or corrupt output without flagging real-world extremes.
MAX_PLAUSIBLE_PCT_PER_DECADE = 20.0     # |median %/decade| must not exceed this
MIN_CI_WIDTH_PCT = 0.01                 # CI on median %/decade must not collapse
MAX_PLAUSIBLE_GAP_AT_30YR = 0.30        # |ref - reg rejection gap| at 30yr cap

# Limitations are echoed into results.json so downstream consumers see them
# alongside the numbers.
LIMITATIONS = [
    "Gauge coverage is a curated subset of GAGES-II, not all CONUS stations.",
    "Hamed-Rao assumes lag-k autocorrelation structure is itself stationary.",
    "Annual peaks are calendar-year maxima, not water-year maxima.",
    "GAGES-II reference classification (Falcone 2011) predates post-2011 land-use change.",
    "Unstratified permutation ignores record-length and HUC2 imbalance.",
    "Results do not attribute non-stationarity to any specific climate or LUC driver.",
]

NWIS_URL_TEMPLATES = [
    "https://nwis.waterdata.usgs.gov/nwis/peak?site_no={site}&agency_cd=USGS&format=rdb",
    "https://waterdata.usgs.gov/nwis/peak?site_no={site}&agency_cd=USGS&format=rdb",
]

# GAGES-II *reference* gauges (Falcone 2011): hydrologically unimpaired,
# curated across CONUS hydrologic regions. Long-record sites only.
USGS_REFERENCE_SITES_STRICT = [
    # Northeast
    "01013500", "01022500", "01030500", "01031500", "01047000",
    "01054200", "01055000", "01064500", "01073000", "01134500",
    "01137500", "01144000", "01162500", "01169000", "01170100",
    "01187300", "01321000", "01334500", "01350000", "01413500",
    "01414500", "01420500", "01435000",
    # Southeast
    "02016000", "02051500", "02064000", "02070000", "02074500",
    "02111000", "02128000", "02143040", "02177000", "02198000",
    "02215100", "02296750", "02310000", "02314500", "02329000",
    "02347500", "02361000",
    # Ohio / Tennessee
    "03010500", "03015500", "03049800", "03069500", "03118500",
    "03144000", "03164000", "03170000", "03186500", "03219500",
    "03281100", "03285000", "03302680", "03320000", "03455000",
    "03460000", "03473000", "03500000", "03504000",
    # Great Lakes
    "04040500", "04056500", "04059500", "04105000",
    "04122500", "04185000",
    # Upper Mississippi / Red-Souris
    "05056000", "05062000", "05082500", "05120500", "05291000",
    "05362000", "05389500", "05412500", "05420500", "05454300",
    "05495500",
    # Missouri
    "06191500", "06214500", "06280300", "06404000", "06441500",
    "06452000", "06622700", "06623800", "06784000", "06814000",
    "06892000",
    # Arkansas-Red / Lower Mississippi
    "07014500", "07019000", "07056000", "07068000",
    "07148400", "07196500",
    # Texas-Gulf
    "08070200", "08082700", "08171000", "08178880", "08189500",
    # Colorado / Great Basin / Rio Grande
    "09066300", "09081600", "09306242", "09430500", "09444500",
    "09484000", "09508500", "10234500", "10396000",
    # Pacific Northwest / California / Idaho
    "11124500", "11186000", "11253310", "11264500", "11266500",
    "11381500", "11475000", "11476500", "11532500",
    "12010000", "12035000", "12054000", "12134500", "12186000",
    "12413000", "12451000",
    "13185000", "13313000", "13317000", "13340600", "13337000",
    "14020000", "14154500", "14222500", "14301000", "14305500",
]

# Additional GAGES-II *non-reference but minimally disturbed* gauges:
# small impervious fraction, no major upstream storage, used only for
# the "permissive" reference-set sensitivity variant.
USGS_REFERENCE_SITES_ADDITIONAL = [
    "01389500", "02083500", "03237500", "04193500",
    "05418500", "06823500", "07263295", "08049500",
    "09379500", "11402000", "12056500", "13317660",
]

# Regulated / impaired comparator gauges (reservoirs, diversions, or
# heavy urban land use upstream), roughly region-matched. Never placed
# in the reference set regardless of variant.
USGS_REGULATED_SITES = [
    "01100000", "01184000", "01357500", "01428500",      # Northeast
    "02035000", "02087500", "02198500", "02326000",      # Southeast
    "03049500", "03086000", "03155500", "03247500",      # Ohio/TN
    "04085427", "04101500", "04121970",                   # Great Lakes
    "05288500", "05331000", "05420690", "05474000",      # UMiss
    "06192500", "06309000", "06468170",                   # Missouri
    "07022000", "07144100", "07263500",                   # ArkRed/LMiss
    "08057410", "08068500", "08211000",                   # TexasGulf
    "09163500", "09402000", "09522000", "10310000",      # Southwest
    "11303500", "11447650", "12462500", "13269000",      # Pacific
    "14103000", "14191000",
]

# HUC2 region membership (first two digits of USGS site_no identify the
# 2-digit hydrologic unit code, which is the coarse hydrographic region).
# We use this for region-stratified permutation of group labels so the
# permutation null preserves climatic composition.
# (No configuration needed — derived from site_no at runtime.)


# ═══════════════════════════════════════════════════════════════
# STATISTICAL HELPERS (domain-agnostic)
# ═══════════════════════════════════════════════════════════════

def normal_cdf(z):
    return 0.5 * (1.0 + math.erf(z / math.sqrt(2.0)))


def mean_val(x):
    return sum(x) / len(x) if x else 0.0


def median_val(x):
    s = sorted(x)
    n = len(s)
    if n == 0:
        return float("nan")
    if n % 2 == 1:
        return s[n // 2]
    return 0.5 * (s[n // 2 - 1] + s[n // 2])


def percentile_val(lst, p):
    s = sorted(lst)
    n = len(s)
    if n == 0:
        return float("nan")
    k = (p / 100.0) * (n - 1)
    lo = int(math.floor(k))
    hi = min(lo + 1, n - 1)
    frac = k - lo
    return s[lo] * (1.0 - frac) + s[hi] * frac


def rank_array(x):
    """Average-tie ranks (1-based). Matches the rank convention used in
    Hamed-Rao (1998)."""
    n = len(x)
    order = sorted(range(n), key=lambda i: x[i])
    ranks = [0.0] * n
    i = 0
    while i < n:
        j = i
        while j + 1 < n and x[order[j + 1]] == x[order[i]]:
            j += 1
        avg = (i + j) / 2.0 + 1.0
        for k in range(i, j + 1):
            ranks[order[k]] = avg
        i = j + 1
    return ranks


def autocorr_rank(ranks, lag):
    """Autocorrelation of ranks at given lag (Pearson on ranks = Spearman
    on original values when no trend is present)."""
    n = len(ranks)
    if lag >= n or lag < 1:
        return 0.0
    mu = mean_val(ranks)
    num = sum((ranks[i] - mu) * (ranks[i + lag] - mu)
              for i in range(n - lag))
    den = sum((r - mu) ** 2 for r in ranks)
    return num / den if den > 0 else 0.0


def mann_kendall_S(x):
    n = len(x)
    S = 0
    for i in range(n - 1):
        xi = x[i]
        for j in range(i + 1, n):
            d = x[j] - xi
            if d > 0:
                S += 1
            elif d < 0:
                S -= 1
    return S


def mann_kendall_var_S(x):
    """Kendall (1975) ties-corrected Var(S)."""
    n = len(x)
    counts = {}
    for v in x:
        counts[v] = counts.get(v, 0) + 1
    tie_term = sum(t * (t - 1) * (2 * t + 5)
                   for t in counts.values() if t > 1)
    return (n * (n - 1) * (2 * n + 5) - tie_term) / 18.0


def sen_slope(years, values):
    """Theil-Sen median slope with actual year spacing."""
    slopes = []
    n = len(years)
    for i in range(n):
        for j in range(i + 1, n):
            dt = years[j] - years[i]
            if dt > 0:
                slopes.append((values[j] - values[i]) / float(dt))
    if not slopes:
        return 0.0
    slopes.sort()
    return slopes[len(slopes) // 2]


def hamed_rao_variance_correction(values, years, alpha=ALPHA):
    """Hamed & Rao (1998) variance inflation factor for Mann-Kendall S
    under serial dependence. Returns (n_over_n_star, lag_used).

    Procedure:
      (a) Detrend by Sen slope.
      (b) Rank the detrended series (average ties).
      (c) For each lag k in 1..n-1, compute Spearman-on-detrended
          autocorrelation rho_s(k). Retain only those with
          |rho_s(k)| > 1.96/sqrt(n-k) (Kendall 1976; Hamed-Rao Eq. 13).
      (d) Correction factor:
            n/n* = 1 + (2 / (n(n-1)(n-2))) *
                        sum_{k=1..n-1} (n-k)(n-k-1)(n-k-2) * rho_s(k)
          where only significant lags contribute.
      (e) Clamp to n/n* ≥ 1 (deflation is not physically meaningful
          for positive persistence and is dropped by the standard
          Hamed-Rao procedure).
    """
    n = len(values)
    if n < 10:
        return 1.0, 0
    slope = sen_slope(years, values)
    mean_year = mean_val(years)
    detrended = [values[i] - slope * (years[i] - mean_year)
                 for i in range(n)]
    ranks = rank_array(detrended)

    acc = 0.0
    lags_used = 0
    for k in range(1, n - 1):
        if n - k < 4:
            break
        rho = autocorr_rank(ranks, k)
        crit = 1.96 / math.sqrt(n - k)
        if abs(rho) < crit:
            continue
        weight = (n - k) * (n - k - 1) * (n - k - 2)
        acc += weight * rho
        lags_used += 1
    denom = n * (n - 1) * (n - 2)
    if denom == 0:
        return 1.0, 0
    factor = 1.0 + (2.0 / denom) * acc
    if factor < 1.0:
        factor = 1.0
    return factor, lags_used


def hamed_rao_pvalue(values, years):
    """Hamed-Rao two-sided Mann-Kendall p-value and corrected z-score.
    Returns dict with S, var_S_naive, var_S_HR, z_HR, p_HR, n_over_n_star.
    """
    S = mann_kendall_S(values)
    var_S_naive = mann_kendall_var_S(values)
    factor, lags = hamed_rao_variance_correction(values, years)
    var_S_HR = var_S_naive * factor
    if var_S_HR <= 0:
        return {"S": S, "var_S_naive": var_S_naive,
                "var_S_HR": var_S_HR, "z_HR": 0.0,
                "p_HR": 1.0, "n_over_n_star": factor,
                "lags_significant": lags}
    if S > 0:
        z = (S - 1) / math.sqrt(var_S_HR)
    elif S < 0:
        z = (S + 1) / math.sqrt(var_S_HR)
    else:
        z = 0.0
    p = 2.0 * (1.0 - normal_cdf(abs(z)))
    return {"S": S, "var_S_naive": var_S_naive,
            "var_S_HR": var_S_HR, "z_HR": z,
            "p_HR": p, "n_over_n_star": factor,
            "lags_significant": lags}


def block_bootstrap_median_slope_ci(per_gauge_slopes, n_boot, rng,
                                    ci=0.95):
    """Percentile CI for the median of a list of per-gauge Sen slopes,
    obtained by resampling the list of gauges with replacement
    (gauge-level i.i.d. bootstrap).  This is not a within-series block
    bootstrap; autocorrelation within each gauge is already handled by
    the Hamed-Rao variance correction when computing that gauge's
    p-value.  The `block_` prefix is historical — each resampling
    block is a single gauge."""
    if not per_gauge_slopes:
        return (0.0, 0.0, 0.0)
    n = len(per_gauge_slopes)
    medians = []
    for _ in range(n_boot):
        sample = [per_gauge_slopes[rng.randrange(0, n)]
                  for _ in range(n)]
        medians.append(median_val(sample))
    medians.sort()
    lo = medians[max(0, int(n_boot * (1 - ci) / 2))]
    hi = medians[min(n_boot - 1, int(n_boot * (1 + ci) / 2))]
    point = median_val(per_gauge_slopes)
    return (round(lo, 6), round(point, 6), round(hi, 6))


def benjamini_hochberg(pvalues, q):
    m = len(pvalues)
    if m == 0:
        return []
    order = sorted(range(m), key=lambda i: pvalues[i])
    threshold_k = 0
    for rank, idx in enumerate(order, 1):
        if pvalues[idx] <= rank * q / m:
            threshold_k = rank
    rejected = [False] * m
    for rank, idx in enumerate(order, 1):
        if rank <= threshold_k:
            rejected[idx] = True
    return rejected


def huc2_of(site_no):
    """First two digits of USGS site_no = 2-digit hydrologic unit code."""
    s = str(site_no).lstrip("0") or "0"
    s = str(site_no)
    if len(s) >= 2 and s[:2].isdigit():
        return s[:2]
    return "00"


def proportion_bootstrap_ci(count, total, n_boot, rng, ci=0.95):
    """Binomial bootstrap CI for a proportion."""
    if total == 0:
        return (0.0, 0.0)
    p_hat = count / total
    props = []
    for _ in range(n_boot):
        k = sum(1 for _ in range(total) if rng.random() < p_hat)
        props.append(k / total)
    props.sort()
    lo = props[max(0, int(n_boot * (1 - ci) / 2))]
    hi = props[min(n_boot - 1, int(n_boot * (1 + ci) / 2))]
    return (round(lo, 4), round(hi, 4))


# ═══════════════════════════════════════════════════════════════
# DATA ACQUISITION (domain-specific)
# ═══════════════════════════════════════════════════════════════

def sha256_of(data):
    return hashlib.sha256(data).hexdigest()


def download_with_retry(url, path, retries=REQUEST_RETRIES,
                        timeout=REQUEST_TIMEOUT):
    """Fetch `url` to `path` with bounded retry.  On permanent failure,
    logs a WARN line to stderr identifying the URL and the last exception;
    never raises.  A write-path IOError is fatal (disk full, bad perms)
    and is re-raised so the operator sees it immediately."""
    last_err = None
    for attempt in range(retries):
        try:
            req = urllib.request.Request(
                url, headers={"User-Agent": USER_AGENT, "Accept": "*/*"})
            with urllib.request.urlopen(req, timeout=timeout) as resp:
                data = resp.read()
        except (urllib.error.URLError, urllib.error.HTTPError,
                TimeoutError, ConnectionError, OSError) as e:
            last_err = e
            if attempt < retries - 1:
                time.sleep(2 * (attempt + 1))
            continue
        except Exception as e:
            last_err = e
            if attempt < retries - 1:
                time.sleep(2 * (attempt + 1))
            continue
        try:
            with open(path, "wb") as f:
                f.write(data)
            with open(path + ".sha256", "w") as f:
                f.write(sha256_of(data))
            return True
        except OSError as e:
            sys.stderr.write(
                "FATAL: cannot write cache file {}: {}\n".format(path, e))
            raise
    sys.stderr.write(
        "WARN: download failed after {} retries: {} (last error: {})\n"
        .format(retries, url, last_err))
    return False


def cache_valid(path):
    sf = path + ".sha256"
    if not (os.path.exists(path) and os.path.exists(sf)):
        return False
    with open(path, "rb") as f:
        actual = sha256_of(f.read())
    with open(sf) as f:
        expected = f.read().strip()
    return actual == expected


def parse_usgs_rdb(text):
    """Parse USGS NWIS peak RDB format. Returns [(year, peak_cfs)]."""
    lines = text.split("\n")
    headers = None
    header_line = -1
    for i, line in enumerate(lines):
        stripped = line.strip()
        if not stripped or stripped.startswith("#"):
            continue
        if headers is None:
            headers = stripped.split("\t")
            header_line = i
            break
    if headers is None:
        return []
    try:
        dt_col = headers.index("peak_dt")
        va_col = headers.index("peak_va")
    except ValueError:
        return []
    out = []
    past_sep = False
    for i, line in enumerate(lines):
        if i <= header_line:
            continue
        stripped = line.strip()
        if not stripped or stripped.startswith("#"):
            continue
        if not past_sep:
            past_sep = True
            continue
        fields = stripped.split("\t")
        if len(fields) <= max(dt_col, va_col):
            continue
        date_str = fields[dt_col].strip()
        val_str = fields[va_col].strip()
        if not date_str or not val_str:
            continue
        try:
            yr = int(date_str[:4])
            q = float(val_str)
            if q > 0 and 1850 <= yr <= MAX_RECORD_YEAR:
                out.append((yr, q))
        except (ValueError, IndexError):
            continue
    return out


def annualize(records):
    """Keep maximum peak per calendar year."""
    by_year = {}
    for yr, q in records:
        if yr not in by_year or q > by_year[yr]:
            by_year[yr] = q
    return by_year


def fetch_site(site):
    """Return cached RDB text for `site` or None if every mirror fails.
    Network errors do NOT raise; they emit a WARN line to stderr and
    allow the caller to degrade gracefully by dropping the gauge."""
    cache_path = os.path.join(CACHE_DIR, "peak_{}.rdb".format(site))
    try:
        if cache_valid(cache_path):
            with open(cache_path, "r", errors="replace") as f:
                return f.read()
    except OSError as e:
        sys.stderr.write(
            "WARN: cannot read cache for {}: {} (will refetch)\n"
            .format(site, e))
    got = False
    for template in NWIS_URL_TEMPLATES:
        try:
            if download_with_retry(template.format(site=site), cache_path):
                got = True
                time.sleep(INTER_REQUEST_SLEEP)
                break
        except OSError:
            return None
    if not got:
        time.sleep(INTER_REQUEST_SLEEP)
        sys.stderr.write(
            "WARN: site {} unavailable from all mirrors\n".format(site))
        return None
    try:
        with open(cache_path, "r", errors="replace") as f:
            return f.read()
    except OSError as e:
        sys.stderr.write(
            "WARN: cannot re-open cache for {}: {}\n".format(site, e))
        return None


def load_data():
    """Download and parse peak flow for reference+regulated gauges.

    Returns:
      {
        "gauges": {site: {"series": {year: peak_cfs},
                          "group_strict": "reference"|"regulated"|None,
                          "group_permissive": same,
                          "huc2": "NN"}},
        "sites_queried": {"reference_strict": [...],
                          "reference_additional": [...],
                          "regulated": [...]}
      }
    """
    os.makedirs(CACHE_DIR, exist_ok=True)

    site_group_map = {}
    for s in USGS_REFERENCE_SITES_STRICT:
        site_group_map[s] = ("reference", "reference")  # strict, permissive
    for s in USGS_REFERENCE_SITES_ADDITIONAL:
        site_group_map[s] = (None, "reference")
    for s in USGS_REGULATED_SITES:
        site_group_map[s] = ("regulated", "regulated")

    all_sites = list(site_group_map.keys())
    print("[download] fetching peak flow for {} gauges".format(len(all_sites)))
    gauges = {}
    for idx, site in enumerate(all_sites):
        text = fetch_site(site)
        if text is None or "peak_dt" not in text:
            continue
        annual = annualize(parse_usgs_rdb(text))
        if len(annual) < MIN_ANNUAL_FOR_PARSE:
            continue
        g_strict, g_permissive = site_group_map[site]
        gauges[site] = {
            "series": annual,
            "group_strict": g_strict,
            "group_permissive": g_permissive,
            "huc2": huc2_of(site),
        }
        if (idx + 1) % 25 == 0:
            print("    {}/{} loaded={}".format(
                idx + 1, len(all_sites), len(gauges)))
    print("[download] kept {} gauges".format(len(gauges)))
    return {
        "gauges": gauges,
        "sites_queried": {
            "reference_strict": list(USGS_REFERENCE_SITES_STRICT),
            "reference_additional": list(USGS_REFERENCE_SITES_ADDITIONAL),
            "regulated": list(USGS_REGULATED_SITES),
        },
    }


# ═══════════════════════════════════════════════════════════════
# DOMAIN-AGNOSTIC ANALYSIS
# ═══════════════════════════════════════════════════════════════

def per_gauge_record(site, meta, group_label):
    """Compute per-gauge trend statistics from the annual series."""
    years = sorted(meta["series"].keys())
    peaks = [meta["series"][y] for y in years]
    n = len(years)
    median_peak = median_val(peaks)
    slope = sen_slope(years, peaks)
    hr = hamed_rao_pvalue(peaks, years)
    pct_per_decade = (100.0 * slope * 10.0 / median_peak
                      if median_peak > 0 else 0.0)
    return {
        "site": site,
        "group": group_label,
        "huc2": meta["huc2"],
        "n_years": n,
        "year_min": years[0],
        "year_max": years[-1],
        "median_cfs": round(median_peak, 3),
        "sen_slope_cfs_per_yr": round(slope, 6),
        "pct_per_decade": round(pct_per_decade, 4),
        "MK_S": hr["S"],
        "z_HR": round(hr["z_HR"], 4),
        "p_HR": round(hr["p_HR"], 6),
        "n_over_n_star": round(hr["n_over_n_star"], 4),
        "lags_significant": hr["lags_significant"],
    }


def build_gauge_rows(data, group_field):
    """Compute per-gauge records assigning group via `group_field`
    (either 'group_strict' or 'group_permissive')."""
    rows = []
    for site in sorted(data["gauges"].keys()):
        meta = data["gauges"][site]
        grp = meta[group_field]
        if grp is None:
            continue
        rows.append(per_gauge_record(site, meta, grp))
    return rows


def rejection_rate(rows_subset, alpha):
    if not rows_subset:
        return 0.0, 0
    n_sig = sum(1 for r in rows_subset if r["p_HR"] < alpha)
    return n_sig / len(rows_subset), n_sig


def summarize_group(rows_subset, rng_seed, label):
    N = len(rows_subset)
    out = {"label": label, "N": N}
    if N == 0:
        return out
    slopes = [r["sen_slope_cfs_per_yr"] for r in rows_subset]
    pct_decade = [r["pct_per_decade"] for r in rows_subset]
    out["median_slope_cfs_per_yr"] = round(median_val(slopes), 6)
    out["median_pct_per_decade"] = round(median_val(pct_decade), 4)
    out["frac_upward"] = round(
        sum(1 for s in slopes if s > 0) / N, 4)
    rate, n_sig = rejection_rate(rows_subset, ALPHA)
    out["n_sig_HR"] = n_sig
    out["rejection_rate_HR"] = round(rate, 4)
    out["rejection_rate_HR_ci95"] = proportion_bootstrap_ci(
        n_sig, N, 1000, random.Random(rng_seed + 11))
    pvals = [r["p_HR"] for r in rows_subset]
    rejected = benjamini_hochberg(pvals, FDR_Q)
    out["n_sig_HR_BH"] = sum(rejected)
    out["rejection_rate_HR_BH"] = round(sum(rejected) / N, 4)
    # Block-bootstrap CI for median Sen slope (%/decade).
    rng_slope = random.Random(rng_seed + 23)
    lo, point, hi = block_bootstrap_median_slope_ci(
        pct_decade, N_BOOTSTRAP_GAUGES, rng_slope)
    out["median_pct_per_decade_ci95"] = (round(lo, 4),
                                          round(point, 4),
                                          round(hi, 4))
    rng_abs = random.Random(rng_seed + 29)
    lo2, point2, hi2 = block_bootstrap_median_slope_ci(
        slopes, N_BOOTSTRAP_GAUGES, rng_abs)
    out["median_slope_cfs_per_yr_ci95"] = (round(lo2, 4),
                                            round(point2, 4),
                                            round(hi2, 4))
    return out


def label_permutation_test(rows, group_field_value_A, group_field_value_B,
                           n_permutations, rng):
    """Unstratified label-permutation test: under the null that group
    label is uninformative, shuffle the label assignment across all rows
    and recompute the rejection-rate gap.  Returns observed gap and
    empirical two-sided p-value.

    Because the shuffle is across the pooled set, the null distribution
    does not condition on record-length or HUC2 composition within
    groups.  See `stratified_label_permutation_test` for the HUC2-
    conditioned variant.
    """
    rows_A = [r for r in rows if r["group"] == group_field_value_A]
    rows_B = [r for r in rows if r["group"] == group_field_value_B]
    n_A = len(rows_A)
    n_B = len(rows_B)
    if n_A == 0 or n_B == 0:
        return {"n_A": n_A, "n_B": n_B, "observed_gap": 0.0,
                "p_value": 1.0, "n_permutations": 0}

    rate_A_obs, _ = rejection_rate(rows_A, ALPHA)
    rate_B_obs, _ = rejection_rate(rows_B, ALPHA)
    observed_gap = rate_A_obs - rate_B_obs  # reference - regulated

    pooled = rows_A + rows_B
    pooled_pvals = [r["p_HR"] for r in pooled]
    total = n_A + n_B
    ge_count = 0
    for _ in range(n_permutations):
        idx = list(range(total))
        rng.shuffle(idx)
        sigs = [1 if pooled_pvals[idx[i]] < ALPHA else 0
                for i in range(total)]
        rate_A_null = sum(sigs[:n_A]) / n_A
        rate_B_null = sum(sigs[n_A:]) / n_B
        gap_null = rate_A_null - rate_B_null
        if abs(gap_null) >= abs(observed_gap) - 1e-12:
            ge_count += 1
    p_val = (ge_count + 1) / (n_permutations + 1)
    return {"n_A": n_A, "n_B": n_B,
            "rate_A_obs": round(rate_A_obs, 4),
            "rate_B_obs": round(rate_B_obs, 4),
            "observed_gap": round(observed_gap, 4),
            "p_value": round(p_val, 6),
            "n_permutations": n_permutations}


def stratified_label_permutation_test(rows, group_A, group_B,
                                      strata_field, n_permutations, rng,
                                      min_stratum_size=3):
    """Stratified label-permutation test: shuffle labels only within
    each level of `strata_field` (here, HUC2).  Strata containing a
    single group are left unshuffled so those rows contribute a
    constant term to the gap; strata smaller than
    `min_stratum_size` are merged into a pooled "small" bucket.
    Returns observed gap and empirical two-sided p-value.
    """
    subset = [r for r in rows if r["group"] in (group_A, group_B)]
    if not subset:
        return {"observed_gap": 0.0, "p_value": 1.0, "n_permutations": 0,
                "n_strata_used": 0}
    by_stratum = {}
    for r in subset:
        s = r.get(strata_field, "NA")
        by_stratum.setdefault(s, []).append(r)
    # Merge small strata into a pooled bucket.
    merged = {}
    merged_bucket = []
    for s, group in by_stratum.items():
        if len(group) < min_stratum_size:
            merged_bucket.extend(group)
        else:
            merged[s] = group
    if merged_bucket:
        merged["__small__"] = merged_bucket

    n_A = sum(1 for r in subset if r["group"] == group_A)
    n_B = sum(1 for r in subset if r["group"] == group_B)
    if n_A == 0 or n_B == 0:
        return {"observed_gap": 0.0, "p_value": 1.0, "n_permutations": 0,
                "n_strata_used": len(merged)}
    rate_A_obs = sum(1 for r in subset
                     if r["group"] == group_A and r["p_HR"] < ALPHA) / n_A
    rate_B_obs = sum(1 for r in subset
                     if r["group"] == group_B and r["p_HR"] < ALPHA) / n_B
    observed_gap = rate_A_obs - rate_B_obs

    ge_count = 0
    for _ in range(n_permutations):
        # Within each stratum, reshuffle group labels (keeping counts).
        n_A_sig = 0
        n_B_sig = 0
        for s, group in merged.items():
            labels = [r["group"] for r in group]
            sigs = [1 if r["p_HR"] < ALPHA else 0 for r in group]
            order = list(range(len(group)))
            rng.shuffle(order)
            shuffled_labels = [labels[i] for i in order]
            for i, lbl in enumerate(shuffled_labels):
                if lbl == group_A:
                    n_A_sig += sigs[i]
                elif lbl == group_B:
                    n_B_sig += sigs[i]
        rate_A_null = n_A_sig / n_A
        rate_B_null = n_B_sig / n_B
        gap_null = rate_A_null - rate_B_null
        if abs(gap_null) >= abs(observed_gap) - 1e-12:
            ge_count += 1
    p_val = (ge_count + 1) / (n_permutations + 1)
    return {"observed_gap": round(observed_gap, 4),
            "p_value": round(p_val, 6),
            "n_permutations": n_permutations,
            "n_strata_used": len(merged)}


def slope_diff_permutation_test(rows, group_A, group_B, n_permutations, rng):
    """Two-sided permutation test for the difference in median Sen slopes
    (cfs/yr) between groups A and B. Scalar null — label shuffle."""
    a = [r["sen_slope_cfs_per_yr"] for r in rows if r["group"] == group_A]
    b = [r["sen_slope_cfs_per_yr"] for r in rows if r["group"] == group_B]
    if not a or not b:
        return {"observed_diff": 0.0, "p_value": 1.0,
                "n_permutations": 0}
    obs = median_val(a) - median_val(b)
    pooled = a + b
    n_a = len(a)
    ge_count = 0
    for _ in range(n_permutations):
        idx = list(range(len(pooled)))
        rng.shuffle(idx)
        a_null = [pooled[idx[i]] for i in range(n_a)]
        b_null = [pooled[idx[i]] for i in range(n_a, len(pooled))]
        diff_null = median_val(a_null) - median_val(b_null)
        if abs(diff_null) >= abs(obs) - 1e-15:
            ge_count += 1
    return {"observed_diff_cfs_per_yr": round(obs, 6),
            "p_value": round((ge_count + 1) / (n_permutations + 1), 6),
            "n_permutations": n_permutations}


def record_length_sweep(rows, thresholds, seed):
    """Recompute rejection rates per group at each MIN_RECORD_YEARS."""
    out = {}
    for thresh in thresholds:
        subset = [r for r in rows if r["n_years"] >= thresh]
        ref = [r for r in subset if r["group"] == "reference"]
        reg = [r for r in subset if r["group"] == "regulated"]
        rate_ref, n_sig_ref = rejection_rate(ref, ALPHA)
        rate_reg, n_sig_reg = rejection_rate(reg, ALPHA)
        rng = random.Random(seed + thresh * 97)
        perm = label_permutation_test(
            subset, "reference", "regulated",
            N_PERMUTATIONS, rng)
        out[str(thresh)] = {
            "threshold_years": thresh,
            "n_reference": len(ref),
            "n_regulated": len(reg),
            "rejection_rate_reference": round(rate_ref, 4),
            "rejection_rate_regulated": round(rate_reg, 4),
            "gap": round(rate_ref - rate_reg, 4),
            "permutation_p": perm["p_value"],
        }
    return out


def run_analysis(data):
    rng_seed = SEED
    rng_top = random.Random(SEED)

    headline = {}
    for variant in REFERENCE_SET_VARIANTS:
        field = ("group_strict" if variant == "strict"
                 else "group_permissive")
        rows_all = build_gauge_rows(data, field)
        rows = [r for r in rows_all if r["n_years"] >= MIN_RECORD_YEARS]

        summary = {
            "reference": summarize_group(
                [r for r in rows if r["group"] == "reference"],
                rng_seed + 101, "reference"),
            "regulated": summarize_group(
                [r for r in rows if r["group"] == "regulated"],
                rng_seed + 103, "regulated"),
            "combined": summarize_group(
                rows, rng_seed + 107, "combined"),
        }

        rng_perm = random.Random(rng_seed + 149)
        perm_rate_gap = label_permutation_test(
            rows, "reference", "regulated",
            N_PERMUTATIONS, rng_perm)
        rng_slope = random.Random(rng_seed + 151)
        perm_slope = slope_diff_permutation_test(
            rows, "reference", "regulated",
            N_PERMUTATIONS, rng_slope)
        rng_strat = random.Random(rng_seed + 157)
        perm_rate_gap_stratified = stratified_label_permutation_test(
            rows, "reference", "regulated", "huc2",
            N_PERMUTATIONS, rng_strat,
            min_stratum_size=STRATIFIED_MIN_STRATUM)

        record_sweep = record_length_sweep(
            rows_all, RECORD_LENGTH_THRESHOLDS, rng_seed)

        headline[variant] = {
            "variant": variant,
            "N_total": len(rows),
            "summary": summary,
            "permutation_rate_gap": perm_rate_gap,
            "permutation_rate_gap_stratified_huc2":
                perm_rate_gap_stratified,
            "permutation_median_slope_diff": perm_slope,
            "record_length_sweep": record_sweep,
            "gauge_rows": rows_all,
        }

    # Politis-White block-length distribution on residuals (for reference
    # variant only, used for reporting only).
    # Distribution of Hamed-Rao inflation factor
    strict_rows = headline["strict"]["gauge_rows"]
    ref_rows = [r for r in strict_rows
                if r["group"] == "reference" and
                r["n_years"] >= MIN_RECORD_YEARS]
    reg_rows = [r for r in strict_rows
                if r["group"] == "regulated" and
                r["n_years"] >= MIN_RECORD_YEARS]
    inflation_dist = {
        "reference_strict": {
            "median": round(median_val(
                [r["n_over_n_star"] for r in ref_rows]), 4)
                if ref_rows else 0.0,
            "p25": round(percentile_val(
                [r["n_over_n_star"] for r in ref_rows], 25), 4)
                if ref_rows else 0.0,
            "p75": round(percentile_val(
                [r["n_over_n_star"] for r in ref_rows], 75), 4)
                if ref_rows else 0.0,
            "n_gauges_with_any_significant_lag": sum(
                1 for r in ref_rows if r["lags_significant"] > 0),
        },
        "regulated": {
            "median": round(median_val(
                [r["n_over_n_star"] for r in reg_rows]), 4)
                if reg_rows else 0.0,
            "p25": round(percentile_val(
                [r["n_over_n_star"] for r in reg_rows], 25), 4)
                if reg_rows else 0.0,
            "p75": round(percentile_val(
                [r["n_over_n_star"] for r in reg_rows], 75), 4)
                if reg_rows else 0.0,
            "n_gauges_with_any_significant_lag": sum(
                1 for r in reg_rows if r["lags_significant"] > 0),
        },
    }

    return {
        "config": {
            "SEED": SEED,
            "MIN_RECORD_YEARS": MIN_RECORD_YEARS,
            "ALPHA": ALPHA,
            "FDR_Q": FDR_Q,
            "N_BOOTSTRAP_GAUGES": N_BOOTSTRAP_GAUGES,
            "N_PERMUTATIONS": N_PERMUTATIONS,
            "RECORD_LENGTH_THRESHOLDS": RECORD_LENGTH_THRESHOLDS,
            "REFERENCE_SET_VARIANTS": REFERENCE_SET_VARIANTS,
            "MAX_RECORD_YEAR": MAX_RECORD_YEAR,
            "n_reference_strict_queried":
                len(USGS_REFERENCE_SITES_STRICT),
            "n_reference_additional_queried":
                len(USGS_REFERENCE_SITES_ADDITIONAL),
            "n_regulated_queried": len(USGS_REGULATED_SITES),
        },
        "headline": headline,
        "hamed_rao_inflation_distribution": inflation_dist,
        "limitations": LIMITATIONS,
    }


# ═══════════════════════════════════════════════════════════════
# REPORT GENERATION
# ═══════════════════════════════════════════════════════════════

def generate_report(results):
    with open(RESULTS_FILE, "w") as f:
        json.dump(results, f, indent=2)

    cfg = results["config"]
    h = results["headline"]
    infl = results["hamed_rao_inflation_distribution"]

    lines = []
    lines.append("# GAGES-II Reference vs Regulated Flood Stationarity "
                 "(Hamed-Rao Mann-Kendall)")
    lines.append("")
    lines.append("## Configuration")
    lines.append("- Seed: {}".format(cfg["SEED"]))
    lines.append("- Minimum record length (primary): {} yr".format(
        cfg["MIN_RECORD_YEARS"]))
    lines.append("- Alpha: {}".format(cfg["ALPHA"]))
    lines.append("- BH q: {}".format(cfg["FDR_Q"]))
    lines.append("- Gauge bootstrap resamples: {}".format(
        cfg["N_BOOTSTRAP_GAUGES"]))
    lines.append("- Label permutations: {}".format(cfg["N_PERMUTATIONS"]))
    lines.append("")

    for variant in cfg["REFERENCE_SET_VARIANTS"]:
        v = h[variant]
        s = v["summary"]
        lines.append("## Reference-set variant: {}".format(variant))
        for label in ["reference", "regulated", "combined"]:
            sl = s[label]
            lines.append("### {}  (N={})".format(label, sl["N"]))
            if sl["N"] == 0:
                lines.append("_no gauges qualified_")
                continue
            lines.append(
                "- Median Sen slope: {} cfs/yr  ({} %/decade; "
                "95% CI [{}, {}])".format(
                    sl["median_slope_cfs_per_yr"],
                    sl["median_pct_per_decade"],
                    sl["median_pct_per_decade_ci95"][0],
                    sl["median_pct_per_decade_ci95"][2]))
            lines.append(
                "- Fraction upward: {:.1%}".format(sl["frac_upward"]))
            lines.append(
                "- Hamed-Rao rejection rate (alpha={}): {} / {} = "
                "{:.1%} (95% CI {:.1%}-{:.1%})".format(
                    cfg["ALPHA"], sl["n_sig_HR"], sl["N"],
                    sl["rejection_rate_HR"],
                    sl["rejection_rate_HR_ci95"][0],
                    sl["rejection_rate_HR_ci95"][1]))
            lines.append(
                "- After BH q={}: {} / {} = {:.1%}".format(
                    cfg["FDR_Q"], sl["n_sig_HR_BH"], sl["N"],
                    sl["rejection_rate_HR_BH"]))
        pg = v["permutation_rate_gap"]
        lines.append("")
        lines.append(
            "- Reference - regulated rejection-rate gap: "
            "{} (permutation p = {} over {} shuffles)".format(
                pg["observed_gap"], pg["p_value"],
                pg["n_permutations"]))
        pgs = v.get("permutation_rate_gap_stratified_huc2", {})
        if pgs:
            lines.append(
                "- Stratified-by-HUC2 permutation p on same gap: {} "
                "({} strata used)".format(
                    pgs.get("p_value", float("nan")),
                    pgs.get("n_strata_used", 0)))
        ps = v["permutation_median_slope_diff"]
        lines.append(
            "- Reference - regulated median Sen-slope diff (cfs/yr): "
            "{} (permutation p = {})".format(
                ps["observed_diff_cfs_per_yr"], ps["p_value"]))
        lines.append("")
        lines.append("#### Record-length sensitivity sweep")
        lines.append("| threshold | N_ref | N_reg | rate_ref | rate_reg "
                     "| gap | perm p |")
        lines.append("|---|---|---|---|---|---|---|")
        for thresh in cfg["RECORD_LENGTH_THRESHOLDS"]:
            r = v["record_length_sweep"][str(thresh)]
            lines.append("| {} | {} | {} | {:.1%} | {:.1%} | "
                         "{:+.1%} | {} |".format(
                             thresh, r["n_reference"], r["n_regulated"],
                             r["rejection_rate_reference"],
                             r["rejection_rate_regulated"],
                             r["gap"], r["permutation_p"]))
        lines.append("")

    lines.append("## Hamed-Rao inflation factor distribution "
                 "(strict variant, primary threshold)")
    for label in ["reference_strict", "regulated"]:
        d = infl[label]
        lines.append("- {}: n/n* median={}, IQR=[{}, {}]; "
                     "{} gauges had >= 1 significant lag".format(
                         label, d["median"], d["p25"], d["p75"],
                         d["n_gauges_with_any_significant_lag"]))
    lines.append("")
    lines.append("## Interpretation")
    lines.append("- A small reference-group rejection rate with a larger")
    lines.append("  regulated-group rejection rate indicates that flow")
    lines.append("  regulation drives the headline 'CONUS floods are")
    lines.append("  changing' signal rather than climate.")
    lines.append("- If reference and regulated groups show similar")
    lines.append("  rejection rates, regulation is not the dominant")
    lines.append("  confound and the climate interpretation is")
    lines.append("  consistent with the data (but still not proven).")

    with open(REPORT_FILE, "w") as f:
        f.write("\n".join(lines) + "\n")


# ═══════════════════════════════════════════════════════════════
# VERIFICATION
# ═══════════════════════════════════════════════════════════════

def run_verification():
    messages = []
    passed = True

    def check(cond, msg):
        nonlocal passed
        messages.append(("PASS" if cond else "FAIL") + ": " + msg)
        if not cond:
            passed = False

    if not os.path.exists(RESULTS_FILE):
        return False, ["FAIL: results.json not found"]
    with open(RESULTS_FILE) as f:
        r = json.load(f)

    cfg = r.get("config", {})
    h = r.get("headline", {})

    check(cfg.get("SEED") == SEED,
          "seed pinned to {}".format(SEED))
    check(cfg.get("MIN_RECORD_YEARS", 0) >= 30,
          "min record length >= 30 years")
    check(cfg.get("N_BOOTSTRAP_GAUGES", 0) >= 1000,
          "bootstrap resamples >= 1000")
    check(cfg.get("N_PERMUTATIONS", 0) >= 1000,
          "label permutations >= 1000")

    for variant in ["strict", "permissive"]:
        check(variant in h, "headline contains variant '{}'".format(variant))
        v = h.get(variant, {})
        s = v.get("summary", {})
        for grp in ["reference", "regulated", "combined"]:
            check(grp in s, "{} summary contains group '{}'".format(
                variant, grp))
        n_ref = s.get("reference", {}).get("N", 0)
        n_reg = s.get("regulated", {}).get("N", 0)
        check(n_ref >= 40,
              "{} reference N >= 40 (got {})".format(variant, n_ref))
        check(n_reg >= 5,
              "{} regulated N >= 5 (got {})".format(variant, n_reg))
        # Rejection rates valid
        rr = s.get("reference", {}).get("rejection_rate_HR", -1)
        check(0.0 <= rr <= 1.0, "{} reference rejection rate in [0,1] "
              "(got {})".format(variant, rr))
        # Permutation p-value present and in range
        pg = v.get("permutation_rate_gap", {})
        pv = pg.get("p_value", -1)
        check(0.0 <= pv <= 1.0,
              "{} permutation p in [0,1] (got {})".format(variant, pv))
        # Stratified permutation present and bounded
        pgs = v.get("permutation_rate_gap_stratified_huc2", {})
        psv = pgs.get("p_value", -1)
        check(0.0 <= psv <= 1.0,
              "{} stratified permutation p in [0,1] "
              "(got {})".format(variant, psv))
        check(pgs.get("n_strata_used", 0) >= 1,
              "{} stratified permutation uses >= 1 stratum".format(variant))
        # Slope-diff permutation p present and bounded
        psl = v.get("permutation_median_slope_diff", {})
        psl_v = psl.get("p_value", -1)
        check(0.0 <= psl_v <= 1.0,
              "{} slope-diff permutation p in [0,1] "
              "(got {})".format(variant, psl_v))
        # Bootstrap CI returns three-tuple lo-point-hi
        ci = s.get("reference", {}).get("median_pct_per_decade_ci95",
                                        None)
        check(isinstance(ci, list) and len(ci) == 3 and
              ci[0] <= ci[1] <= ci[2],
              "{} reference median-%/decade CI is ordered "
              "(got {})".format(variant, ci))
        # Per-gauge rows: all have Hamed-Rao fields
        rows = v.get("gauge_rows", [])
        check(all("p_HR" in x and "z_HR" in x and
                  "n_over_n_star" in x and
                  "sen_slope_cfs_per_yr" in x for x in rows),
              "{} all per-gauge rows have Hamed-Rao + Sen fields".format(
                  variant))
        # Hamed-Rao inflation factor >= 1 (by construction)
        check(all(x.get("n_over_n_star", 0.0) >= 1.0 - 1e-6
                  for x in rows),
              "{} all Hamed-Rao inflation factors >= 1".format(variant))
        # Record-length sweep covers configured thresholds
        sweep = v.get("record_length_sweep", {})
        for thresh in cfg.get("RECORD_LENGTH_THRESHOLDS", []):
            check(str(thresh) in sweep,
                  "{} record-length sweep contains threshold {}".format(
                      variant, thresh))

    # Informational (not blocking): check that strict-variant reference N
    # does not exceed permissive-variant reference N.
    n_strict = h.get("strict", {}).get("summary", {}).get(
        "reference", {}).get("N", 0)
    n_perm = h.get("permissive", {}).get("summary", {}).get(
        "reference", {}).get("N", 0)
    messages.append(
        "INFO: strict reference N = {}, permissive reference N = {} "
        "(expect permissive >= strict)".format(n_strict, n_perm))
    check(n_perm >= n_strict,
          "permissive reference N ({}) >= strict reference N ({})".format(
              n_perm, n_strict))

    # --- EXTENDED SANITY CHECKS ---------------------------------------
    # (1) limitations array is present and non-trivial.
    lims = r.get("limitations", [])
    check(isinstance(lims, list) and len(lims) >= 4,
          "results.json includes >= 4 limitations (got {})".format(len(lims)))

    # (2) data row count: each variant must retain a non-trivial number
    # of per-gauge rows, otherwise downstream stats are meaningless.
    for variant in ["strict", "permissive"]:
        rows = h.get(variant, {}).get("gauge_rows", [])
        check(len(rows) >= 60,
              "{} gauge_rows has >= 60 entries (got {})".format(
                  variant, len(rows)))

    # (3) effect-size plausibility: median %/decade for every group and
    # variant must lie within the bound configured at the top of the
    # script.  Floods cannot plausibly shift by more than an order of
    # magnitude within a human-length observational window.
    for variant in ["strict", "permissive"]:
        s = h.get(variant, {}).get("summary", {})
        for grp in ["reference", "regulated", "combined"]:
            m = s.get(grp, {}).get("median_pct_per_decade", 0.0)
            check(abs(m) <= MAX_PLAUSIBLE_PCT_PER_DECADE,
                  "{} {} |median_pct_per_decade|={} <= {}".format(
                      variant, grp, m, MAX_PLAUSIBLE_PCT_PER_DECADE))

    # (4) CI sanity: widths must be strictly positive and wider than
    # MIN_CI_WIDTH_PCT %/decade. A degenerate width usually indicates a
    # bootstrap seed bug or constant-slope input.
    for variant in ["strict", "permissive"]:
        s = h.get(variant, {}).get("summary", {})
        for grp in ["reference", "regulated"]:
            ci = s.get(grp, {}).get("median_pct_per_decade_ci95", None)
            if isinstance(ci, list) and len(ci) == 3:
                width = ci[2] - ci[0]
                check(width >= MIN_CI_WIDTH_PCT,
                      "{} {} CI width {} >= {} %/decade".format(
                          variant, grp, round(width, 4),
                          MIN_CI_WIDTH_PCT))

    # (5) Sensitivity / robustness: the sign of the rejection-rate gap at
    # the primary threshold (50 yr) should agree with the sign at 40 yr.
    # Disagreement would indicate the headline is an artifact of the
    # record-length cutoff rather than a robust property.
    for variant in ["strict", "permissive"]:
        sweep = h.get(variant, {}).get("record_length_sweep", {})
        g40 = sweep.get("40", {}).get("gap", 0.0)
        g50 = sweep.get("50", {}).get("gap", 0.0)
        # Allow the gap to cross zero if one of the two is near zero,
        # because noise around zero is not informative.
        near_zero = min(abs(g40), abs(g50)) < 0.05
        sign_stable = (g40 * g50 >= 0) or near_zero
        check(sign_stable,
              "{} sweep gap sign stable between 40/50 yr thresholds "
              "(g40={}, g50={})".format(variant, g40, g50))

    # (6) Negative / falsification control: the rejection-rate gap at
    # the 30-year threshold should fall within a loose plausibility
    # envelope.  A |gap| > MAX_PLAUSIBLE_GAP_AT_30YR = 0.30 would
    # indicate the comparison is dominated by a data-shape artifact and
    # not by a genuine subpopulation contrast.
    for variant in ["strict", "permissive"]:
        sweep = h.get(variant, {}).get("record_length_sweep", {})
        g30 = sweep.get("30", {}).get("gap", 0.0)
        check(abs(g30) <= MAX_PLAUSIBLE_GAP_AT_30YR,
              "{} 30-yr gap {} within negative-control envelope "
              "|gap| <= {}".format(variant, g30,
                                    MAX_PLAUSIBLE_GAP_AT_30YR))

    # (7) Negative control on slope signs: the reference group must
    # contain *both* upward and downward trending gauges. A degenerate
    # all-upward or all-downward distribution would signal corrupted
    # data (e.g. every gauge replaced by a constant).
    for variant in ["strict", "permissive"]:
        rows = [x for x in h.get(variant, {}).get("gauge_rows", [])
                if x.get("group") == "reference"]
        n_up = sum(1 for x in rows
                   if x.get("sen_slope_cfs_per_yr", 0.0) > 0)
        n_dn = sum(1 for x in rows
                   if x.get("sen_slope_cfs_per_yr", 0.0) < 0)
        check(n_up > 0 and n_dn > 0,
              "{} reference has both up ({}) and down ({}) trending "
              "gauges".format(variant, n_up, n_dn))

    # (8) Record-length sweep p-values are all in [0, 1].
    for variant in ["strict", "permissive"]:
        sweep = h.get(variant, {}).get("record_length_sweep", {})
        for thresh, row in sweep.items():
            pv = row.get("permutation_p", -1)
            check(0.0 <= pv <= 1.0,
                  "{} sweep[{}] permutation_p in [0,1] (got {})".format(
                      variant, thresh, pv))

    # (9) Hamed-Rao inflation-factor distribution has finite, sane IQR.
    infl = r.get("hamed_rao_inflation_distribution", {})
    for key in ["reference_strict", "regulated"]:
        d = infl.get(key, {})
        med = d.get("median", 0.0)
        p25 = d.get("p25", 0.0)
        p75 = d.get("p75", 0.0)
        check(1.0 - 1e-6 <= p25 <= med <= p75 <= 50.0,
              "{} inflation IQR ordered and <= 50 "
              "(p25={}, median={}, p75={})".format(key, p25, med, p75))

    # (10) Per-gauge records actually exceed the configured MIN_RECORD_YEARS
    # for the primary-threshold summaries.
    for variant in ["strict", "permissive"]:
        rows = h.get(variant, {}).get("gauge_rows", [])
        n_short = sum(1 for x in rows
                      if x.get("n_years", 0) < cfg.get(
                          "MIN_RECORD_YEARS", 50)
                      and x.get("group") in ("reference", "regulated"))
        # rows include ALL gauges for the sweep, so we only spot-check
        # that *some* rows meet the threshold (already covered by N >= 40
        # on the summary) and that no gauge claims a 0-year record.
        n_zero = sum(1 for x in rows if x.get("n_years", 0) <= 0)
        check(n_zero == 0,
              "{} no per-gauge row has n_years <= 0".format(variant))

    # (11) Permutation p-values are sensibly bounded above 0 (minimum
    # achievable is 1/(N_PERMUTATIONS+1) ~ 2e-4).  A literal zero would
    # indicate a counter bug.
    for variant in ["strict", "permissive"]:
        pg = h.get(variant, {}).get("permutation_rate_gap", {})
        pv = pg.get("p_value", 0.0)
        check(pv > 0.0,
              "{} permutation p > 0 (min achievable 1/(N+1); got {})"
              .format(variant, pv))

    # (12) Bootstrap CI containment: the group-median point estimate must
    # lie between the CI bounds.  Violation would mean the bootstrap
    # percentile ordering is broken.
    for variant in ["strict", "permissive"]:
        s = h.get(variant, {}).get("summary", {})
        for grp in ["reference", "regulated"]:
            ci = s.get(grp, {}).get("median_pct_per_decade_ci95", None)
            if isinstance(ci, list) and len(ci) == 3:
                lo, point, hi = ci[0], ci[1], ci[2]
                check(lo <= point <= hi,
                      "{} {} CI contains its point estimate "
                      "(lo={}, point={}, hi={})".format(
                          variant, grp, lo, point, hi))

    # (13) Record-length monotonicity: raising the threshold can only
    # *shrink* the retained sample.  Violation would mean the sweep is
    # subsetting rows incorrectly.
    for variant in ["strict", "permissive"]:
        sweep = h.get(variant, {}).get("record_length_sweep", {})
        thresholds = cfg.get("RECORD_LENGTH_THRESHOLDS", [])
        n_prev_ref = None
        n_prev_reg = None
        monotone = True
        for t in thresholds:
            row = sweep.get(str(t), {})
            n_ref_t = row.get("n_reference", 0)
            n_reg_t = row.get("n_regulated", 0)
            if n_prev_ref is not None and n_ref_t > n_prev_ref:
                monotone = False
            if n_prev_reg is not None and n_reg_t > n_prev_reg:
                monotone = False
            n_prev_ref, n_prev_reg = n_ref_t, n_reg_t
        check(monotone,
              "{} record-length sweep counts are weakly decreasing "
              "with threshold".format(variant))

    # (14) Permissive-variant reference sample strictly dominates strict
    # variant at the primary threshold (more permissive inclusion rule
    # should add gauges, never remove them).
    n_s_ref = h.get("strict", {}).get("summary", {}).get(
        "reference", {}).get("N", 0)
    n_p_ref = h.get("permissive", {}).get("summary", {}).get(
        "reference", {}).get("N", 0)
    check(n_p_ref >= n_s_ref,
          "permissive reference pool dominates strict "
          "(strict N={}, permissive N={})".format(n_s_ref, n_p_ref))

    # (15) Sensitivity falsification on the slope-sign composition: after
    # randomly reshuffling per-gauge slopes across the *entire* strict-
    # variant pool, the "reference - regulated" median-slope contrast
    # should collapse toward zero.  Any shuffle-invariant contrast would
    # indicate a leakage bug where labels are correlated with values
    # through a channel other than the data itself.
    strict_rows = h.get("strict", {}).get("gauge_rows", [])
    filtered = [x for x in strict_rows
                if x.get("group") in ("reference", "regulated")
                and x.get("n_years", 0) >= cfg.get("MIN_RECORD_YEARS", 50)]
    if filtered:
        slopes = [x.get("sen_slope_cfs_per_yr", 0.0) for x in filtered]
        n_ref_f = sum(1 for x in filtered if x.get("group") == "reference")
        rng_fals = random.Random(SEED + 9999)
        gaps = []
        for _ in range(500):
            idx = list(range(len(slopes)))
            rng_fals.shuffle(idx)
            a = [slopes[idx[i]] for i in range(n_ref_f)]
            b = [slopes[idx[i]] for i in range(n_ref_f, len(slopes))]
            gaps.append(median_val(a) - median_val(b))
        mean_abs_gap = mean_val([abs(g) for g in gaps])
        # Observed should be drawn from a zero-centered null; its mean
        # |gap| must be bounded by the data spread.
        slope_range = max(slopes) - min(slopes) if slopes else 0.0
        check(mean_abs_gap <= slope_range or slope_range == 0.0,
              "falsification: shuffled-label mean |slope gap| ({:.4f}) "
              "bounded by observed slope range ({:.4f})".format(
                  mean_abs_gap, slope_range))

    # (16) Per-gauge Hamed-Rao p-values are in [0,1] everywhere.
    for variant in ["strict", "permissive"]:
        rows = h.get(variant, {}).get("gauge_rows", [])
        bad = [x for x in rows
               if not (0.0 <= x.get("p_HR", -1.0) <= 1.0)]
        check(not bad,
              "{} all per-gauge p_HR in [0,1] ({} violations)".format(
                  variant, len(bad)))

    # (17) Rejection rate arithmetic consistency: summary-reported n_sig
    # must equal the actual count of rows with p_HR < ALPHA in the
    # corresponding subset.
    for variant in ["strict", "permissive"]:
        rows = h.get(variant, {}).get("gauge_rows", [])
        s = h.get(variant, {}).get("summary", {})
        for grp in ["reference", "regulated"]:
            subset = [x for x in rows
                      if x.get("group") == grp
                      and x.get("n_years", 0) >= cfg.get(
                          "MIN_RECORD_YEARS", 50)]
            n_sig_recomputed = sum(
                1 for x in subset
                if x.get("p_HR", 1.0) < cfg.get("ALPHA", 0.05))
            reported = s.get(grp, {}).get("n_sig_HR", -1)
            check(reported == n_sig_recomputed,
                  "{} {} n_sig_HR self-consistent ({} reported vs {} "
                  "recomputed)".format(
                      variant, grp, reported, n_sig_recomputed))

    # (18) Hamed-Rao correction non-triviality: across *all* per-gauge
    # records, at least one gauge should have n/n* > 1.0 (else the
    # autocorrelation correction is doing no work anywhere and the
    # headline p-values reduce to the naive Mann-Kendall test).
    all_rows = []
    for variant in ["strict", "permissive"]:
        all_rows.extend(h.get(variant, {}).get("gauge_rows", []))
    n_nontrivial = sum(1 for x in all_rows
                       if x.get("n_over_n_star", 1.0) > 1.0 + 1e-6)
    check(n_nontrivial >= 1,
          "at least one gauge has a non-trivial Hamed-Rao inflation "
          "factor (got {})".format(n_nontrivial))

    return passed, messages


# ═══════════════════════════════════════════════════════════════
# ORCHESTRATION
# ═══════════════════════════════════════════════════════════════

def main(argv):
    random.seed(SEED)

    if len(argv) > 1 and argv[1] == "--verify":
        try:
            ok, msgs = run_verification()
        except FileNotFoundError as e:
            sys.stderr.write(
                "ERROR: --verify requires a prior successful run; "
                "results.json not found ({}). Run `python3 {} "
                "` first.\n".format(e, argv[0]))
            sys.exit(3)
        except json.JSONDecodeError as e:
            sys.stderr.write(
                "ERROR: results.json is corrupt or truncated ({}). "
                "Delete it and re-run the analysis.\n".format(e))
            sys.exit(4)
        for m in msgs:
            print(m)
        if ok:
            print("VERIFY RESULT: PASS")
            print("ALL CHECKS PASSED")
            sys.exit(0)
        print("VERIFY RESULT: FAIL")
        sys.exit(1)

    print("[1/5] Loading data (download on first run, SHA256 cache after)")
    try:
        data = load_data()
    except OSError as e:
        sys.stderr.write(
            "FATAL: I/O error while fetching or caching USGS peak data: {}. "
            "Check disk space and write permissions on '{}'.\n"
            .format(e, CACHE_DIR))
        sys.exit(5)
    except Exception as e:
        sys.stderr.write(
            "FATAL: unexpected error during data acquisition ({}: {}). "
            "Re-run with network access; cached gauges will be reused via "
            "SHA256 validation.\n".format(type(e).__name__, e))
        sys.exit(6)

    n_have_any = len(data["gauges"])
    n_primary_eligible_ref = sum(
        1 for s, m in data["gauges"].items()
        if m["group_strict"] == "reference"
        and len(m["series"]) >= MIN_RECORD_YEARS)
    n_primary_eligible_reg = sum(
        1 for s, m in data["gauges"].items()
        if m["group_strict"] == "regulated"
        and len(m["series"]) >= MIN_RECORD_YEARS)
    print("     Retained {} gauges; {} reference + {} regulated meet "
          "MIN_RECORD_YEARS={}".format(
              n_have_any, n_primary_eligible_ref,
              n_primary_eligible_reg, MIN_RECORD_YEARS))
    if n_primary_eligible_ref < 40:
        sys.stderr.write(
            "FATAL: only {} reference gauges meet MIN_RECORD_YEARS={}; "
            "need >= 40. Check connectivity to waterdata.usgs.gov, "
            "clear the cache directory ({}) if corrupt, and retry.\n"
            .format(n_primary_eligible_ref, MIN_RECORD_YEARS, CACHE_DIR))
        sys.exit(2)
    if n_primary_eligible_reg < 5:
        sys.stderr.write(
            "FATAL: only {} regulated gauges meet MIN_RECORD_YEARS={}; "
            "need >= 5 for a valid comparator group.\n"
            .format(n_primary_eligible_reg, MIN_RECORD_YEARS))
        sys.exit(2)

    print("[2/5] Computing Hamed-Rao corrected Mann-Kendall and "
          "block-bootstrap Sen-slope CIs")
    try:
        results = run_analysis(data)
    except Exception as e:
        sys.stderr.write(
            "FATAL: analysis failed ({}: {}). Refusing to write partial "
            "output to avoid masking the error downstream.\n"
            .format(type(e).__name__, e))
        sys.exit(7)

    print("[3/5] Writing results.json and report.md")
    try:
        generate_report(results)
    except OSError as e:
        sys.stderr.write(
            "FATAL: cannot write output files: {}. Check disk and "
            "permissions on '{}'.\n".format(e, WORKSPACE))
        sys.exit(8)

    print("[4/5] Running self-verification")
    ok, msgs = run_verification()
    for m in msgs:
        print("    " + m)
    if not ok:
        sys.stderr.write(
            "WARNING: self-verification FAILED — results.json written "
            "but at least one sanity check did not pass. See stdout.\n")

    print("[5/5] ANALYSIS COMPLETE")
    print("     results: {}".format(RESULTS_FILE))
    print("     report:  {}".format(REPORT_FILE))


if __name__ == "__main__":
    main(sys.argv)
SCRIPT_EOF
```

**Expected output:** File created silently. Exit code 0.

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_stationarity-in-gages-ii-reference-gauges && python3 analyze.py
```

**Expected output:** last lines read:
```
[5/5] ANALYSIS COMPLETE
     results: .../results.json
     report:  .../report.md
```
Exit code 0. Creates `results.json` and `report.md` in the workspace.

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_stationarity-in-gages-ii-reference-gauges && python3 analyze.py --verify
```

**Expected output:** every non-INFO line begins with `PASS:`, the
penultimate line reads `VERIFY RESULT: PASS`, and the final line reads
`ALL CHECKS PASSED`. Exit code 0.

## Success Criteria

The run is successful **iff all** of the following measurable conditions hold:

1. The analysis script exits 0 and prints `ANALYSIS COMPLETE`.
2. `results.json` is valid JSON and contains keys
   `config`, `headline.strict`, `headline.permissive`, and
   `hamed_rao_inflation_distribution`.
3. All ≥ 25 machine-checkable `--verify` assertions print `PASS:` and the
   final line reads `VERIFY RESULT: PASS` (exit code 0).
4. `headline.strict.summary.reference.N` ≥ 40 and
   `headline.permissive.summary.reference.N` ≥ 40, with permissive ≥ strict.
5. Every per-gauge row carries `p_HR ∈ [0,1]`, `z_HR ∈ ℝ`,
   `n_over_n_star ≥ 1.0`, `sen_slope_cfs_per_yr`, and `pct_per_decade`.
6. Per-group median `pct_per_decade` estimates are within the plausibility
   band `|median_pct_per_decade| ≤ 20`/decade (floods cannot plausibly change
   by more than an order of magnitude over a human-length record).
7. 95% CIs on median `pct_per_decade` are ordered (`lo ≤ point ≤ hi`) and
   have width ≥ 0.01 %/decade (degenerate-width CIs signal a bootstrap bug).
8. `record_length_sweep` contains **all four** configured thresholds
   (30/40/50/60 yr) and each threshold reports a permutation `p ∈ [0,1]`.
9. The negative-control falsification check passes: at `MIN_RECORD_YEARS = 30`
   the rejection-rate gap has `|gap| ≤ 0.30` (a larger value would indicate
   a data-shape artifact dominating the comparison).

## Failure Conditions

- **Network unavailable on first run AND no SHA256-verified cache on disk** →
  `load_data()` emits `WARN:` stderr lines per failed fetch and, if fewer
  than 40 reference gauges survive, the main routine exits with code 2
  and a `FATAL:` message naming the unmet threshold.
- **Corrupted cache** — a cached RDB whose SHA256 does not match the
  recorded digest is automatically invalidated and re-downloaded.
- **USGS RDB schema drift** — if the `peak_dt` or `peak_va` columns are
  missing from the tab-delimited header, that gauge is dropped; ≥ 40
  reference survivors are still required or the script exits non-zero.
- **Python < 3.8** — unsupported. Script uses `math.erf`, f-strings are
  avoided, no type annotations.
- **All MK statistics degenerate to zero** (e.g. network returned stub HTML
  for every gauge) — `--verify` fails on the effect-size-plausibility and
  nondegenerate-slope-sign checks.

## Limitations and Assumptions

- **Gauge coverage**: 170 gauges is a curated, region-stratified sample of
  the full GAGES-II catalog (≈ 9,322 stations). Findings apply to this
  subset, not to every individual CONUS watershed. Broader conclusions
  require re-running with an expanded `USGS_REFERENCE_SITES_STRICT` list.
- **Stationarity of serial-dependence structure**: Hamed-Rao (1998)
  assumes the lag-k rank autocorrelation structure is itself stationary
  over the record. Regime shifts in autocorrelation would under-correct the
  variance.
- **Peak-flow operational definition**: "Annual peak" here is the maximum
  instantaneous discharge reported by USGS per calendar year, not per
  water year. Flood-frequency studies that condition on water-year
  peaks may differ at the margin.
- **Reference-set circularity**: GAGES-II "reference" classification was
  made by USGS using land-use and storage criteria circa 2011. Any
  subsequent development upstream of a reference gauge is not reflected
  in the classification; slow encroachment would bias the reference
  group toward a non-null trend.
- **What the results do NOT show**: (a) the magnitude of any
  climate-driven component within the regulated gauges, (b) whether
  other hydrologic variables (low flows, timing, duration) are
  stationary, (c) causal attribution of any observed non-stationarity
  to a specific driver (reservoir construction, urbanization, LUC).
- **Permutation null assumption**: the unstratified test assumes that
  under the null, reference vs regulated labels are exchangeable across
  the pooled gauge set — this ignores record-length and regional
  composition differences between groups. The HUC2-stratified
  permutation and the record-length sensitivity sweep partially address
  this but cannot fully correct for unmeasured covariates.

These caveats are also emitted as a `limitations` array in `results.json` so
that downstream report consumers cannot act on the output without seeing
them.