watershed.js
172 lines · in natura-src-1.0.0.zip as natura/watershed.js · a module · its entry in the module reference · its test: watershed-test.js
// watershed.js -- Natura GROUNDWATER + DRAINAGE (the subsurface layer the water CYCLE (hydro.js) leaves out).
// hydro.js moves water through its phases above ground; ocean.js routes coastal rivers; but nowhere does the
// land itself have a DRAINAGE NETWORK carved by where the water must go, nor a hidden WATER TABLE beneath it.
// Here the terrain is read once and solved: every hollow is filled so water never dead-ends (priority-flood),
// flow is routed downhill (D8) and ACCUMULATED, so the branching tree of streams-to-rivers-to-sea emerges as
// the memory of where the water went. Beneath it sits the AQUIFER -- a smoothed water table that RISES with
// recharge (couple it to hydro's moisture or weather's rain) and FALLS with warmth/drought; where it breaks
// the surface, SPRINGS flow; where it sinks away, the ground goes dry. So drought and flush here are not a
// dialled knob but an emergent consequence of the coupled climate. Deterministic, node-testable. UMD.
(function (root, factory) {
if (typeof module !== 'undefined' && module.exports) module.exports = factory();
else root.Watershed = factory();
})(typeof self !== 'undefined' ? self : this, function () {
'use strict';
const clamp = (v, a, b) => v < a ? a : v > b ? b : v;
// config: { N, elevation(x,y)|elevationFrom, seaLevel, riverFrac, tableDepth, tableBlur, rechargeFrom,
// recharge, rechargeK, drawdownK, leakK, temperature|weather, maxDtDays, name }
function makeWatershed(config) {
config = config || {};
const N = config.N || 48;
const SEA = config.seaLevel != null ? config.seaLevel : 0.34;
const riverFrac = config.riverFrac != null ? config.riverFrac : 0.006; // a cell is a stream once this frac of the land drains through it
const tableDepth = config.tableDepth != null ? config.tableDepth : 0.02; // how far the table sits below its regional potential
const tableBlur = config.tableBlur != null ? config.tableBlur : 6; // groundwater is a SMOOTHED replica of the topography
const rechargeK = config.rechargeK != null ? config.rechargeK : 0.010; // rain lifts the table
const drawdownK = config.drawdownK != null ? config.drawdownK : 0.008; // warmth/drought draws it down
const leakK = config.leakK != null ? config.leakK : 0.03; // the table relaxes toward its balance
const maxOff = config.maxOff != null ? config.maxOff : 0.06; // clamp on the seasonal excursion
const maxDd = config.maxDtDays != null ? config.maxDtDays : 4;
const rechargeFrom = config.rechargeFrom || null; // a system whose moisture/liquid field feeds the table
const weatherName = config.weather || null; // a system whose temp field drives drawdown
const resolveEveryDays = config.resolveEveryDays != null ? config.resolveEveryDays : 50; // re-carve the drainage as terrain erodes (slow)
const rechargeCfg = config.recharge != null ? config.recharge : 0.5; // fallback recharge 0..1 when uncoupled
const tempCfg = config.temperature != null ? config.temperature : 15; // fallback temp when uncoupled
const elevFrom = config.elevationFrom || null; // read elevation from a coupled terrain field's key
const elevFn = config.elevation || ((x, y) => 0.08 + 0.82 * (0.5 - 0.5 * y)); // default: a plane tilted N-high to S-sea
const gwDiffuse = config.gwDiffuse != null ? config.gwDiffuse : 0.12; // groundwater spreads sideways through the ground
const drawFrom = config.draw || null; // OPTIONAL: a system exposing a `draw` field -- HUMAN extraction (wells) pumping the table DOWN
const drawK = config.drawK != null ? config.drawK : 1.0; // extraction gain: how hard a unit of `draw` lowers the local table
const gx = (i) => (i / (N - 1)) * 2 - 1;
const scal = (c, x, y) => typeof c === 'function' ? c(x, y) : c;
// ---- a binary min-heap (for priority-flood) ----
function Heap(cap) { const h = new Float64Array(cap + 1), idx = new Int32Array(cap + 1); let n = 0;
return {
get size() { return n; },
push(hv, ii) { n++; let c = n; h[c] = hv; idx[c] = ii; while (c > 1) { const p = c >> 1; if (h[p] <= h[c]) break; const th = h[p]; h[p] = h[c]; h[c] = th; const ti = idx[p]; idx[p] = idx[c]; idx[c] = ti; c = p; } },
pop() { const ri = idx[1], rh = h[1]; h[1] = h[n]; idx[1] = idx[n]; n--; let c = 1; while (true) { let l = c * 2, r = l + 1, s = c; if (l <= n && h[l] < h[s]) s = l; if (r <= n && h[r] < h[s]) s = r; if (s === c) break; const th = h[s]; h[s] = h[c]; h[c] = th; const ti = idx[s]; idx[s] = idx[c]; idx[c] = ti; c = s; } return { i: ri, h: rh }; }
};
}
function boxBlur(src, w, hh, r) { const tmp = new Float32Array(w * hh), out = new Float32Array(w * hh), win = r * 2 + 1;
for (let y = 0; y < hh; y++) { const row = y * w; let s = 0; for (let k = -r; k <= r; k++) s += src[row + clamp(k, 0, w - 1)]; for (let x = 0; x < w; x++) { tmp[row + x] = s / win; s += src[row + Math.min(w - 1, x + r + 1)] - src[row + Math.max(0, x - r)]; } }
for (let x = 0; x < w; x++) { let s = 0; for (let k = -r; k <= r; k++) s += tmp[clamp(k, 0, hh - 1) * w + x]; for (let y = 0; y < hh; y++) { out[y * w + x] = s / win; s += tmp[Math.min(hh - 1, y + r + 1) * w + x] - tmp[Math.max(0, y - r) * w + x]; } }
return out; }
const NB = [[1, 0], [-1, 0], [0, 1], [0, -1], [1, 1], [1, -1], [-1, 1], [-1, -1]];
// solve the static hydrology of a heightfield: filled surface, D8 down-pointers, flow accumulation, lakes
function solve(elev) {
const M = N * N, filled = new Float32Array(M), closed = new Uint8Array(M), EPS = 1e-5, heap = Heap(M);
for (let x = 0; x < N; x++) { const a = x, b = (N - 1) * N + x; closed[a] = closed[b] = 1; heap.push(elev[a], a); heap.push(elev[b], b); }
for (let y = 0; y < N; y++) { const a = y * N, b = y * N + N - 1; if (!closed[a]) { closed[a] = 1; heap.push(elev[a], a); } if (!closed[b]) { closed[b] = 1; heap.push(elev[b], b); } }
while (heap.size > 0) { const c = heap.pop(), cx = c.i % N, cy = (c.i / N) | 0; filled[c.i] = c.h;
for (const d of NB) { const nx = cx + d[0], ny = cy + d[1]; if (nx < 0 || ny < 0 || nx >= N || ny >= N) continue; const ni = ny * N + nx; if (closed[ni]) continue; closed[ni] = 1; const nh = Math.max(elev[ni], c.h + EPS); filled[ni] = nh; heap.push(nh, ni); } }
const down = new Int32Array(M).fill(-1);
for (let y = 0; y < N; y++) for (let x = 0; x < N; x++) { const i = y * N + x; if (elev[i] < SEA) continue; let best = filled[i], bi = -1;
for (const d of NB) { const nx = x + d[0], ny = y + d[1]; if (nx < 0 || ny < 0 || nx >= N || ny >= N) continue; const ni = ny * N + nx; if (filled[ni] < best) { best = filled[ni]; bi = ni; } } down[i] = bi; }
const acc = new Float32Array(M).fill(1);
const order = []; for (let i = 0; i < M; i++) if (elev[i] >= SEA) order.push(i);
order.sort((a, b) => filled[b] - filled[a]);
for (const i of order) { const dn = down[i]; if (dn >= 0) acc[dn] += acc[i]; }
const lake = new Uint8Array(M); for (let i = 0; i < M; i++) if (elev[i] >= SEA && filled[i] - elev[i] > 0.004) lake[i] = 1;
const base = boxBlur(filled, N, N, tableBlur); for (let i = 0; i < M; i++) base[i] -= tableDepth; // regional groundwater potential
let land = 0; for (let i = 0; i < M; i++) if (elev[i] >= SEA) land++;
return { filled, down, acc, lake, base, land };
}
function meanField(f, elev) { if (!f) return null; let s = 0, n = 0; for (let i = 0; i < f.length; i++) { if (elev && elev[i] < SEA) continue; s += f[i]; n++; } return n ? s / n : 0; }
return {
name: config.name || 'watershed',
domain: 'physical',
scale: 'landscape',
params: config,
init() {
const M = N * N, elev = new Float32Array(M);
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) elev[j * N + i] = elevFn(gx(i), gx(j));
const sol = solve(elev);
const riverThresh = Math.max(3, sol.land * riverFrac);
return { N, SEA, elev, riverThresh, adj: new Float32Array(M), ...sol }; // adj = the SPATIAL table excursion per cell
},
step(st, ctx, dt) {
const dd = clamp(dt / 86400, 0, maxDd); if (dd <= 0) return st;
// COUPLE TO A FORMING TERRAIN: rivers are the geology's own agents (terrain.js erodes the land) --
// so re-derive the drainage + aquifer from terrain's EVOLVING elevation, occasionally (erosion is slow).
if (elevFrom && ctx && ctx.fields && ctx.fields[elevFrom] && ctx.fields[elevFrom].elevation) {
st._sinceSolve = (st._sinceSolve || 0) + dd;
if (!st._solvedCoupled || st._sinceSolve >= resolveEveryDays) {
const te = ctx.fields[elevFrom].elevation, M = N * N;
for (let i = 0; i < M; i++) st.elev[i] = te[i];
const sol = solve(st.elev);
st.filled = sol.filled; st.down = sol.down; st.acc = sol.acc; st.lake = sol.lake; st.base = sol.base; st.land = sol.land;
st.riverThresh = Math.max(3, sol.land * riverFrac);
st._solvedCoupled = true; st._sinceSolve = 0;
}
}
// the table is a SPATIAL field now: it rises where rain + rivers recharge it, falls where it is warm --
// so springs appear + dry in the right PLACES, not everywhere at once.
const M = N * N;
const cf = rechargeFrom && ctx && ctx.fields && ctx.fields[rechargeFrom] ? ctx.fields[rechargeFrom] : null;
const moist = cf && cf.moisture ? cf.moisture : null, liqF = cf && cf.liquid ? cf.liquid : null;
const wf = weatherName && ctx && ctx.fields && ctx.fields[weatherName] ? ctx.fields[weatherName] : null;
const tempF = wf && wf.temp ? wf.temp : null;
const drawF = drawFrom && ctx && ctx.fields && ctx.fields[drawFrom] && ctx.fields[drawFrom].draw ? ctx.fields[drawFrom].draw : null; // HUMAN extraction
const adj = st.adj;
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) { const k = j * N + i; if (st.elev[k] < SEA) { adj[k] = 0; continue; }
let rc = moist ? clamp(moist[k], 0, 1.5) : liqF ? clamp(liqF[k] / 0.5, 0, 1.5) : scal(rechargeCfg, gx(i), gx(j));
if (st.acc[k] >= st.riverThresh) rc += 0.4; // a stream RECHARGES the aquifer it runs over (an influent reach)
const temp = tempF ? tempF[k] : scal(tempCfg, gx(i), gx(j));
const warmN = clamp((temp - 2) / 25, 0, 1.4);
const draw = drawF ? drawK * drawF[k] : 0; // wells pump the aquifer DOWN, locally -- the human hand on the water
adj[k] += (rechargeK * rc - drawdownK * warmN - draw - leakK * adj[k]) * dd;
}
// groundwater SPREADS sideways through the ground (a wet region lifts its neighbours; capped for stability)
const GD = Math.min(0.2, gwDiffuse * dd);
if (GD > 0) { const d = st._adjd || (st._adjd = new Float32Array(M)); d.fill(0);
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) { const k = j * N + i; if (st.elev[k] < SEA) continue;
if (i < N - 1 && st.elev[k + 1] >= SEA) { const f = GD * (adj[k] - adj[k + 1]); d[k] -= f; d[k + 1] += f; }
if (j < N - 1 && st.elev[k + N] >= SEA) { const f = GD * (adj[k] - adj[k + N]); d[k] -= f; d[k + N] += f; } }
for (let k = 0; k < M; k++) adj[k] += d[k]; }
for (let k = 0; k < M; k++) { if (adj[k] > maxOff) adj[k] = maxOff; else if (adj[k] < -maxOff) adj[k] = -maxOff; }
return st;
},
// OUT: the drainage (flowAccum, riverMask) + the aquifer (waterTable, depthToWater) + springs. Other systems
// read `depthToWater` (vegetation thrives where it is shallow) and `springMask` (life gathers at a spring).
fields(st) {
const M = st.N * st.N;
const wt = st._wt || (st._wt = new Float32Array(M));
const depth = st._depth || (st._depth = new Float32Array(M));
const river = st._river || (st._river = new Uint8Array(M));
const spring = st._spring || (st._spring = new Uint8Array(M));
let springCount = 0, riverCount = 0, maxAcc = 0, adjSum = 0, adjN = 0;
for (let i = 0; i < M; i++) {
if (st.elev[i] < SEA) { wt[i] = SEA; depth[i] = 0; river[i] = 0; spring[i] = 0; continue; }
wt[i] = st.base[i] + st.adj[i]; // the SPATIAL water table
depth[i] = st.filled[i] - wt[i];
river[i] = st.acc[i] >= st.riverThresh ? 1 : 0; if (river[i]) riverCount++;
if (st.acc[i] > maxAcc) maxAcc = st.acc[i];
spring[i] = (!river[i] && !st.lake[i] && depth[i] <= 0.004) ? 1 : 0; if (spring[i]) springCount++;
adjSum += st.adj[i]; adjN++;
}
return { N: st.N, elev: st.elev, filled: st.filled, flowAccum: st.acc, riverMask: river, waterTable: wt,
depthToWater: depth, springMask: spring, lake: st.lake, adj: st.adj, offset: adjN ? adjSum / adjN : 0,
springCount, riverCount, maxAcc, land: st.land };
},
// EVENTS: the coupled climate tips the country between flush and drought (on the mean table excursion)
events(st) {
let s = 0, n = 0; for (let i = 0; i < st.N * st.N; i++) { if (st.elev[i] < SEA) continue; s += st.adj[i]; n++; }
const m = n ? s / n : 0, out = [];
if (m <= -maxOff * 0.8) out.push({ type: 'drought', offset: +m.toFixed(4) });
else if (m >= maxOff * 0.8) out.push({ type: 'flush', offset: +m.toFixed(4) });
return out;
}
};
}
return { makeWatershed };
});