cryosphere.js
145 lines · in natura-src-1.0.0.zip as natura/cryosphere.js · a module · its entry in the module reference · its test: cryosphere-test.js
// cryosphere.js -- Natura CRYOSPHERE (the ICE domain, SYSTEMS.md: the missing cold force under the physical
// stack). Land ice as a per-cell mass balance: snow ACCUMULATES where the air is below freezing and ice
// ABLATES (melts) where it is above, by a degree-year model. The whole thing is bound together by the
// ICE-ALBEDO FEEDBACK -- white ice reflects sunlight, so the more of the land is iced the COLDER the world
// gets, which grows more ice: a positive loop that makes the climate BISTABLE. That is the emergence this
// system exists to show -- an ICE AGE is not a smooth response to cooling but a TIPPING: the ice stays near
// zero as the world cools, then past a threshold it runs away to a glaciated planet, and to melt it back you
// must warm well PAST where it first froze (hysteresis). The ice also LOCKS UP water (the sea withdraws) and
// GRINDS the bedrock it sits on (glacial carving), so it couples out to the ocean and the terrain.
//
// Reads a deep-time temperature forcing (a scalar, or a fn of the clock) + optional STELLAR insolation (a
// fainter sun -> a colder world -> more ice) + an elevation sampler (mountains glaciate first). Deterministic,
// node-testable. UMD (node test / window.Cryosphere).
(function (root, factory) {
if (typeof module !== 'undefined' && module.exports) module.exports = factory();
else root.Cryosphere = factory();
})(typeof self !== 'undefined' ? self : this, function () {
'use strict';
const clamp = (v, a, b) => v < a ? a : v > b ? b : v;
const YEAR = 365 * 86400;
// make a Cryosphere system. config:
// { N, elevation(x,y), seaLevel, baseTempC, lapseC, freezeC, meltC, accumRate, meltRate,
// feedbackK (albedo-cooling strength, degC per unit ice-fraction), iceThresh, glaciateAt, deglaciateAt,
// insolSensitivity, stellar, seaPerVolume, carveRate, forcing (number | fn(clock)->degC), maxDtYears }
function makeCryosphere(config) {
config = config || {};
const N = config.N || 64;
const seaLevel = config.seaLevel != null ? config.seaLevel : 0;
const baseTempC = config.baseTempC != null ? config.baseTempC : 14; // sea-level temp at zero forcing, no ice
const lapseC = config.lapseC != null ? config.lapseC : 26; // degC colder from sea level up to elev = 1
const freezeC = config.freezeC != null ? config.freezeC : 0; // snow accumulates below this
const meltC = config.meltC != null ? config.meltC : 0; // ice melts above this
const accumRate = config.accumRate != null ? config.accumRate : 0.55; // ice units / year where cold
const meltRate = config.meltRate != null ? config.meltRate : 0.42; // ice units / (degC * year) where warm
const feedbackK = config.feedbackK != null ? config.feedbackK : 13; // the ICE-ALBEDO feedback: degC cooling at full ice
const iceThresh = config.iceThresh != null ? config.iceThresh : 0.5; // thickness that counts a cell as glaciated
const glaciateAt = config.glaciateAt != null ? config.glaciateAt : 0.30; // ice-fraction that fires 'glaciation'
const deglaciateAt = config.deglaciateAt != null ? config.deglaciateAt : 0.12;
const latGradient = config.latGradient != null ? config.latGradient : 0; // degC colder toward the poles (|y|=1); 0 = uniform (analytic test)
const insolSensitivity = config.insolSensitivity != null ? config.insolSensitivity : 40;
const stellarName = config.stellar || 'stellar';
const maxIce = config.maxIce != null ? config.maxIce : 3.5; // ice sheets saturate (a few km) -- accumulation caps here, so melt can undo it
const maxSeaDrop = config.maxSeaDrop != null ? config.maxSeaDrop : 0.08; // sea-level drop (world elev units) at full glaciation
const carveRate = config.carveRate != null ? config.carveRate : 0.004; // bedrock ground away / year under thick ice
const flowK = config.flowK != null ? config.flowK : 0; // GLACIAL FLOW: ice creeps downhill (0 = off; the analytic emergence test leaves it off)
const albedoLand = config.albedoLand != null ? config.albedoLand : 0.22;
const albedoIce = config.albedoIce != null ? config.albedoIce : 0.62;
const forcingCfg = config.forcing != null ? config.forcing : 0; // extra degC applied everywhere (a driver ramps it)
const maxDtYears = config.maxDtYears != null ? config.maxDtYears : 200; // per-step cap; huge for deep time
const elevFn = config.elevation || ((x, y) => 0.2);
const gx = (i) => (i / (N - 1)) * 2 - 1;
const forcingOf = (clock) => (typeof forcingCfg === 'function' ? forcingCfg(clock) : forcingCfg);
return {
name: config.name || 'cryosphere',
domain: 'physical',
scale: 'geological',
params: config,
init() {
const elev = new Float32Array(N * N);
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) elev[j * N + i] = elevFn(gx(i), gx(j));
return {
N, elev,
ice: new Float32Array(N * N),
carve: new Float32Array(N * N),
iceFraction: 0, volume: 0, meanAlbedo: albedoLand, seaLevelDrop: 0,
snowElev: 1, forcing: forcingOf({ t: 0, years: () => 0 }), event: null
};
},
step(st, ctx, dt) {
const dy = clamp(dt / YEAR, 0, maxDtYears); if (dy <= 0) return st;
const forcing = forcingOf(ctx.clock);
const stel = ctx.fields[stellarName];
const insolOffset = stel && stel.insolation != null ? insolSensitivity * (stel.insolation - 1) : 0;
const albedoCooling = feedbackK * st.iceFraction; // FEEDBACK read from the PREVIOUS fraction (one-step delay -> stable)
const { elev, ice, carve } = st;
// 1. MASS BALANCE: snow accumulates below freezing, ice melts above (degree-year), per cell
for (let k = 0; k < N * N; k++) {
if (elev[k] < seaLevel) { ice[k] = 0; continue; } // land ice only (sea ice is a later rung)
const yAbs = Math.abs((((k / N) | 0) / (N - 1)) * 2 - 1); // 0 at the equator .. 1 at the poles
const tCell = baseTempC + forcing + insolOffset - lapseC * Math.max(0, elev[k] - seaLevel) - latGradient * yAbs - albedoCooling;
let d = 0;
if (tCell < freezeC) d += accumRate * dy * clamp((freezeC - tCell) / 10, 0.15, 1); // colder air, heavier snow (capped)
if (tCell > meltC) d -= meltRate * (tCell - meltC) * dy; // degree-year melt
ice[k] = clamp(ice[k] + d, 0, maxIce); // saturates at maxIce (a bounded sheet, meltable)
}
// 2. GLACIAL FLOW (optional): ice creeps downhill on its OWN surface (bed + ice) -- mass-conserving,
// stability-capped. Ice pools in valleys and flows past the snowline to a melting terminus.
if (flowK > 0) {
const surf = st._surf || (st._surf = new Float32Array(N * N));
const dh = st._dh || (st._dh = new Float32Array(N * N));
for (let k = 0; k < N * N; k++) { surf[k] = (elev[k] < seaLevel ? seaLevel : elev[k]) + ice[k]; dh[k] = 0; }
const mf = Math.min(0.22, flowK * dy); // max fraction of a cell's ice that leaves per step (stable)
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) {
const k = j * N + i; if (ice[k] <= 0.01 || elev[k] < seaLevel) continue;
const nb = [i > 0 ? k - 1 : -1, i < N - 1 ? k + 1 : -1, j > 0 ? k - N : -1, j < N - 1 ? k + N : -1];
const w = [0, 0, 0, 0]; let tot = 0;
for (let n = 0; n < 4; n++) { const kk = nb[n]; if (kk < 0 || elev[kk] < seaLevel) continue; const dd = surf[k] - surf[kk]; if (dd > 0) { w[n] = dd; tot += dd; } }
if (tot <= 0) continue;
const out = mf * ice[k];
for (let n = 0; n < 4; n++) { if (w[n] <= 0) continue; const amt = out * (w[n] / tot); dh[k] -= amt; dh[nb[n]] += amt; }
}
for (let k = 0; k < N * N; k++) if (elev[k] >= seaLevel) ice[k] = clamp(ice[k] + dh[k], 0, maxIce);
}
// 3. COUNT + CARVE (after flow, so carving concentrates where the glaciers actually are -> valleys)
let iced = 0, land = 0, vol = 0, snowElev = 1;
for (let k = 0; k < N * N; k++) {
if (elev[k] < seaLevel) continue; land++;
if (ice[k] > iceThresh) {
iced++; vol += ice[k];
if (elev[k] < snowElev) snowElev = elev[k]; // lowest iced land = the snowline
carve[k] += carveRate * Math.min(ice[k], 4) * dy; // thick moving ice grinds the bed
}
}
const frac = land ? iced / land : 0;
st.event = null;
if (frac >= glaciateAt && st.iceFraction < glaciateAt) st.event = 'glaciation';
else if (frac <= deglaciateAt && st.iceFraction > deglaciateAt) st.event = 'deglaciation';
st.iceFraction = frac;
st.volume = vol;
st.seaLevelDrop = land ? (vol / (land * maxIce)) * maxSeaDrop : 0; // normalized: 0 (no ice) .. maxSeaDrop (full)
st.meanAlbedo = albedoLand + (albedoIce - albedoLand) * frac;
st.snowElev = iced ? snowElev : 1;
st.forcing = forcing;
return st;
},
// OUT: ice thickness grid + the summary scalars other systems read (sea withdrawal, albedo, snowline, carving)
fields(st) {
return {
N: st.N, ice: st.ice, carve: st.carve,
iceFraction: st.iceFraction, meanAlbedo: st.meanAlbedo,
seaLevelDrop: st.seaLevelDrop, snowElev: st.snowElev, volume: st.volume, forcing: st.forcing
};
},
events(st) { return st.event ? [{ kind: st.event, iceFraction: st.iceFraction }] : []; }
};
}
return { makeCryosphere };
});