# Methodology — Petrel Karst & Depressions Layer v0.1

## 1. Depression depth

For each 10-degree tile, the bare-earth DEM is read **with a 256-pixel halo
(~23 km) on every side**, depressions are filled to their lowest pour point,
and the result is cropped back to the tile core.

Filling uses the Wang & Liu (2006) priority-flood algorithm via `pyflwdir`
(Deltares, MIT). Every cell is raised to the elevation of the lowest point over
which water could escape to the grid boundary; the shipped value is

    depression_depth = filled_surface - original_surface

so a pixel that already drains reads exactly 0, and a pixel at the floor of a
basin reads the height of its rim above it, measured at the *pour point*, not
at the highest surrounding ground.

### Why the halo, and why 256 pixels specifically

A depression clipped by a tile edge is filled to the edge rather than to its
real pour point, and reports the wrong depth. The halo width was **measured,
not chosen**, by filling the same terrain under increasing windows and
comparing against a 512-pixel reference across six karst regions (Florida,
Dinaric, Causses, South China, Veneto, Mammoth Cave):

| halo | worst region error (mean / p99) |
|---|---|
| 64 px (~5.8 km) | South China cone karst: 3.76 m / 52.4 m |
| 128 px (~11.5 km) | South China cone karst: 1.72 m / 41.5 m |
| 256 px (~23 km) | 0.000 m in every region tested |

Florida alone converged at 64 px. Relying on it would have shipped 52-metre
errors in Chinese cone karst, where enclosed depressions are kilometres across
— one region is a sample, not a measurement.

The halo's effect on the product was then verified directly: the same tile
built at halo 256 against halo 0 differs in **0.00000** of interior pixels
(max 0.03 m) and by up to **41.58 m** in the edge frame. The largest difference
sits at 0 px from the tile edge and reads 41.58 m with the halo against 0.00 m
without — without it, a 41-metre depression on a tile boundary ships as
perfectly flat.

Unlike flow accumulation, this converges. A depression's pour point is local,
so a sufficiently wide window contains it; accumulation depends on the whole
upstream basin and cannot be fixed by any halo.

### Why bare earth

The input is GEDTM30 bare earth, not a surface model. Buildings and tree canopy
in a surface DEM create closed depressions that are not terrain — a courtyard,
a clearing in a forest — and those artefacts are indistinguishable from
dolines once filled.

## 2. Bedrock solubility

`karst_solubility` is an ordinal reclass of GLiM lithology, ordered by how
readily the rock dissolves: 0 insoluble, 1 possible (mixed sediments),
2 carbonate, 3 evaporite. Evaporites are kept separate from carbonates
deliberately — gypsum and halite dissolve orders of magnitude faster and
produce catastrophic collapse rather than slow doline development.

It is reclassed from the lithology tile on the **same grid as the DEM**, so the
two bands are co-registered by construction rather than by assumption.

## 3. What was deliberately NOT done

**No probability band, and no model.** See README.md and model_card.md: the
lawful training corpus is 103 verified positives, and the 4,232 unverified
Florida reports were measured to carry no usable signal (AUC 0.520 against
confirmed non-sinkholes). Fitting thirteen features to ~200 points would
produce a surface that looks like a product and means nothing.

**No solubility mask on the depth band.** The obvious inference — a closed
depression on carbonate is dissolution, on granite it is a kettle hole — was
built and then rejected on measurement: masking depth to solubility >= 2 moves
the Florida AUC from 0.649 [0.58, 0.71] to 0.520 [0.48, 0.56]. GLiM marks only
~20% of the labelled Florida points as soluble because the surficial geology
there is sand over carbonate. The two bands ship separately so a user can
combine them knowingly, in a region where they have evidence it helps.

**No aggregation of depression statistics.** Depression density, mean depth and
count were tested alongside depth and added nothing (0.49-0.59 individually,
and combining them scored below depth alone). Only depth ships.

## 4. Encoding

Depth is capped at 3000 m and encoded uint16 with a physical scale in the band
tags. The ceiling is set by measurement, not by the karst range: real karst
reaches 100 m (the Viñales mogotes in western Cuba measure 111 m in this build,
verified against the source DEM rather than assumed to be a spike), but fill
depth is *pour-point height above ground*, so large endorheic basins run to a
measured maximum of 3736 m. 3000 m clips one pixel globally and preserves the
validated Florida AUC; lower ceilings clip entire basin floors and higher ones
cost skill by adding quantisation ties. model_card.md carries the full table.
Solubility is uint8 with 255 as nodata. Land/ocean masking uses the canonical
GSHHG-f full-resolution shoreline.

## References

Wang, L. & Liu, H. (2006). An efficient method for identifying and filling
surface depressions in digital elevation models. *International Journal of
Geographical Information Science* 20(2), 193-213.