soil.js
758 lines · in natura-src-1.0.0.zip as natura/soil.js · a module · its entry in the module reference · its test: soil-test.js
// soil.js -- Natura SOIL: pedogenesis and the matter cycle, the living skin between rock and life. Most
// worlds treat "fertility" as one number painted on the ground; this treats the SOIL as a real body that
// GROWS from bedrock, RECORDS the life that lived on it, MOVES downhill as mass, and feeds the plants that
// in turn build it. It is the slow verb of the land: stone becomes grit becomes soil becomes food becomes
// leaf becomes soil, and a hillside sorts itself into bare ridges and fat valley floors because the matter
// literally runs downhill. Nothing here is a scalar knob -- every cell carries a column of real stuff whose
// mass is conserved, so the emergences (a soil CATENA, fertile bottoms, a land that remembers its biology,
// a mass-balanced loop, and the wound a fire leaves on bare ground) FALL OUT of the mechanism rather than
// being written in.
//
// The mechanisms, all coupled on an N x N terrain grid:
// 1. WEATHERING -- bedrock -> regolith -> mineral soil. PHYSICAL (freeze-thaw, driven by temperature
// swings across 0C) makes coarse grit (sand/silt); CHEMICAL (dissolution/oxidation, driven by
// moisture x warmth) makes clay and frees a trickle of mineral nutrient. Rate gated by climate + slope.
// 2. HORIZONS / DEPTH -- a real soil column: a LITTER layer (fresh dead matter), a HUMUS store (the dark
// decomposed carbon), a MINERAL soil depth, over WEATHERED regolith, over BEDROCK.
// 3. ORGANIC CYCLE -- litter in (from the vegetation/carcass driver) -> DECOMPOSITION (moisture+temp) ->
// humus -> MINERALIZATION -> plant-available NITROGEN and PHOSPHORUS, each with uptake + leaching.
// 4. MOISTURE + TEXTURE -- water-holding capacity rises with CLAY and humus; sandy/steep ground DRAINS and
// leaches. Texture itself is an emergent product of which weathering won.
// 5. EROSION + DEPOSITION -- slope + surface water strip the topsoil (humus, fine mineral, and the
// nutrients riding in them) off ridges and lay them down in the valleys. Mass MOVES, never vanishes.
// 6. VEGETATION FEEDBACK (two-way) -- soil depth + fertility + moisture GATE plant growth; plant growth and
// death BUILD soil (litter above, roots mixed below). The land and its life co-produce each other.
//
// Ships with its OWN self-contained drivers (a slope heightfield, a moisture+temperature climate, and a
// soil-responsive vegetation/litter source) so it RUNS and TESTS standalone -- but it reads ctx.fields for
// 'terrain', 'weather', and 'vegetation' when present, so it COUPLES to the real systems unchanged.
//
// A Natura system spec (init/step/fields/events), deterministic by construction: randomness comes ONLY from
// init(rng) and ctx.rng, time is dt (seconds) and ctx.clock. UMD (node test / window.Soil).
(function (root, factory) {
if (typeof module !== 'undefined' && module.exports) module.exports = factory();
else root.Soil = factory();
})(typeof self !== 'undefined' ? self : this, function () {
'use strict';
const clamp = (v, a, b) => v < a ? a : v > b ? b : v;
const smooth = (v) => v <= 0 ? 0 : v >= 1 ? 1 : v * v * (3 - 2 * v); // smoothstep 0..1
// make a Soil system. config below; every driver is overridable so the soil can be fed by the real
// terrain/weather/vegetation systems instead of its built-in stand-ins.
function makeSoil(config) {
config = config || {};
const N = config.N || 32; // grid resolution (N x N cells)
// --- which upstream systems to read when they are present (else the built-in drivers run) ------------
const terrainName = config.terrain || 'terrain';
const weatherName = config.weather || 'weather';
const vegName = config.vegetation || 'vegetation';
// --- WEATHERING rates (per sim-year; pedogenesis is slow, so soil depth reaches O(1) over centuries) --
const physRate = config.physRate != null ? config.physRate : 0.012; // freeze-thaw bedrock->regolith (depth-equiv/yr at full drive)
const chemRate = config.chemRate != null ? config.chemRate : 0.02; // chemical regolith ripening (per yr at full drive)
const regoToSoil = config.regoToSoil != null ? config.regoToSoil : 0.03; // slow mechanical regolith -> mineral soil (per yr)
const nutFromRock = config.nutFromRock != null ? config.nutFromRock : 0.08; // mineral nutrient freed per unit chem weathering
const pFromRock = config.pFromRock != null ? config.pFromRock : 0.5; // phosphorus is rock-derived, so a substantial share comes from weathering
// --- ORGANIC cycle rates (per sim-year) -------------------------------------------------------------
const decompRate = config.decompRate != null ? config.decompRate : 8.0; // litter -> humus at full moisture+warmth
const humifyFrac = config.humifyFrac != null ? config.humifyFrac : 0.35; // fraction of decomposed matter that becomes stable humus (rest respired)
const mineralizeRate = config.mineralizeRate != null ? config.mineralizeRate : 0.12; // organic N/P -> plant-available (per yr at full drive)
// --- MOISTURE + TEXTURE -----------------------------------------------------------------------------
const rainK = config.rainK != null ? config.rainK : 4.0; // rain -> soil water recharge (a moist climate refills the profile each year)
const baseWHC = config.baseWHC != null ? config.baseWHC : 0.15; // water-holding of pure sand (fraction)
const clayWHC = config.clayWHC != null ? config.clayWHC : 0.55; // extra WHC at pure clay
const humusWHC = config.humusWHC != null ? config.humusWHC : 0.6;// humus greatly boosts water holding
const leachRate = config.leachRate != null ? config.leachRate : 0.9; // nutrient loss with drainage (per yr at full drainage)
// --- EROSION + DEPOSITION ---------------------------------------------------------------------------
const erodeK = config.erodeK != null ? config.erodeK : 0.02; // stream-power topsoil detachment coeff
const flowM = config.flowM != null ? config.flowM : 0.5; // discharge exponent
const slopeN = config.slopeN != null ? config.slopeN : 1.1; // slope exponent
const coverProtect = config.coverProtect != null ? config.coverProtect : 0.9; // vegetation cover fraction that resists erosion
const settleFrac = config.settleFrac != null ? config.settleFrac : 0.7; // fraction of through-flowing load that settles per cell of travel
// --- VEGETATION FEEDBACK (the built-in driver; skipped when a real vegetation field is present) -------
const sterile = !!config.sterile; // a lifeless world: no plants at all (for isolating weathering vs biology)
const bioRate = config.bioRate != null ? config.bioRate : 3.0; // logistic plant growth rate (per yr)
const bioRespire = config.bioRespire != null ? config.bioRespire : 0.5; // baseline biomass turnover (per yr) -> litter+roots
const rootFrac = config.rootFrac != null ? config.rootFrac : 0.35; // share of dead biomass entering soil as roots (rest is surface litter)
const uptakeN = config.uptakeN != null ? config.uptakeN : 0.4; // N drawn per unit biomass grown
const uptakeP = config.uptakeP != null ? config.uptakeP : 0.04; // P drawn per unit biomass grown (plants need far less P than N)
const nHalf = config.nHalf != null ? config.nHalf : 0.15; // N half-saturation for growth
const pHalf = config.pHalf != null ? config.pHalf : 0.015; // P half-saturation for growth
const depthHalf = config.depthHalf != null ? config.depthHalf : 0.25; // soil-depth half-saturation for rooting
const moistHalf = config.moistHalf != null ? config.moistHalf : 0.25; // moisture half-saturation for growth
// --- built-in DRIVERS (a terrain heightfield, a climate, a litter source) ---------------------------
// A HILLSLOPE draining into an INTERIOR VALLEY FLOOR: the top rows are a high ridge, the ground falls
// away to a low flat valley in the lower-middle, and the very bottom rows rise slightly into a lip so
// the valley is a genuine interior sink (eroded soil that reaches it stays, rather than exporting off
// the grid edge). A single low outlet notch in the bottom edge lets the trickle that overtops the
// valley leave the world. Used only when no external terrain field is present. Height in [0,1].
const heightFn = config.height || function (i, j) {
const f = j / (N - 1); // 0 at top ridge .. 1 at bottom edge
const valleyEnd = 0.86; // the valley floor bottoms out at ~86% down
let h;
if (f <= valleyEnd) {
// a smooth CONCAVE hillslope: steep near the ridge, flattening as it approaches the valley floor
// (the classic soil-catena profile). u in [0,1] along the hillslope; h uses u^1.7 so the gradient
// eases toward the bottom, giving a genuinely low-slope valley floor rather than an abrupt step.
const u = f / valleyEnd;
h = 0.32 + 0.68 * Math.pow(1 - u, 1.7);
} else {
// a gentle lip rising below the valley, damming the basin so deposited soil is trapped
const v = (f - valleyEnd) / (1 - valleyEnd);
h = 0.32 + 0.5 * (v * v);
}
// one outlet notch: a single cell of the bottom edge is cut low so the basin drains through ONE point
if (j === N - 1 && i === (N >> 1)) h = 0.28;
return h;
};
// Climate: a mean temperature and a diurnal/seasonal SWING (freeze-thaw needs the swing to cross 0C),
// and a moisture supply. Overridable; defaults are a cool temperate montane climate.
const climate = config.climate || {};
const meanTemp = climate.meanTemp != null ? climate.meanTemp : 6; // deg C annual mean (cool -> freeze-thaw active)
const tempSwing = climate.tempSwing != null ? climate.tempSwing : 14; // deg C amplitude of the seasonal swing
const diurnalSwing = climate.diurnalSwing != null ? climate.diurnalSwing : 8; // deg C day-night amplitude
const moistSupply = climate.moisture != null ? climate.moisture : 0.5; // baseline moisture supply 0..1
return {
name: config.name || 'soil',
domain: 'geology',
scale: 'pedological',
params: config,
// --- STATE: every cell carries a real column of matter ------------------------------------------
init(rng, params, clock) {
const n2 = N * N;
const st = {
N, t: 0,
// terrain / geometry
bedrock: new Float32Array(n2), // bedrock surface elevation (the unweathered base)
slope: new Float32Array(n2),
flow: new Float32Array(n2), // drainage-area proxy (D8 accumulation)
land: new Uint8Array(n2),
// the soil COLUMN (all in mass-per-column units)
regolith: new Float32Array(n2), // weathered rock not yet true soil
mineral: new Float32Array(n2), // mineral soil depth (the B/C horizon body)
humus: new Float32Array(n2), // stable organic carbon (the dark A horizon)
litter: new Float32Array(n2), // fresh undecomposed dead matter (the O horizon)
// TEXTURE -- fractions of the mineral fraction (sand+silt+clay track separately, sum ~ mineral)
sand: new Float32Array(n2), silt: new Float32Array(n2), clay: new Float32Array(n2),
// NUTRIENTS -- plant-available pools (organic N/P live inside litter+humus via their content)
availN: new Float32Array(n2), availP: new Float32Array(n2),
organicN: new Float32Array(n2), organicP: new Float32Array(n2), // N/P locked in litter+humus
plantN: new Float32Array(n2), plantP: new Float32Array(n2), // N/P taken up and held in living tissue
// WATER
water: new Float32Array(n2), // current soil water store
whc: new Float32Array(n2), // water-holding capacity (derived from texture+humus)
moistIdx: new Float32Array(n2), // saturation fraction 0..1 (the biologically relevant wetness)
// VEGETATION (built-in driver state; overwritten by a real veg field when present)
bio: new Float32Array(n2), // standing plant biomass
cleared: null, // optional Uint8 mask: cells held bare (fire/clearing/grazing)
// per-step climate scratch
temp: new Float32Array(n2),
// per-cell CUMULATIVE erosion/deposition ledgers -- the record of where the land lost soil and
// where it gained it, so the catena can be read straight off the geomorphic budget.
cumEroded: new Float32Array(n2), // total topsoil mass detached and carried OUT of this cell
cumDeposited: new Float32Array(n2), // total mass settled INTO this cell from upslope
// book-keeping / conservation ledgers (cumulative). Per-element and split by pathway so the
// closed-loop mass balance is EXACTLY checkable: for N and for P,
// final stock == (rock source + atmospheric source) - (leached + exported)
// and nothing is created or destroyed anywhere else.
ledger: {
weatheredRock: 0, // bedrock mass converted to regolith (rock -> soil column, a mass source)
rockN: 0, rockP: 0, // N, P freed from rock by chemical weathering + continued dissolution (SOURCE)
atmN: 0, // N from the atmosphere: baseline deposition + biological fixation (SOURCE)
leachedN: 0, leachedP: 0, // N, P flushed out of the column by drainage (SINK)
erodedN: 0, erodedP: 0, // N, P carried off the grid's outlet riding on eroded soil (SINK)
erodedOut: 0, // soil MASS carried off the grid's outlet (mass SINK)
respiredC: 0 // organic carbon respired to the air (CO2; a carbon sink, tracked)
},
_scratch: {
mob: new Float32Array(n2), mobH: new Float32Array(n2), mobN: new Float32Array(n2),
mobP: new Float32Array(n2), mobSand: new Float32Array(n2), mobSilt: new Float32Array(n2),
mobClay: new Float32Array(n2), order: new Int32Array(n2)
},
evs: [],
_haveTerrain: false, _haveVeg: false
};
// lay the bedrock from the driver (or the coupled terrain when it arrives in step 0)
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) {
const k = j * N + i;
st.bedrock[k] = heightFn(i, j);
st.land[k] = 1;
// a whisper of noise so the field is not perfectly degenerate (deterministic, from init rng)
st.bedrock[k] += (rng() - 0.5) * 0.01;
}
// remember the ORIGINAL terrain relief -- the landscape's inherent steepness, the driver of where
// erosion bites. (The live bedrock lowers as it weathers; this preserved copy is the clean
// reference for "how steep is this ground," undistorted by deposition roughening the live surface.)
st.bedrock0 = st.bedrock.slice();
st.initSlope = new Float32Array(n2);
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) {
const il = Math.max(0, i - 1), ir = Math.min(N - 1, i + 1), jd = Math.max(0, j - 1), ju = Math.min(N - 1, j + 1);
const dx = (st.bedrock0[j * N + ir] - st.bedrock0[j * N + il]) * 0.5 * N;
const dz = (st.bedrock0[ju * N + i] - st.bedrock0[jd * N + i]) * 0.5 * N;
st.initSlope[j * N + i] = Math.hypot(dx, dz);
}
computeSlopeFlow(st);
return st;
},
step(st, ctx, dt) {
st.evs.length = 0;
const dy = clamp((dt || 0) / 31536000, 0, 5); // sim-YEARS this step (pedogenesis is slow)
if (dy <= 0) return st;
st.t += dt;
// ---- 0. INGEST external drivers when present (else the built-ins supply the same fields) -------
ingestTerrain(st, ctx);
const clim = sampleClimate(st, ctx, dt); // fills st.temp, returns {moist[], phys[], warm[]}
ingestVegetation(st, ctx, dy, clim); // sets st.bio + returns litter/root inputs via st
const n2 = N * N;
const { bedrock, regolith, mineral, humus, litter, sand, silt, clay,
availN, availP, organicN, organicP, water, whc, moistIdx, bio, land } = st;
const L = st.ledger;
// ================================================================================================
// 1. WEATHERING -- bedrock into regolith into mineral soil.
// PHYSICAL: freeze-thaw. Its drive is how hard the temperature swing crosses 0C (a cell that
// never freezes, or never thaws, does not shatter). Coarse product (sand/silt).
// CHEMICAL: dissolution/oxidation. Its drive is moisture x warmth (a warm wet soil ripens; a
// frozen or bone-dry one does not). Fine product (clay) and a trickle of freed nutrient.
// Both slow as the soil mantle thickens (a deep soil shields its own bedrock -- the classic
// humped/â– declining weathering-vs-depth curve), and physical gets a slope kick (steep faces
// shed their debris and keep exposing fresh rock).
// ================================================================================================
for (let k = 0; k < n2; k++) {
if (!land[k]) continue;
const phys = clim.phys[k]; // 0..1 freeze-thaw intensity this step
const warm = clim.warm[k]; // 0..1 warmth for chemistry
const moist = moistIdx[k]; // previous-step wetness drives chemistry
const mantle = regolith[k] + mineral[k] + humus[k];
const shield = Math.exp(-mantle / 0.6); // thick soil shields the rock beneath
const slopeKick = 1 + 1.5 * clamp(st.slope[k], 0, 1); // steep ground keeps exposing fresh rock
// PHYSICAL: bedrock -> regolith. Freeze-thaw shatters rock into coarse fragments; regolith is
// untextured "proto-mineral" (texture is only assigned when it becomes true soil, below).
const dPhys = physRate * phys * shield * slopeKick * dy;
const takenP = Math.min(dPhys, Math.max(0, bedrock[k] - (config.floor != null ? config.floor : -1)));
if (takenP > 0) { bedrock[k] -= takenP; regolith[k] += takenP; L.weatheredRock += takenP; }
// CHEMICAL: regolith ripens into mineral soil, and dissolution mints CLAY (fine, clay-rich
// product) while freeing a trickle of rock nutrient. This is the fine-textured pathway.
const dChem = chemRate * warm * moist * shield * dy;
const takenC = Math.min(dChem, regolith[k]);
if (takenC > 0) {
regolith[k] -= takenC; mineral[k] += takenC;
addTexture(st, k, takenC, 0.15, 0.30, 0.55); // chemistry -> clay-biased soil
const freedN = nutFromRock * takenC; // freed mineral N
const freedP = nutFromRock * takenC * pFromRock; // P is disproportionately rock-sourced
availN[k] += freedN; availP[k] += freedP;
L.rockN += freedN; L.rockP += freedP;
}
// MECHANICAL comminution: regolith -> mineral even without full chemistry, staying COARSE
// (sand/silt-rich). This is the pathway that dominates in cold/dry places where chemistry stalls,
// so those soils come out sandy while warm wet soils come out clayey -- texture records climate.
const takenM = Math.min(regoToSoil * dy * regolith[k], regolith[k]);
if (takenM > 0) {
regolith[k] -= takenM; mineral[k] += takenM;
addTexture(st, k, takenM, 0.55, 0.35, 0.10); // comminution -> coarse soil
}
// CONTINUED mineral dissolution: primary grains in the mineral body keep slowly releasing nutrient
// (especially phosphorus) long after the regolith is gone -- the base-level supply that keeps a
// mature soil from fully starving. Warm-wet soils release faster; a cold or dry soil hoards its
// rock nutrient. This is why deeply weathered ancient soils drift toward P limitation but do not
// hit zero. Ledgered as a rock source.
const dSolute = 0.006 * warm * moist * mineral[k] * dy;
if (dSolute > 0) {
const fP = dSolute * pFromRock, fN = dSolute * (1 - pFromRock) * 0.3;
availP[k] += fP; availN[k] += fN;
L.rockP += fP; L.rockN += fN;
}
}
// ================================================================================================
// 4a. TEXTURE-DERIVED WATER-HOLDING + this step's soil MOISTURE (needed before the organic cycle,
// which is moisture-gated). The water-holding CAPACITY per unit of the biologically-active
// rooting zone rises with clay fraction and humus (a clayey, humus-rich soil holds far more
// water than bare sand); the store recharges from climate moisture and DRAINS faster on sandy,
// steep ground. moistIdx is the saturation of that rooting-zone capacity -- the wetness biology
// and chemistry actually feel -- so a wet climate stays wet regardless of how deep the soil is.
// ================================================================================================
st._drain = st._drain || new Float32Array(n2);
for (let k = 0; k < n2; k++) {
if (!land[k]) continue;
const minMass = Math.max(mineral[k], 1e-6);
const clayF = clamp(clay[k] / minMass, 0, 1);
const sandF = clamp(sand[k] / minMass, 0, 1);
// capacity of the rooting zone: quality (per-unit water-holding, texture+humus driven) times the
// rooting depth available, but the depth that matters is capped (plants root the topsoil, not
// bedrock). A deep soil holds MORE total water but the saturation the roots feel saturates.
const quality = baseWHC + (clayWHC - baseWHC) * clayF + humusWHC * clamp(humus[k] / (humus[k] + 0.3), 0, 1);
const rootDepth = clamp((mineral[k] + humus[k]) / 0.5, 0, 1); // fraction of a full rooting zone present
const cap = Math.max(quality * (0.2 + 0.8 * rootDepth), 1e-3);
const supply = clim.moist[k] * rainK * cap * dy; // recharge scales to the capacity it fills
const drainFrac = clamp((0.15 + 0.4 * sandF + 0.4 * clamp(st.slope[k] / 2, 0, 1)) * dy, 0, 0.9);
let w = water[k] + supply;
let over = Math.max(0, w - cap); w -= over; // runoff of anything above field capacity
const drained = w * drainFrac; w -= drained; // slow gravitational drainage
water[k] = w;
whc[k] = cap;
moistIdx[k] = clamp(w / cap, 0, 1);
st._drain[k] = clamp((over + drained) / cap, 0, 1); // fractional through-flow -> nutrient leaching
}
// ================================================================================================
// 3. ORGANIC CYCLE -- litter decomposes to humus (moisture+temp gated), humus mineralizes to
// available N & P, plants take up nutrients (in the veg feedback below), and drainage leaches.
// Organic N & P ride inside litter and humus; when litter decomposes its nutrient content splits
// between the new humus and the freed available pool; when humus mineralizes it frees more.
// ================================================================================================
// Nutrients ride in the organicN/organicP pools (placed there at plant death); litter and humus are
// the CARBON bodies. So decomposition moves carbon litter->humus (respiring the rest), and it is
// MINERALIZATION of the organic nutrient pools that frees plant-available N & P. This keeps the
// matter conserved: nothing is minted from a mass-times-content constant.
for (let k = 0; k < n2; k++) {
if (!land[k]) continue;
const warm = clim.warm[k], moist = moistIdx[k];
const decompDrive = clamp(warm * (0.2 + 0.8 * moist), 0, 1); // decomposition needs warmth AND water
// DECOMPOSE litter carbon -> humus carbon (rest respired to the air as CO2, a tracked sink)
const dLit = clamp(decompRate * decompDrive * dy, 0, 0.98) * litter[k];
if (dLit > 0) {
const toHumus = dLit * humifyFrac;
litter[k] -= dLit; humus[k] += toHumus;
L.respiredC += dLit - toHumus;
}
// MINERALIZE: the organic nutrient pools break down into plant-available N & P at a rate set by the
// same warmth+moisture drive. A fraction of the locked organic N/P is freed each step; the carbon
// of the humus that hosted it is respired away in proportion.
const minFrac = clamp(mineralizeRate * decompDrive * dy, 0, 0.5);
if (minFrac > 0) {
const nRel = organicN[k] * minFrac, pRel = organicP[k] * minFrac;
organicN[k] -= nRel; organicP[k] -= pRel;
availN[k] += nRel; availP[k] += pRel;
const humLost = Math.min(humus[k], humus[k] * minFrac); // humus carbon leaves as the organics mineralize
humus[k] -= humLost; L.respiredC += humLost;
}
// LEACHING -- drainage flushes a fraction of the MOBILE available pool downward/out of the column.
// Nitrate is highly mobile; phosphate adsorbs to clay and leaches far less.
const flush = clamp((st._drain ? st._drain[k] : 0) * leachRate * dy, 0, 0.9);
if (flush > 0) {
const lN = availN[k] * flush, lP = availP[k] * flush * 0.08; // phosphate adsorbs to clay: nearly immobile
availN[k] -= lN; availP[k] -= lP;
L.leachedN += lN; L.leachedP += lP;
}
}
// ================================================================================================
// 6. VEGETATION FEEDBACK (built-in driver) -- soil GATES growth (depth, moisture, N, P all
// co-limit, Liebig-style by the minimum-ish product), and plant turnover BUILDS soil (surface
// litter + belowground roots, each carrying nutrients back). When a real vegetation field is
// coupled, we instead read its biomass and take litter from its death, but still run uptake.
// ================================================================================================
if (!st._haveVeg) {
const plantN = st.plantN, plantP = st.plantP, cleared = st.cleared;
for (let k = 0; k < n2; k++) {
if (!land[k]) continue;
const depth = mineral[k] + humus[k];
const gDepth = depth / (depth + depthHalf); // SOIL DEPTH gates rooting
const gMoist = moistIdx[k] / (moistIdx[k] + moistHalf);// MOISTURE gates growth
const gN = availN[k] / (availN[k] + nHalf); // NITROGEN availability
const gP = availP[k] / (availP[k] + pHalf); // PHOSPHORUS availability
const warm = clim.warm[k]; // warmth gates growth
const nutLim = Math.min(gN, gP); // Liebig: the scarcer nutrient limits
let habitat = clamp(gDepth * gMoist * warm, 0, 1); // the PHYSICAL habitat (depth+water+warmth)
// CLEARING / DISTURBANCE: a cell held bare (fire, grazing, tillage) grows nothing -- no
// recruitment, no fixation, no cover. This is the fire/clearing lever the (e) emergence pulls:
// strip the plants and the naked topsoil is left to erode and leach with nothing to hold it.
if (cleared && cleared[k]) habitat = 0;
if (sterile) habitat = 0; // a lifeless world grows nothing anywhere
// carrying capacity is CO-LIMITED by the habitat and the scarcer nutrient. On deep, moist,
// fertile ground K is high; on thin bare rock it is near zero.
const K = clamp(habitat * (0.15 + 0.85 * nutLim), 0, 1);
const b = bio[k];
// baseline atmospheric N DEPOSITION -- rain and dust bring a trickle of nitrogen everywhere,
// independent of life. This is what lets the very first PIONEERS establish on raw mineral soil
// (they then fix their own N and the succession bootstraps). A real input, so it is ledgered.
const nDep = 0.004 * habitat * dy;
availN[k] += nDep; L.atmN += nDep;
// RECRUITMENT: pioneers colonize bare-but-habitable ground -- gated by the PHYSICAL habitat, NOT
// by fertility, because early colonizers fix their own nitrogen. This breaks the no-plant ->
// no-litter -> no-nutrient -> no-plant deadlock the way primary succession really does.
const recruit = 0.03 * habitat * dy;
// demand-limited logistic growth, then throttled by whether the nutrients are actually there
const potGrowth = Math.max(0, bioRate * b * (1 - b / Math.max(K, 1e-3)) * dy) + recruit;
const wantN = potGrowth * uptakeN, wantP = potGrowth * uptakeP;
const gotN = Math.min(wantN, availN[k]), gotP = Math.min(wantP, availP[k]);
const nutFrac = wantN > 1e-12 || wantP > 1e-12
? Math.min(wantN > 1e-12 ? gotN / wantN : 1, wantP > 1e-12 ? gotP / wantP : 1) : 1;
const growth = potGrowth * nutFrac; // realized growth: limited by the scarcer supply
// UPTAKE: available -> plant tissue pools (mass moves, is not created)
const takeN = growth * uptakeN, takeP = growth * uptakeP;
availN[k] -= takeN; availP[k] -= takeP; plantN[k] += takeN; plantP[k] += takeP;
const nb = clamp(b + growth, 0, 1.5);
// TURNOVER / DEATH -> litter (surface) + roots (belowground). Dead tissue carries its held N & P
// OUT of the plant pool (conserved): the humified share locks into organicN/P, the rest frees to
// available -- so the nutrients a plant took up cycle back through the soil, never vanish.
const death = bioRespire * nb * dy;
const deadFrac = nb > 1e-9 ? death / nb : 0;
const rootPart = death * rootFrac, surfPart = death - rootPart;
litter[k] += surfPart; humus[k] += rootPart * humifyFrac;
L.respiredC += rootPart * (1 - humifyFrac); // root carbon partly respired in humifying
const relN = plantN[k] * deadFrac, relP = plantP[k] * deadFrac;
plantN[k] -= relN; plantP[k] -= relP;
organicN[k] += relN * humifyFrac; organicP[k] += relP * humifyFrac;
availN[k] += relN * (1 - humifyFrac); availP[k] += relP * (1 - humifyFrac);
// biological N FIXATION: the one honest external N source (legumes/microbes pull N from the air).
// Scaled to living cover and ledgered so mass balance can subtract it as an input.
const nFix = 0.03 * nb * dy;
availN[k] += nFix; L.atmN += nFix;
bio[k] = nb - death;
}
} else {
// COUPLED path: a real vegetation field drives biomass; we run the same conserving uptake/return,
// taking gross growth and death from the coupled field's change (computed in ingestVegetation).
const plantN = st.plantN, plantP = st.plantP;
for (let k = 0; k < n2; k++) {
if (!land[k]) continue;
const grown = Math.max(0, st._vegGrowth ? st._vegGrowth[k] : 0);
const died = Math.max(0, st._vegDeath ? st._vegDeath[k] : 0);
// UPTAKE -> plant pools (mass moves, bounded by availability)
const takeN = Math.min(grown * uptakeN, availN[k]), takeP = Math.min(grown * uptakeP, availP[k]);
availN[k] -= takeN; availP[k] -= takeP; plantN[k] += takeN; plantP[k] += takeP;
// DEATH -> litter/humus carbon + nutrient return from the plant pool (conserved)
const bnow = Math.max(bio[k], 1e-9), deadFrac = clamp(died / bnow, 0, 1);
const rootPart = died * rootFrac, surfPart = died - rootPart;
litter[k] += surfPart; humus[k] += rootPart * humifyFrac; L.respiredC += rootPart * (1 - humifyFrac);
const relN = plantN[k] * deadFrac, relP = plantP[k] * deadFrac;
plantN[k] -= relN; plantP[k] -= relP;
organicN[k] += relN * humifyFrac; organicP[k] += relP * humifyFrac;
availN[k] += relN * (1 - humifyFrac); availP[k] += relP * (1 - humifyFrac);
const nFix = 0.03 * bio[k] * dy; availN[k] += nFix; L.atmN += nFix;
}
}
// ================================================================================================
// 5. EROSION + DEPOSITION -- the topsoil (litter, humus, fine mineral = silt+clay, and the nutrients
// riding in them) is detached by stream power (discharge^m * slope^n), throttled by vegetation
// cover, then ROUTED down the D8 flow network. A fraction settles in each downhill cell it
// crosses (deposition), the rest continues, and whatever leaves the lowest outlets is ledgered as
// export. Mass is conserved: what leaves a ridge lands in a valley; the only true loss is the
// grid-edge outlet, and that is counted.
// ================================================================================================
erodeAndDeposit(st, dy);
// ================================================================================================
// Recompute geometry (soil built up / stripped changes the SURFACE, hence slope + flow) + events.
// ================================================================================================
computeSlopeFlow(st);
// keep organicN/P from drifting negative through rounding
for (let k = 0; k < n2; k++) { if (organicN[k] < 0) organicN[k] = 0; if (organicP[k] < 0) organicP[k] = 0; if (availN[k] < 0) availN[k] = 0; if (availP[k] < 0) availP[k] = 0; }
// EVENTS -- named soil-scale happenings a coupled world (or a renderer) can react to
let maxDepth = 0, sumHumus = 0;
for (let k = 0; k < n2; k++) { if (!land[k]) continue; const d = mineral[k] + humus[k]; if (d > maxDepth) maxDepth = d; sumHumus += humus[k]; }
if (st._erodedThisStep > 0.02) st.evs.push({ kind: 'erosion-pulse', mass: st._erodedThisStep });
if (st._exportedThisStep > 0.01) st.evs.push({ kind: 'sediment-export', mass: st._exportedThisStep });
const meanHumus = sumHumus / n2;
if (!st._maturedFired && maxDepth > 0.5) { st.evs.push({ kind: 'soil-matured', meanHumus, maxDepth }); st._maturedFired = true; }
return st;
},
// --- fields(): everything a coupled system (or a renderer) might read ------------------------------
fields(st) {
const n2 = st.N * st.N;
// derived, allocation-light summaries recomputed lazily
return {
N: st.N, land: st.land,
bedrock: st.bedrock, regolith: st.regolith, mineral: st.mineral, humus: st.humus, litter: st.litter,
soilDepth: derived(st, 'soilDepth'), // mineral + humus (the true soil body)
surface: derived(st, 'surface'), // bedrock + full column (the topographic surface)
sand: st.sand, silt: st.silt, clay: st.clay,
availN: st.availN, availP: st.availP, organicN: st.organicN, organicP: st.organicP,
fertility: derived(st, 'fertility'), // a coupling-friendly 0..1 scalar FROM the real pools
moisture: st.moistIdx, whc: st.whc, water: st.water,
slope: st.slope, flow: st.flow,
initSlope: st.initSlope, // the inherent terrain steepness (undistorted reference)
biomass: st.bio,
ledger: st.ledger
};
},
events(st) { return st.evs; }
};
// ==================================================================================================
// helpers (closures over N and the config)
// ==================================================================================================
// texture accounting: add `m` of new mineral grains with the given sand/silt/clay split
function addTexture(st, k, m, fs, fi, fc) {
const s = fs + fi + fc;
st.sand[k] += m * fs / s; st.silt[k] += m * fi / s; st.clay[k] += m * fc / s;
}
// recompute slope + D8 flow from the current SURFACE (bedrock + full soil column)
function computeSlopeFlow(st) {
const n2 = N * N, H = new Float32Array(n2);
for (let k = 0; k < n2; k++) H[k] = st.bedrock[k] + st.regolith[k] + st.mineral[k] + st.humus[k] + st.litter[k];
st._surf = H;
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) {
const il = Math.max(0, i - 1), ir = Math.min(N - 1, i + 1), jd = Math.max(0, j - 1), ju = Math.min(N - 1, j + 1);
const dx = (H[j * N + ir] - H[j * N + il]) * 0.5 * N, dz = (H[ju * N + i] - H[jd * N + i]) * 0.5 * N;
st.slope[j * N + i] = Math.hypot(dx, dz);
}
// D8 flow accumulation on the surface (high-to-low), each cell routes to its lowest neighbour
const order = st._scratch.order;
for (let k = 0; k < n2; k++) order[k] = k;
// insertion into an array + sort by height desc
const arr = Array.prototype.slice.call(order);
arr.sort((a, b) => H[b] - H[a]);
const flow = st.flow; for (let k = 0; k < n2; k++) flow[k] = 1;
st._downhill = st._downhill || new Int32Array(n2);
for (const k of arr) {
const i = k % N, j = (k / N) | 0; let lo = k, lh = H[k];
for (const d of [[1, 0], [-1, 0], [0, 1], [0, -1], [1, 1], [1, -1], [-1, 1], [-1, -1]]) {
const ni = i + d[0], nj = j + d[1]; if (ni < 0 || ni >= N || nj < 0 || nj >= N) continue;
const m = nj * N + ni; if (H[m] < lh) { lh = H[m]; lo = m; }
}
st._downhill[k] = lo === k ? -1 : lo; // -1 = a pit/outlet
if (lo !== k) flow[lo] += flow[k];
}
st._order = arr;
}
// EROSION + DEPOSITION over the D8 network, routing detached topsoil downhill and settling it.
function erodeAndDeposit(st, dy) {
const n2 = N * N, sc = st._scratch;
const { humus, mineral, litter, sand, silt, clay, availN, availP, organicN, organicP, land, bio } = st;
const mob = sc.mob, mobH = sc.mobH, mobN = sc.mobN, mobP = sc.mobP,
mobSand = sc.mobSand, mobSilt = sc.mobSilt, mobClay = sc.mobClay;
for (let k = 0; k < n2; k++) { mob[k] = 0; mobH[k] = 0; mobN[k] = 0; mobP[k] = 0; mobSand[k] = 0; mobSilt[k] = 0; mobClay[k] = 0; }
// 1) DETACH: stream power strips the topsoil, vegetation cover protects it
let totalDetached = 0;
for (let k = 0; k < n2; k++) {
if (!land[k]) continue;
const cover = clamp(bio[k] / 0.5, 0, 1) * coverProtect; // dense plants armour the ground
const stream = erodeK * Math.pow(st.flow[k], flowM) * Math.pow(st.slope[k], slopeN) * (1 - cover) * dy;
if (stream <= 0) continue;
// detach from the top down: litter first, then humus, then FINE mineral (silt+clay), sparing sand
const topsoil = litter[k] + humus[k] + silt[k] + clay[k];
const det = Math.min(stream, topsoil * 0.5); // never strip more than half the topsoil in a step
if (det <= 0) continue;
// apportion the detachment across the top layers by availability
let rem = det;
const takeLit = Math.min(litter[k], rem); litter[k] -= takeLit; rem -= takeLit; mob[k] += takeLit;
const takeHum = Math.min(humus[k], rem); humus[k] -= takeHum; rem -= takeHum; mobH[k] += takeHum; mob[k] += takeHum;
// fine mineral (silt+clay) rides next, in proportion
const fine = silt[k] + clay[k];
if (fine > 1e-9 && rem > 0) {
const takeFine = Math.min(fine, rem);
const fSilt = takeFine * silt[k] / fine, fClay = takeFine * clay[k] / fine;
silt[k] -= fSilt; clay[k] -= fClay; mineral[k] -= takeFine;
mobSilt[k] += fSilt; mobClay[k] += fClay; mob[k] += takeFine; rem -= takeFine;
}
// nutrients ride with the eroded humus + fine mineral (they are adsorbed to them)
const stripFrac = (takeHum + (det - takeLit - takeHum)) / Math.max(humus[k] + takeHum + silt[k] + clay[k] + (det - takeLit - takeHum), 1e-9);
const carryN = availN[k] * clamp(stripFrac, 0, 0.9) + organicN[k] * clamp(takeHum / Math.max(humus[k] + takeHum, 1e-9), 0, 0.9);
const carryP = availP[k] * clamp(stripFrac, 0, 0.9) + organicP[k] * clamp(takeHum / Math.max(humus[k] + takeHum, 1e-9), 0, 0.9);
const cN = Math.min(availN[k] + organicN[k], carryN), cP = Math.min(availP[k] + organicP[k], carryP);
// pull carried nutrients proportionally from avail then organic
pullNutrient(st, k, 'availN', 'organicN', cN); mobN[k] += cN;
pullNutrient(st, k, 'availP', 'organicP', cP); mobP[k] += cP;
totalDetached += mob[k];
st.cumEroded[k] += mob[k]; // this cell's cumulative soil loss
}
// 2) ROUTE downhill in high-to-low order. A cell's detachment leaves it entirely (it never re-settles
// at its own origin); the load that ARRIVES from upslope settles a fraction here (deposition) and
// carries the rest onward. So a steep, high-flow cell is a net LOSER (it detaches and passes
// through) and a flat valley cell is a net GAINER (much arrives, little leaves, most settles).
// Whatever reaches the single low OUTLET notch leaves the world and is ledgered as export.
const order = st._order, downhill = st._downhill;
// `inb*` = load that has arrived at a cell from upslope and is available to settle or pass on.
const inb = sc.mob2 || (sc.mob2 = new Float32Array(n2)), inbH = sc.mobH2 || (sc.mobH2 = new Float32Array(n2)),
inbN = sc.mobN2 || (sc.mobN2 = new Float32Array(n2)), inbP = sc.mobP2 || (sc.mobP2 = new Float32Array(n2)),
inbSa = sc.mobSa2 || (sc.mobSa2 = new Float32Array(n2)), inbSi = sc.mobSi2 || (sc.mobSi2 = new Float32Array(n2)),
inbCl = sc.mobCl2 || (sc.mobCl2 = new Float32Array(n2));
for (let k = 0; k < n2; k++) { inb[k] = mob[k]; inbH[k] = mobH[k]; inbN[k] = mobN[k]; inbP[k] = mobP[k]; inbSa[k] = mobSand[k]; inbSi[k] = mobSilt[k]; inbCl[k] = mobClay[k]; }
let exported = 0;
for (const k of order) {
const dest = downhill[k];
// deposition: a cell settles settleFrac of the load available to it, EXCEPT the share it detached
// itself this step (that mass is leaving downslope, not settling at its source).
const self = mob[k]; // this cell's own detachment (already inside inb[k])
const throughFrac = inb[k] > 1e-12 ? clamp((inb[k] - self) / inb[k], 0, 1) : 0; // fraction that came from upslope
const settle = settleFrac * throughFrac; // only through-flow settles here
const sMob = inb[k] * settle, sH = inbH[k] * settle, sN = inbN[k] * settle, sP = inbP[k] * settle,
sSa = inbSa[k] * settle, sSi = inbSi[k] * settle, sCl = inbCl[k] * settle;
if (sMob > 0 || sH > 0) { depositHere(st, k, sMob, sH, sN, sP, sSa, sSi, sCl); st.cumDeposited[k] += sMob; }
const oMob = inb[k] - sMob, oH = inbH[k] - sH, oN = inbN[k] - sN, oP = inbP[k] - sP,
oSa = inbSa[k] - sSa, oSi = inbSi[k] - sSi, oCl = inbCl[k] - sCl;
if (dest < 0) {
// the outlet notch (or an interior pit): interior pits keep the mass (settle it all), the true
// low edge-outlet exports it out of the world.
const i = k % N, j = (k / N) | 0;
const edge = (i === 0 || j === 0 || i === N - 1 || j === N - 1);
if (edge) {
exported += oMob; st.ledger.erodedOut += oMob;
st.ledger.erodedN += oN; st.ledger.erodedP += oP; // nutrients riding eroded soil off the grid
} else {
depositHere(st, k, oMob, oH, oN, oP, oSa, oSi, oCl); st.cumDeposited[k] += oMob;
}
} else {
inb[dest] += oMob; inbH[dest] += oH; inbN[dest] += oN; inbP[dest] += oP;
inbSa[dest] += oSa; inbSi[dest] += oSi; inbCl[dest] += oCl;
}
}
st._erodedThisStep = totalDetached;
st._exportedThisStep = exported;
}
// deposit a parcel of eroded material into cell k (add mass back into the layers it came from)
function depositHere(st, k, mass, hum, nN, nP, sa, si, cl) {
if (mass <= 0 && hum <= 0) return;
st.humus[k] += hum;
st.litter[k] += Math.max(0, mass - hum - sa - si - cl); // whatever was litter
st.sand[k] += sa; st.silt[k] += si; st.clay[k] += cl;
st.mineral[k] += sa + si + cl;
// deposited nutrients: land in available + organic (mostly organic, adsorbed to the humus/clay)
st.availN[k] += nN * 0.5; st.organicN[k] += nN * 0.5;
st.availP[k] += nP * 0.5; st.organicP[k] += nP * 0.5;
}
// pull `amount` of a nutrient from a cell, drawing available first then organic
function pullNutrient(st, k, availKey, orgKey, amount) {
let rem = amount;
const a = Math.min(st[availKey][k], rem); st[availKey][k] -= a; rem -= a;
const o = Math.min(st[orgKey][k], rem); st[orgKey][k] -= o; rem -= o;
}
// ---- driver ingestion ----------------------------------------------------------------------------
function ingestTerrain(st, ctx) {
const tf = ctx && ctx.fields ? ctx.fields[terrainName] : null;
if (tf && tf.elevation && !st._haveTerrain) {
// adopt the real terrain's bedrock (resampled to our N) once
const TN = tf.N || N;
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) {
const si = Math.round(i / (N - 1) * (TN - 1)), sj = Math.round(j / (N - 1) * (TN - 1));
st.bedrock[j * N + i] = tf.elevation[sj * TN + si];
st.land[j * N + i] = tf.land ? tf.land[sj * TN + si] : 1;
}
st._haveTerrain = true;
// re-baseline the inherent-steepness reference to the coupled terrain
st.bedrock0 = st.bedrock.slice();
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) {
const il = Math.max(0, i - 1), ir = Math.min(N - 1, i + 1), jd = Math.max(0, j - 1), ju = Math.min(N - 1, j + 1);
const dx = (st.bedrock0[j * N + ir] - st.bedrock0[j * N + il]) * 0.5 * N;
const dz = (st.bedrock0[ju * N + i] - st.bedrock0[jd * N + i]) * 0.5 * N;
st.initSlope[j * N + i] = Math.hypot(dx, dz);
}
computeSlopeFlow(st);
}
}
// build this step's climate fields (temp, moisture, freeze-thaw intensity, warmth) either from the
// coupled weather system or from the built-in climate model. Deterministic; uses ctx.clock for season.
//
// Pedogenesis steps on sim-YEARS, so the built-in climate models the ANNUAL climate, not an instant:
// - temp is the annual MEAN (with a small seasonal wobble that only matters if a caller steps
// sub-annually), cooled with elevation by a lapse rate;
// - WARMTH for chemistry/biology is the GROWING-SEASON average -- the year's warm half does the
// ripening even where the annual mean is cold, so a montane soil still develops, only slower;
// - FREEZE-THAW is driven by whether the year's temperature RANGE (mean +/- seasonal +/- diurnal
// amplitude) STRADDLES 0C. A cell always frozen (deep cold) or always thawed (warm lowland) gets
// few cycles; a cell whose swing crosses freezing shatters rock hardest. This is the correct
// annual-timestep form of freeze-thaw and does not depend on WHERE in the year a step lands.
function sampleClimate(st, ctx, dt) {
const n2 = N * N;
const moist = new Float32Array(n2), phys = new Float32Array(n2), warm = new Float32Array(n2);
const wf = ctx && ctx.fields ? ctx.fields[weatherName] : null;
const clock = ctx && ctx.clock ? ctx.clock : null;
const dyStep = (dt || 0) / 31536000;
// seasonal wobble only bites for sub-annual steps; for year+ steps it averages out to ~0
const seasonT = dyStep < 0.5 && clock ? Math.sin((clock.yearPhase() - 0.25) * 2 * Math.PI) * (1 - clamp(dyStep * 2, 0, 1)) : 0;
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) {
const k = j * N + i;
let Tmean, m;
if (wf && wf.temp) {
const WN = wf.N || N; const si = Math.round(i / (N - 1) * (WN - 1)), sj = Math.round(j / (N - 1) * (WN - 1));
Tmean = wf.temp[sj * WN + si];
m = wf.moisture ? clamp(wf.moisture[sj * WN + si], 0, 1) : (wf.rain ? clamp(wf.rain[sj * WN + si], 0, 1) : moistSupply);
} else {
// built-in climate: annual mean, cooler with elevation (bedrock height as a proxy for altitude)
const elev = st.bedrock[k];
Tmean = meanTemp + tempSwing * seasonT - elev * 6; // lapse rate: high ground is colder
m = moistSupply * (0.6 + 0.4 * (1 - elev)); // valleys a touch wetter
}
st.temp[k] = Tmean;
moist[k] = clamp(m, 0, 1);
// FREEZE-THAW: the annual range is [Tmean - (seasonal+diurnal), Tmean + (seasonal+diurnal)]. The
// drive peaks when that range is CENTRED on 0C (equal freezing and thawing) and falls off as the
// mean pulls fully above or below freezing.
const amp = tempSwing + diurnalSwing;
let ftl = 0;
if (Tmean - amp < 0 && Tmean + amp > 0) { // the year crosses freezing at all
ftl = smooth(1 - Math.abs(Tmean) / amp); // strongest at the 0C isotherm
}
phys[k] = clamp(ftl, 0, 1);
// WARMTH for chemistry/biology: the GROWING-SEASON average, not the annual mean. Approximate the
// mean of the warm half of the year (mean + a fraction of the seasonal amplitude), ramped 0..1.
const Twarmhalf = Tmean + 0.4 * tempSwing;
warm[k] = clamp(Twarmhalf / 22, 0, 1);
}
return { moist, phys, warm };
}
// vegetation feedback ingestion: if a real veg field is present use its biomass and derive growth/death
// from its change; otherwise the built-in driver (in step) runs. Sets st.bio and st._vegGrowth/_vegDeath.
function ingestVegetation(st, ctx, dy, clim) {
const vf = ctx && ctx.fields ? ctx.fields[vegName] : null;
if (vf && vf.biomass) {
const VN = vf.N || N, n2 = N * N;
st._vegGrowth = st._vegGrowth || new Float32Array(n2);
st._vegDeath = st._vegDeath || new Float32Array(n2);
st._prevVeg = st._prevVeg || new Float32Array(n2);
for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) {
const k = j * N + i;
const si = Math.round(i / (N - 1) * (VN - 1)), sj = Math.round(j / (N - 1) * (VN - 1));
const b = vf.biomass[sj * VN + si];
const prev = st._prevVeg[k];
const dB = b - prev;
st._vegGrowth[k] = Math.max(0, dB) + bioRespire * b * dy; // gross growth ~ net change + turnover
st._vegDeath[k] = bioRespire * b * dy + Math.max(0, -dB);
st.bio[k] = b; st._prevVeg[k] = b;
}
st._haveVeg = true;
} else {
st._haveVeg = false;
}
}
}
// derived fields computed on demand from state (kept out of the hot path)
function derived(st, which) {
const n2 = st.N * st.N, out = new Float32Array(n2);
if (which === 'soilDepth') { for (let k = 0; k < n2; k++) out[k] = st.mineral[k] + st.humus[k]; }
else if (which === 'surface') { for (let k = 0; k < n2; k++) out[k] = st.bedrock[k] + st.regolith[k] + st.mineral[k] + st.humus[k] + st.litter[k]; }
else if (which === 'fertility') {
// a coupling-friendly 0..1 fertility FROM the real pools: depth, moisture, and the scarcer nutrient
for (let k = 0; k < n2; k++) {
const depth = st.mineral[k] + st.humus[k];
const gDepth = depth / (depth + 0.25);
const gN = st.availN[k] / (st.availN[k] + 0.3);
const gP = st.availP[k] / (st.availP[k] + 0.04);
const gM = st.moistIdx[k];
out[k] = clamp(gDepth * Math.min(gN, gP) * (0.4 + 0.6 * gM), 0, 1);
}
}
return out;
}
return { makeSoil };
});