cryosphere-test.js
133 lines · in natura-src-1.0.0.zip as natura/cryosphere-test.js · a test
// cryosphere-test.js -- assertions for the CRYOSPHERE (land ice, the ice-albedo feedback, ice-age hysteresis).
const N = require('./natura.js');
const { makeCryosphere } = require('./cryosphere.js');
const { makeStellar } = require('./stellar.js');
const { World, UNITS } = N;
const GRID = 48, YEAR = UNITS.YEAR;
let pass = 0, fail = 0;
const ok = (c, m) => { if (c) pass++; else { fail++; console.log(' FAIL: ' + m); } };
// a central mountain falling to lowland at the rim -- a broad spread of elevations, all above sea level
const EL = (x, y) => 0.55 * Math.max(0, 1 - Math.hypot(x, y) / 1.3) + 0.05;
function build(forcing, seed, extra) {
return World({ seed: seed == null ? 3 : seed, systems: [
makeCryosphere(Object.assign({ N: GRID, elevation: EL, seaLevel: 0, forcing }, extra))
] });
}
// step a world until a test on its cryosphere field returns true (or a cap), return the field snapshot then
function until(w, test, cap) {
cap = cap || 4000;
for (let i = 0; i < cap; i++) { w.step(YEAR); const f = w.fields().cryosphere; if (test(f, w)) return f; }
return w.fields().cryosphere;
}
console.log('cryosphere assertions:');
// 1. DETERMINISM -- same seed + same dt sequence -> identical ice grid
{
const a = build(-4); const b = build(-4);
for (let i = 0; i < 300; i++) { a.step(YEAR); b.step(YEAR); }
const ia = a.fields().cryosphere.ice, ib = b.fields().cryosphere.ice;
let same = true; for (let k = 0; k < GRID * GRID; k++) if (ia[k] !== ib[k]) { same = false; break; }
ok(same, 'same seed -> identical ice field (reproducible)');
}
// 2. BOUNDEDNESS / FINITENESS -- ice >= 0 finite; iceFraction in [0,1]; albedo between land and ice
{
const w = build((clk) => 4 - 0.03 * clk.years()); // sweep warm -> cold
let okAll = true, fr = [], al = [];
for (let i = 0; i < 500; i++) { w.step(YEAR); const f = w.fields().cryosphere;
for (let k = 0; k < GRID * GRID; k++) if (!(f.ice[k] >= 0 && isFinite(f.ice[k]))) okAll = false;
fr.push(f.iceFraction); al.push(f.meanAlbedo);
}
ok(okAll, 'ice thickness stays finite and non-negative everywhere');
ok(fr.every(x => x >= 0 && x <= 1), 'ice fraction stays in [0,1]');
ok(al.every(x => x >= 0.22 - 1e-6 && x <= 0.62 + 1e-6), 'mean albedo stays between land and ice albedo');
}
// 3. MOUNTAINS GLACIATE FIRST -- at a partial-ice climate the iced cells sit HIGHER than the bare ones
{
const w = build(-2);
const f = until(w, (f) => f.iceFraction > 0.10 && f.iceFraction < 0.55, 2000);
const ice = f.ice; let ihi = 0, in_ = 0, lho = 0, ln = 0;
for (let j = 0; j < GRID; j++) for (let i = 0; i < GRID; i++) {
const k = j * GRID + i, e = EL((i / (GRID - 1)) * 2 - 1, (j / (GRID - 1)) * 2 - 1);
if (ice[k] > 0.5) { ihi += e; in_++; } else { lho += e; ln++; }
}
ok(in_ > 0 && ln > 0 && ihi / in_ > lho / ln + 0.05,
'ice sits on the high ground first (iced mean elev ' + (ihi / in_).toFixed(2) + ' vs bare ' + (lho / ln).toFixed(2) + ')');
}
// 4. THE EMERGENCE -- ICE-ALBEDO HYSTERESIS (an ice age is a TIPPING, not a smooth response)
// Cool the world and find the forcing where ice takes over; then, from a glaciated world, warm it and find
// the forcing where the ice lets go. A bistable system must be warmed PAST where it froze: F_deglac > F_glac.
{
// cooling sweep: forcing falls slowly from +6; record the forcing at 10% / 50% / 90% ice cover
const cool = build((clk) => 6 - 0.02 * clk.years());
let F10 = null, F50 = null, F90 = null;
for (let i = 0; i < 2000; i++) { cool.step(YEAR); const f = cool.fields().cryosphere;
if (F10 == null && f.iceFraction >= 0.10) F10 = f.forcing;
if (F50 == null && f.iceFraction >= 0.50) F50 = f.forcing;
if (F90 == null && f.iceFraction >= 0.90) { F90 = f.forcing; break; }
}
// warming sweep: hold very cold to a glaciated world for 500yr, then warm; record the forcing at deglaciation
const warm = build((clk) => { const y = clk.years(); return y < 500 ? -16 : -16 + 0.02 * (y - 500); });
let Fdeglac = null;
for (let i = 0; i < 3000; i++) { warm.step(YEAR); const f = warm.fields().cryosphere;
if (warm.clock.years() > 520 && f.iceFraction <= 0.5) { Fdeglac = f.forcing; break; }
}
ok(F50 != null && Fdeglac != null, 'the world both glaciates on cooling and deglaciates on warming');
ok(F50 != null && Fdeglac != null && Fdeglac > F50 + 1.5,
'HYSTERESIS: must warm past where it froze (glaciate at ' + (F50 || 0).toFixed(1) + 'C forcing, deglaciate at ' + (Fdeglac || 0).toFixed(1) + 'C)');
ok(F10 != null && F90 != null && Math.abs(F10 - F90) < 3.5,
'the flip is a SHARP tip: 10%->90% ice over ' + Math.abs((F10 || 0) - (F90 || 0)).toFixed(1) + 'C (a smooth response would need the full ~14C lapse range)');
}
// 5. STELLAR COUPLING -- a fainter (young) sun gives a colder world -> more ice at the same forcing
{
const buildS = (age0Gyr) => World({ seed: 3, systems: [
makeStellar({ age0Gyr }),
makeCryosphere({ N: GRID, elevation: EL, seaLevel: 0, forcing: 0, stellar: 'stellar' })
] });
const faint = buildS(0.3); for (let i = 0; i < 400; i++) faint.step(YEAR);
const bright = buildS(4.6); for (let i = 0; i < 400; i++) bright.step(YEAR);
const ff = faint.fields().cryosphere.iceFraction, bf = bright.fields().cryosphere.iceFraction;
ok(ff > bf + 0.05, 'a fainter sun ices more of the world (' + (ff * 100).toFixed(0) + '% faint vs ' + (bf * 100).toFixed(0) + '% present)');
}
// 6. OUT-COUPLINGS -- ice locks up water (the sea withdraws) and grinds the bed (carving accumulates)
{
const w = build(-10); for (let i = 0; i < 300; i++) w.step(YEAR);
const f = w.fields().cryosphere;
ok(f.seaLevelDrop > 0 && f.iceFraction > 0.3, 'a glaciated world drops the sea level (drop ' + f.seaLevelDrop.toFixed(2) + ' at ice ' + (f.iceFraction * 100).toFixed(0) + '%)');
let carved = 0; for (let k = 0; k < GRID * GRID; k++) if (f.carve[k] > 0) carved++;
ok(carved > 0, 'thick ice grinds the bedrock (' + carved + ' cells carved)');
// and a warm world holds no ice
const warm = build(20); for (let i = 0; i < 200; i++) warm.step(YEAR);
ok(warm.fields().cryosphere.iceFraction === 0, 'a hot world holds no ice (fraction 0)');
}
// 7. GLACIAL FLOW: ice creeps downhill, conserving mass (opt-in flowK; default 0 leaves the above untouched)
{
const ramp = (x, y) => 0.5 * (1 - (x + 1) / 2) + 0.05; // high at x=-1 (col 0), low at x=+1 (last col)
const w = World({ seed: 3, systems: [ makeCryosphere({ N: GRID, elevation: ramp, seaLevel: 0, accumRate: 0, meltRate: 0, flowK: 0.05, maxIce: 20, forcing: 0 }) ] });
const ice = w.get('cryosphere').ice;
for (let j = 0; j < GRID; j++) for (let i = 0; i < GRID; i++) if (i < GRID * 0.3) ice[j * GRID + i] = 3; // a lump of ice on the high side
const sum = (a) => { let s = 0; for (let k = 0; k < a.length; k++) s += a[k]; return s; };
const comX = (a) => { let m = 0, mx = 0; for (let j = 0; j < GRID; j++) for (let i = 0; i < GRID; i++) { const v = a[j * GRID + i]; m += v; mx += v * i; } return m ? mx / m : 0; };
const before = sum(ice), cx0 = comX(ice);
for (let s = 0; s < 200; s++) w.step(YEAR);
const after = sum(ice), cx1 = comX(ice);
ok(Math.abs(after - before) < before * 0.002, 'glacial flow conserves ice mass (' + before.toFixed(1) + ' -> ' + after.toFixed(1) + ')');
ok(cx1 > cx0 + 1.5, 'ice flows DOWNHILL (centre of mass moves down-slope, col ' + cx0.toFixed(1) + ' -> ' + cx1.toFixed(1) + ')');
// and with NO flow the lump stays put
const w0 = World({ seed: 3, systems: [ makeCryosphere({ N: GRID, elevation: ramp, seaLevel: 0, accumRate: 0, meltRate: 0, flowK: 0, maxIce: 20, forcing: 0 }) ] });
const ice0 = w0.get('cryosphere').ice;
for (let j = 0; j < GRID; j++) for (let i = 0; i < GRID; i++) if (i < GRID * 0.3) ice0[j * GRID + i] = 3;
const cxA = comX(ice0); for (let s = 0; s < 200; s++) w0.step(YEAR); const cxB = comX(ice0);
ok(Math.abs(cxB - cxA) < 0.5, 'with flow OFF the ice stays where it formed (no drift: ' + cxA.toFixed(1) + ' -> ' + cxB.toFixed(1) + ')');
}
console.log(`\n${pass} passed, ${fail} failed`);
process.exit(fail ? 1 : 0);