diff --git a/server/examples/bench_layer1.rs b/server/examples/bench_layer1.rs
new file mode 100644
index 000000000..5e4183d44
--- /dev/null
+++ b/server/examples/bench_layer1.rs
@@ -0,0 +1,102 @@
+//! Benchmark: per-body Layer-1 topography generation (#953).
+//!
+//! Times `run_layer1` (D8 drainage → terrain analysis → 7 feature tags →
+//! sub-biome classification) on a synthetic canonical 512×256 heightmap — the
+//! cost of "populating the top-level map for a body." Run with:
+//!
+//!   cargo run --release --example bench_layer1
+//!
+//! The committed systems.db carries no heightmaps yet (the importer isn't in
+//! regen-db), so this uses a multi-octave synthetic heightmap representative of
+//! a continental body. Heightmap *content* affects timing (river/attractor
+//! counts), so treat this as an order-of-magnitude figure against the D-208
+//! ~50ms target, not a per-body guarantee.
+
+use std::time::Instant;
+
+use settled_reach_server::atlas::heightmap::{BodyHeightmap, GRID_H, GRID_W};
+use settled_reach_server::atlas::layer1::run_layer1;
+
+/// Deterministic multi-octave sine "noise" → continents, coasts, drainage.
+fn synth_heightmap(w: u32, h: u32) -> Vec<f32> {
+    use std::f32::consts::{PI, TAU};
+    let mut data = vec![0f32; (w * h) as usize];
+    for r in 0..h {
+        for c in 0..w {
+            let x = c as f32 / w as f32 * TAU;
+            let y = r as f32 / h as f32 * PI;
+            let e = 0.5
+                + 0.25 * (x * 3.0).sin() * (y * 2.0).sin()
+                + 0.15 * (x * 7.0).cos() * (y * 5.0).sin()
+                + 0.08 * (x * 13.0).sin() * (y * 11.0).cos()
+                + 0.05 * (x * 23.0).cos() * (y * 19.0).sin();
+            data[(r * w + c) as usize] = e.clamp(0.0, 1.0);
+        }
+    }
+    data
+}
+
+fn main() {
+    let (w, h) = (GRID_W, GRID_H);
+    let hm = BodyHeightmap {
+        body_id: "Bench".into(),
+        width: w,
+        height: h,
+        data: synth_heightmap(w, h),
+        sea_level: 0.40,
+    };
+
+    // Warm-up + output characterization.
+    let o = run_layer1(&hm);
+    let mut by_type: std::collections::BTreeMap<String, usize> = std::collections::BTreeMap::new();
+    for a in &o.attractors {
+        *by_type.entry(format!("{:?}", a.attractor_type)).or_default() += 1;
+    }
+    println!("=== Layer-1 output for one 512×256 body ===");
+    println!(
+        "  river cells: {}   mouths: {}   confluences: {}   basins: {}",
+        o.river_network.river_cells.len(),
+        o.river_network.mouths.len(),
+        o.river_network.confluences.len(),
+        o.drainage_basins.len(),
+    );
+    println!("  attractors: {}  {:?}", o.attractors.len(), by_type);
+
+    // Phase breakdown (median of N) to locate the cost.
+    use settled_reach_server::atlas::drainage;
+    use settled_reach_server::atlas::features::{extract_attractors, TerrainAnalysis};
+    let n = 30;
+    let med = |mut v: Vec<f64>| {
+        v.sort_by(|a, b| a.partial_cmp(b).unwrap());
+        v[v.len() / 2]
+    };
+    let mut t_drain = vec![];
+    let mut t_terr = vec![];
+    let mut t_feat = vec![];
+    let mut t_total = vec![];
+    for _ in 0..n {
+        let t = Instant::now();
+        let dr = drainage::analyze(&hm.data, hm.width, hm.height, hm.sea_level);
+        t_drain.push(t.elapsed().as_secs_f64() * 1000.0);
+
+        let t = Instant::now();
+        let ta = TerrainAnalysis::analyze(&hm, &dr);
+        t_terr.push(t.elapsed().as_secs_f64() * 1000.0);
+
+        let t = Instant::now();
+        let raw = extract_attractors(&hm, &dr, &ta);
+        t_feat.push(t.elapsed().as_secs_f64() * 1000.0);
+        std::hint::black_box(&raw);
+
+        let t = Instant::now();
+        let out = run_layer1(&hm);
+        std::hint::black_box(&out);
+        t_total.push(t.elapsed().as_secs_f64() * 1000.0);
+    }
+    println!("\n=== phase breakdown (median of {n}, 512×256, release) ===");
+    println!("  D8 drainage (D-208):     {:.2}ms", med(t_drain));
+    println!("  terrain analysis:        {:.2}ms", med(t_terr));
+    println!("  feature extraction:      {:.2}ms", med(t_feat));
+    println!("  run_layer1 (total):      {:.2}ms", med(t_total));
+    println!("  D-208 target: ~50ms/body");
+}