Numerical core — siim._core

The fastscape-free numerical kernels shared by the 1D profile model and the 2D landscape model. Most use NumPy/Numba; spectral flexure uses SciPy. They need no fastscape, xsimlab, or Matplotlib, so the numerics remain testable without the optional adapter or plotting layer. The 1D and 2D models build one GlacialParams record and an integer law_code and call the step skeletons directly.

Per-run parameter record — siim._core.params

Frozen per-run parameter record for the law_code skeleton kernels.

One GlacialParams bundles every per-run physical/law scalar the skeletons consume, so a kernel call passes one record instead of ~16 positional floats (and no per-law 0.0 padding). Built once per run by the model, unpacked at the top of each skeleton; the per-law dispatchers still receive plain scalars.

IMPORTANT — this type MUST be defined here (a stable, importable module) and imported everywhere; never reconstruct an equivalent namedtuple in __main__ or per-call. numba keys a @njit(cache=True) function’s on-disk cache on its argument types, and a namedtuple’s type identity comes from its defining module + qualname. A __main__- or locally-defined twin is a different numba type, so the cache silently misses and recompiles every process (same failure class as the closure/exec patterns the Step-0 spike ruled out; verified directly with NUMBA_DEBUG_CACHE).

class siim._core.params.GlacialParams(Ko, Co, ce, n, nu, m, mu, cg, alpha_g, lambda_p, lambda_c, tau_c, coulomb_clamp, rho_g_g, hc_over_H, D_H)

Bases: tuple

Order is load-bearing: the skeletons tuple-unpack p positionally. Every field defaults to 0.0 so the production builders keyword-construct only their law’s real fields (the inactive-law padding stays implicit; audit N19).

Co

Alias for field number 1

D_H

Alias for field number 15

Ko

Alias for field number 0

alpha_g

Alias for field number 8

ce

Alias for field number 2

cg

Alias for field number 7

coulomb_clamp

Alias for field number 12

hc_over_H

Alias for field number 14

lambda_c

Alias for field number 10

lambda_p

Alias for field number 9

m

Alias for field number 5

mu

Alias for field number 6

n

Alias for field number 3

nu

Alias for field number 4

rho_g_g

Alias for field number 13

tau_c

Alias for field number 11

Step skeletons (law_code switch) — siim._core.skeleton

Mode-A/B step skeletons (law_code switch) for the siim numerical core.

One @njit(cache=True) skeleton per (mode x routing) for 2D, plus the 1D joint-walk skeleton. Each takes a small-int law_code and reaches the per-law physics through the shared dispatchers in siim._core.solvers / siim._core.eroders. The 2D model (siim.fastscape.processes) and the 1D model (siim.siim1d) call these skeletons directly with a GlacialParams record + an integer law_code. numpy/numba only – no model/fastscape imports, so the numerical core stays importable without the fastscape stack.

siim._core.skeleton._flot_factor(zs, H, bl, hc_over_H, ramp)

Waterline-flotation factor f in [0, 1] multiplying glacial erosion (the effective-pressure ramp described in docs/guides/concepts.md):

f = clip((zs - bl) / (ramp*hc_over_H*H), 0, 1)

with f = 0 EXACTLY for zs <= bl (the anti-runaway backstop: a fully afloat column does no glacial erosion and, via E_c <= 0, carves nothing). ramp = gamma, the ramp width in ice-column heights (delta = gamma*hc*H). ramp <= 0 — or a degenerate delta (thin ice) — falls back to the hard binary gate: f = 1 for zs >= bl, 0 below, bit-for-bit the pre-ramp behavior. Callers apply it only where the gate is on and the cell is icy.

siim._core.skeleton._implicit_border_step(zb0, U, E, dt, H, hc_over_H, bl, ramp)

Closed-form backward-Euler step of the border-bed budget dzb/dt = U - f(zb)*E, with f the flotation ramp evaluated at the NEW bed:

g(z) = z - zb0 - (U - f(z)*E)*dt
f(z) = clip((z + hc*H - bl) / (ramp*hc*H), 0, 1)

E (the arrival-slope border erosion rate), H and the slope are frozen per step, so f is piecewise linear in z and g is strictly increasing (g’ = 1 + f’E*dt >= 1): exactly one of three branches is consistent — fully grounded (f = 1), fully afloat (f = 0), or the linear ramp. The iterate approaches the flotation-draft equilibrium zb = bl - hc*H + delta*U/E (delta = ramp*hc*H) MONOTONICALLY — it cannot overshoot at any dt (the reason the border budget is dt-robust where the explicit form dug km/step). ramp <= 0 is the binary gate, whose implicit solution is the Filippov sliding mode: the bed sticks at the flotation manifold z = bl - hc*H (no chatter possible). Verified closed-form-exact: residual < 2e-11 over 20 000 random draws (probe + build-time check).

siim._core.skeleton._lake_fill_sfr_2d(z_flat, stack, rec)

In-place monotone fill: walk stack outlet-first; any cell below its receiver gets raised to receiver level. The 2D analog of lake_fill_1d — relies on the basin-corrected receivers from fastscape’s flowrouting (lake-interior cells point toward the spillway), so this single pass fills closed depressions to the spillway elevation. Boundary cells (rec == self) are untouched.

siim._core.skeleton._lake_fill_1d(zb, didx_l, didx_r, nx)

In-place monotone fill on the 1D bed: walk each flank from the divide outward and raise any node below its downstream (already-filled) receiver to it — the 1D analog of _lake_fill_sfr_2d(). Every closed basin spills at its downstream rim. zb is MODIFIED IN PLACE. njit’d because it is on the default mode-B 1D hot path (~33% of the step; audit m55).

siim._core.skeleton._glac_fast_solve_modeA_sfr(z_flat, ice_flux, water_flux, H_flat, law_code, p, dt, lengths, stack, rec)

SFR + mode A (z-tracking), law-agnostic skeleton.

One @njit(cache=True) kernel for all three sliding laws (see the mode-B skeleton _glac_fast_solve_modeB_sfr() for the dispatch design). Mode A tracks the ice surface z_flat directly, so H is a per-node local solve from the surface slope S = (z_i - z_r)/L — exactly the from_slope branch of _modeb_closure() (hc_over_H is irrelevant when the surface is tracked, so 1.0 is passed) — and there is no lake-fill, border-bed budget, surface_out or hc_over_H. The law enters only at the H-closure and the erosion step, both dispatched on law_code; the thin wrappers pass their LAW_* code and 0.0 for the inactive-law constants (which reach only the unused branch). z_flat and H_flat are MODIFIED IN PLACE.

Parameters:

• z_flat (ndarray) – Ice surface elevation, flattened — MODIFIED IN PLACE.

• ice_flux (ndarray) – Accumulated ice / water flux (m^3/yr), flattened.

• water_flux (ndarray) – Accumulated ice / water flux (m^3/yr), flattened.

• H_flat (ndarray) – Ice thickness (m), flattened — MODIFIED IN PLACE.

• law_code (int) – LAW_EFFEXP / LAW_POWER / LAW_COULOMB (see _core.solvers).

• p (GlacialParams) – Packed law/physics constants (see siim._core.params). Notable field: cg = alpha_g * kt * (2*Ac/5) * (rho_g*g)^3 [m^-3 yr^-1] (kt absorbed).

• lengths (ndarray) – Distance from each node to its receiver (m).

• stack (ndarray) – Topological sort (upstream → downstream).

• rec (ndarray) – Receiver indices.

siim._core.skeleton._glac_fast_solve_modeA_dinf(z_flat, ice_flux, water_flux, H_flat, law_code, p, dt, stack, nb_receivers, receivers, weights, lengths)

D-inf + mode A (z-tracking), law-agnostic skeleton.

The D-inf twin of _glac_fast_solve_modeA_sfr(): H is solved per node from the weighted-mean per-cell slope S = Σ_k w_k·max(0, (z_i-z_rk)/L_k) over the cell’s D-inf receivers — that single effective slope feeds the same from_slope branch of _modeb_closure() (hc_over_H irrelevant, 1.0 passed). Erosion uses the D-inf eroders via _erode_modeb_dinf() on (z_flat, zo). z_flat and H_flat are MODIFIED IN PLACE.

siim._core.skeleton._glac_fast_solve_modeB_sfr(zb_flat, ice_flux, water_flux, H_flat, surface_out, law_code, p, dt, lengths, stack, rec, ny, nx, dx_cell, dy_cell, border_bed_uplift, bl=0.0, gate=True, ramp=0.1, wrap_y=False, wrap_x=False, parallel_erode=False)

SFR + mode B (zb-tracking joint walk + lake-fill), law-agnostic skeleton.

One @njit(cache=True) kernel for all three sliding laws; the per-law physics is reached through the law_code dispatch at the two sites that vary — the joint-walk H-closure (_modeb_closure()) and the erosion step (_erode_modeb_sfr()). Everything else (the outflow base-level ice BC, the two-view waterline lake-fill and the interior flotation ramp) is law-agnostic. The thin public wrappers pass their law’s LAW_* code and 0.0 for the inactive-law constants, which reach only the unused dispatch branch.

Base-level ice border = OUTFLOW: the domain edge is an arbitrary cut through a continuing glacier. A through-flowing border gets zero-gradient thickness H_border = H_dominant_donor (the max-ice-flux interior donor), and its bed keeps eroding by the IMPLICIT BORDER BUDGET dzb/dt = U - f*E — E the glacial law on the interior ARRIVAL slope, f the flotation ramp — integrated by the closed-form backward-Euler step _implicit_border_step(), so the bed approaches the flotation-draft equilibrium zb* = bl - hc*H + delta*U/E monotonically at any dt. The waterline-flotation gate (gate/ramp) is the rho_i = rho_w effective-pressure law, applied to the interior erosion delta via _flot_factor() and inside the border step (ramp = 0 is the hard binary gate / its Filippov sliding mode).

Parameters:

• zb_flat (ndarray) – Bedrock, flat — uplift already applied; MODIFIED IN PLACE.

• ice_flux (ndarray) – Accumulated ice / water flux (m^3/yr), flattened. The interior flotation ramp, H-closure, erosion, the dominant-donor pick and the border icy/ice-free switch all read ice_flux.

• water_flux (ndarray) – Accumulated ice / water flux (m^3/yr), flattened. The interior flotation ramp, H-closure, erosion, the dominant-donor pick and the border icy/ice-free switch all read ice_flux.

• H_flat (ndarray) – Ice thickness; previous-step value on input, MODIFIED IN PLACE.

• surface_out (ndarray) – Filled with zb + hc_over_H*H_new (the new ice surface) on return.

• law_code (int) – LAW_EFFEXP / LAW_POWER / LAW_COULOMB (see _core.solvers).

• p (GlacialParams) – Packed law/physics constants (see siim._core.params). Notable fields: D_H (H diffusivity [m^2/yr]; 0 disables, sub-stepped CFL) and hc_over_H (centerline-to-mean depth ratio; zs = zb + hc_over_H*H).

• border_bed_uplift (ndarray) – Flat (nn,) uplift rate (m/yr); only border (self-receiving) cells are read — the U in the icy border budget and the rate of the ice-free post-glacial recovery toward bl.

• bl (float or ndarray) – Base level: the per-step water-line (Dirichlet) datum. Replaces the literal 0 at every waterline site — the ice-free-border erosion-view floor + lake-fill seed, the border recovery threshold, and the flotation reference. A scalar is the single domain-wide datum (default 0.0, bit-for-bit with the historical hard-coded datum); a flat (nn,) array carries each base-level BORDER node’s own side datum (per-side bl), and every interior node then inherits the datum of the outlet its basin drains to (step 0a).

• gate (bool) – Waterline-flotation gate (default constants.FLOTATION_GATE = True): the rho_i = rho_w effective-pressure law. Interior: scales the erosion delta. Border: the physical bound inside the implicit budget. Off is for diagnostics only — it un-bounds the border (f == 1, measured runaway).

• ramp (float) – Flotation-ramp width gamma (default constants.FLOTATION_RAMP = 0.1): glacial erosion is scaled by f = clip((zs - bl)/(gamma*hc*H), 0, 1) instead of the hard on/off switch. 0 = the hard binary gate (interior bit-for-bit; at the border its implicit solution is the flotation sliding mode). Safe ceiling 0.2 (see constants.FLOTATION_RAMP). Only active when gate is on.

• parallel_erode (bool) – Run the erosion step level-scheduled in parallel (topological levels of the flow graph; _core.routing._levels_sfr + eroders._erode_modeb_sfr_levels). BIT-FOR-BIT with the serial eroder at any thread count (disjoint writes, identical per-node arithmetic; pinned by test_parallel_erode). Default False.

siim._core.skeleton._dinf_modeB_recv(zb_flat, H_flat, inode, nb_receivers, receivers, weights, lengths, hc_over_H)

Effective receiver for the joint walk: a = zb_i − Z̄ and L̄ from the weighted receiver surfaces (receivers already solved — stack order). Returns (a, L_eff, ok).

siim._core.skeleton._dinf_modeB_filled_view(zb_flat, H_flat, receivers, hc_over_H, ny, nx, wrap_y, wrap_x, bl_node)

Erosion working view: z’ = zb + hc_over_H*H with ICE-FREE borders floored at their own water line bl_node[i], depression-filled by the flat (eps = 0) priority flood seeded at the borders — the 2D generalization of the SFR lake-fill stack walk. An ICY (outflow) border keeps its true surface (a free outflow, not still base water). Returns the filled view.

siim._core.skeleton._glac_fast_solve_modeB_dinf(zb_flat, ice_flux, water_flux, H_flat, surface_out, law_code, p, dt, stack, nb_receivers, receivers, weights, lengths, ny, nx, dx_cell, dy_cell, border_bed_uplift, wrap_y, wrap_x, bl=0.0, gate=True, ramp=0.1, parallel_erode=False)

D-inf + mode B, law-agnostic skeleton.

The D-inf twin of _glac_fast_solve_modeB_sfr() and the routing twin of _glac_fast_solve_modeA_dinf(): same six-step structure (joint walk, diffuse, filled view, erode, implicit border budget, output surface) over the D-inf graph. The joint walk visits the donor-first stack in REVERSE (receivers-first) and collapses the weighted receiver surfaces to a single effective receiver (_dinf_modeB_recv(), a = zb_i - Z̄ and L̄), so the per-law H-closure is the SAME _modeb_closure() used by SFR; erosion is the D-inf dispatch _erode_modeb_dinf(). See the D-inf mode-B block comment above. Thin wrappers pass their LAW_* code and 0.0 for the inactive-law constants. zb_flat / H_flat are MODIFIED IN PLACE and surface_out is filled.

OUTFLOW base-level ice BC (see the SFR twin): a through-flowing border gets zero-gradient thickness H_border = H_dominant_donor (max ice flux into it) and its bed erodes by the IMPLICIT BORDER BUDGET on the interior arrival slope (_implicit_border_step(), ramp-bounded at the flotation draft). The flotation gate/ramp is the rho_i = rho_w effective-pressure law (ramp = gamma; 0 = the hard binary gate / sliding mode). bl is the per-step water-line datum — a scalar (default 0.0 = bit-for-bit) or a flat (nn,) per-side border datum resolved to the interior by the dominant-receiver basin-outlet walk (see the SFR twin). parallel_erode: level-scheduled parallel erosion step, bit-for-bit with the serial eroder (see the SFR twin). Default False.

siim._core.skeleton._diag_walk(zb, Qg, H_out, law_code, p, dx, didx_l, didx_r, nx)

Mode-B 1D joint walk (H + surface), law-agnostic skeleton.

The 1D twin of the 2D mode-B skeletons. Walk both sides of the divide from the base-level outlets inward, solving H and z_s = zb + hc_over_H*H jointly per node. The per-law H-closure is the SAME _modeb_closure the 2D kernels use — the from-a branch with cell length L = dx at interior nodes. A through-flowing base-level outlet is an OUTFLOW border with zero-gradient thickness H_outlet = H_interior_neighbour, resolved by a 2-pass walk (the neighbour is solved after the outlet, so pass 1 uses the lagged outlet H as a bounded provisional and pass 2 the corrected one). Qg = 0 keeps H = 0. Only the closure dispatch varies by law; the outflow outlet BC and the Qg<=0 terminus are law-agnostic and live here. H_out is MODIFIED IN PLACE (its incoming value is the previous-step provisional); the thin wrappers pass their LAW_* code and 0.0 for the inactive-law constants.

The 1D border bed (the implicit arrival-slope budget) lives in siim.siim1d._erode_border_bed_1d — this walk solves only H.

Scalar Newton solvers & ice-thickness closures — siim._core.solvers

Shared scalar Newton solvers and ice-thickness closures for the siim numerical core.

Single canonical home for the per-law scalar kernels that the 1D profile model (siim.siim1d) and the 2D step skeletons (siim._core.skeleton) both consume. These were previously maintained as duplicate copies in the 1D and 2D modules; they are verified identical in executable logic and unified here so a fix lands once. numpy/numba only – no model, fastscape or matplotlib imports, so the numerical core stays importable without the fastscape stack.

All public entry points are @numba.njit(cache=True) and topology-agnostic: scalar (or small fixed-size) inputs, called per-node by the skeleton kernels. The _*_func / _dz*_func helpers are the per-law residuals and their derivatives used by the bracketed Newton iterations.

siim._core.solvers._solve_ice_thickness_power_analytical(D, lambda_p)

Closed-form root of x^3 + lambda_p^2 x^2 - lambda_p^2 D^4 = 0 with x = H^2. Equivalent to solving H*(1 + (H/lambda_p)^2)^(1/4) = D, but with no iteration. Picks the unique positive real root via Cardano (single-root regime) or trig form (three-real-roots regime).

siim._core.solvers._solve_ice_thickness_coulomb(D, a, lambda_c, clamp, tol=1e-12, max_iter=100)

Solve H*(H + lambda_c/(1 - a*H^3))^(1/5) = D for H in (0, H_max). a = (rho_g*g*S / tau_c)^3, pole at H = H_max = a^(-1/3) (tau = tau_c).

Uses Newton + line-search bisection with a residual-based exit criterion (|f| < tol * D). A step-based criterion is unsafe here because near the pole fp diverges, so dH -> 0 even when f is still large, and Newton will exit prematurely with a spurious H. Residual-based guarantees we actually satisfy the equation.

clamp: minimum relative gap from the pole (1e-12 default) that the line search maintains. Prevents 1-a*H^3 underflowing to 0 in double precision.

siim._core.solvers._solver_nonlinear_dinf(zo_i, n_rec, receivers_i, weights_i, lengths_i, z_flat, A, p, dt, epsilon=0.001, max_iter=50)

Multi-receiver D-inf counterpart of _solver_fluvial: one bracketed implicit erosion Newton summed over a cell’s receivers.

F(z) = z - zo_i + sum_k w_k * (dt / L_k^p) * A * max(0, z - z_rk)^p

F is strictly increasing with its root in [min_k z_rk, zo_i]; the step bisects whenever Newton leaves the bracket. Plain (unbracketed) Newton 2-cycles for p < 1 and silently returns zo_i (zero erosion) – the same failure the scalar _solver_fluvial guards against, one receiver-dimension up. Shared by all three *_erode_2d_dinf eroders: (A, p) = (Ko*Qf^m, n) for the fluvial branch of every law, and the eff-exp / power glacial branches pass their own (A, p). Returns the updated node elevation z_i.

Allocation-free (audit N31): the D-inf pack carries at most 2 receivers per cell (receivers is (n, 2)), so the per-receiver prefactors live in scalars instead of per-call heap arrays – same arithmetic, same order, bit-for-bit with the array form it replaced.

siim._core.solvers._F_and_dF_coulomb(zik, zio, zr, Qg, A_pre, ell, t, cg, rho_g_g, tau_c, lambda_c, dx, clamp)

Residual F(zik) and analytical total derivative dF/dz for the regularized Coulomb erosion law.

The rheology factor R*(H+R)^(-3/5) is evaluated via the H-eq substitution:

H + R = D^5/H^5   =>   R*(H+R)^(-3/5) = (D/H)^2 * (1 - H^6/D^5)

which bypasses the catastrophic cancellation in 1 - (rho_g*g*H*S/tau_c)^3. The H^6/D^5 intermediate also gives R = H*(1-H^6/D^5)/(H^6/D^5) without recomputing 1-y, so the derivative stays precision-clean.

Derivation: differentiate H^5(H+R) = D^5 implicitly, using dD/dS = -3D/(5S). The log-slope gamma = (S/H) dH/dS comes out as:

gamma = -3(q+1) / (6q + 5 - 2y),   q = H*lambda_c/R^2,  y = a*H^3

with gamma -> -1 at the pole (H -> H_max) and gamma -> -1/2 in the zero-sliding regime. Then dM/dS = (dM/dD)(dD/dS) + (dM/dH)(dH/dS) using the expanded form M = D^2/H^2 - H^4/D^3 for the partials.

siim._core.solvers._solver_glacial_coulomb(zio, zr, Qg, A_pre, ell, t, cg, rho_g_g, tau_c, lambda_c, dx, clamp, epsilon=1e-08, max_iter=50)

Scalar Newton on zik, using an analytical total derivative (see _F_and_dF_coulomb) plus Armijo |F|-decrease backtracking. One H-solve per accepted Newton step; the Armijo call also returns dF, which is carried into the next iteration.

siim._core.solvers._diag_solve_power_H_newton(a, K_p, lambda_p2, tol=1e-12, max_iter=80)

Solve H^4 * (H^2 + lambda_p^2) * (H + a)^3 = K_p for H > max(0, -a). K_p = Q_g * dx^3 / cg, lambda_p2 = lambda_p^2. f is strictly increasing on the admissible interval, so Newton + bisect-on-bad-step converges.

siim._core.solvers._diag_solve_coulomb_H_newton(a, K_c, lambda_c, beta, clamp, tol=1e-12, max_iter=80)

Solve H^5 * (H+a)^3 * (H + lambda_c/(1-phi^3)) = K_c for H in admissible domain. K_c = Q_g * dx^3 / cg, beta = rho_g*g/(tau_c*dx), phi = beta*H*(H+a), pole at phi=1 (tau = tau_c). Domain: H > max(0, -a), phi < 1-clamp.

siim._core.solvers._modeb_closure_effexp(from_slope, x, qi, L, hc_over_H, cg, lambda_p)

Eff-exp H-closure. from_slope: x = surface slope (L unused); else x = a with cell length L.

siim._core.solvers._modeb_closure_power(from_slope, x, qi, L, hc_over_H, cg, lambda_p)

Power H-closure. from_slope: x = surface slope (L unused); else x = a with cell length L.

siim._core.solvers._modeb_closure_coulomb(from_slope, x, qi, L, hc_over_H, cg, lambda_c, tau_c, rho_g_g, clamp)

Coulomb H-closure. from_slope: x = surface slope (L unused); else x = a with cell length L.

siim._core.solvers._modeb_closure(law_code, from_slope, x, qi, L, hc_over_H, cg, lambda_p, lambda_c, tau_c, rho_g_g, clamp)

Dispatch the mode-B joint-walk H-closure on law_code.

The skeleton calls this at the three closure sub-sites (interior, border-normal, border-degenerate); the inactive law’s constants are passed as 0.0 by the thin wrappers and reach only the unused branch.

2D erosion loops — siim._core.eroders

Per-law 2D erosion loops (SFR + D-inf) for the siim numerical core.

Walk the flow graph and apply the implicit erosion step per node, delegating the per-node closures to the shared solvers in siim._core.solvers. Consolidated into the numerical core in the pre-v1.0 rewrite. numpy/numba only – no model/fastscape imports.

siim._core.eroders._linear_erode_2d(z, zo, Qf, Qg, Kf, Kg, m, mu, stack, rec)

Linear (n=nu=1) implicit erosion on a directed graph.

siim._core.eroders._nonlinear_erode_2d(z, zo, Qf, Qg, Kf, Kg, m, mu, n, nu, stack, rec)

Nonlinear implicit erosion on a directed graph.

siim._core.eroders._power_erode_2d(z, zo, Qf, Qg, H, Kf, Kg_prefactor, m, n, t, lambda_p, stack, rec)

Power sliding law erosion on a directed graph.

siim._core.eroders._linear_erode_2d_dinf(z, zo, Qf, Qg, dt, Ko, Co, m, mu, stack, nb_receivers, receivers, weights, lengths)

Linear (n=nu=1) implicit erosion on a multi-flow graph.

For ν=n=1, the implicit equation collapses to:

z (1 + K) = zo + Σ_k w_k · Gi_k · z_rk

where K = Σ_k w_k · Gi_k. Closed-form, no Newton needed.

siim._core.eroders._nonlinear_erode_2d_dinf(z, zo, Qf, Qg, dt, Ko, Co, m, mu, n, nu, stack, nb_receivers, receivers, weights, lengths, epsilon=0.001, max_iter=50)

Multi-receiver implicit Newton on z for the nonlinear case.

F(z) = z - zo + Σ_k w_k · (dt/L_k^p) · A · max(0, z - z_rk)^p dF/dz = 1 + Σ_k w_k · p · (dt/L_k^p) · A · max(0, z - z_rk)^(p-1) where (p, A) = (ν, Co·Q_g^μ) when glacial, else (n, Ko·Q_f^m).

siim._core.eroders._power_erode_2d_dinf(z, zo, Qf, Qg, H, dt, Ko, ce, m, n, t, lambda_p, cg, alpha_g, stack, nb_receivers, receivers, weights, lengths, epsilon=0.001, max_iter=50)

Power sliding law D-inf erosion. Same structure as _nonlinear_erode_2d_dinf but with the per-cell H-dependent G_o prefactor for the glacial branch.

siim._core.eroders._coulomb_erode_2d(z, zo, Qf, Qg, Kf, A_const_nodal, m, n, ell, t, cg, rho_g_g, tau_c, lambda_c, clamp, lengths, stack, rec)

Regularized Coulomb sliding law erosion on a directed graph (with kt absorbed in cg). A_const_nodal[i] = (dt / lengths[i]^t) * ce * (cg^(2/5)/alpha_g)^ell; the per-node Qg^(3*ell/5) and (H, R)-dependent mass factor are applied inside.

siim._core.eroders._F_coulomb_dinf_residual(zik, zo_i, n_rec, receivers_i, weights_i, lengths_i, z_flat, Qg, base_A, ell, t, exp_Q, cg, rho_g_g, tau_c, lambda_c, clamp)

Compute the multi-receiver coulomb residual F(zik) and dF/dz.

Mass factor M and H are computed from the cell’s weighted-mean slope S_eff = Σ_k w_k · max(0, (zik - z_rk)/L_k); the derivative is lagged in S (we drop the dM/dS coupling and rely on Newton line-search for robustness). Returns (F, dF). If no downhill receiver exists, returns (zik-zo, 1).

siim._core.eroders._solver_glacial_coulomb_dinf(zo_i, n_rec, receivers_i, weights_i, lengths_i, z_flat, Qg, base_A_with_dt, ell, t, exp_Q, cg, rho_g_g, tau_c, lambda_c, clamp, epsilon=1e-08, max_iter=50)

Multi-receiver coulomb Newton with Armijo |F|-decrease backtracking.

siim._core.eroders._coulomb_erode_2d_dinf(z, zo, Qf, Qg, dt, Ko, ce, m, n, ell, t, cg, rho_g_g, tau_c, lambda_c, clamp, alpha_g, stack, nb_receivers, receivers, weights, lengths)

Coulomb sliding law erosion on a multi-receiver graph.

Per-cell H is computed from the weighted-mean slope; erosion is solved by multi-receiver Newton with lagged-S Jacobian and line-search backtracking. Fluvial fallback shares the bracketed _solver_nonlinear_dinf with the eff-exp/power D-inf eroders.

siim._core.eroders._modeb_border_erosion(law_code, qi, slope, Hi, Co, mu, nu, ce, cg, alpha_g, lambda_p)

Dispatch the mode-B border-bed glacial erosion rate on law_code.

Inactive-law constants are passed as 0.0 by the thin wrappers and reach only the unused branch.

siim._core.eroders._erode_modeb_sfr(law_code, z_filled, z_pre, water_flux, ice_flux, H_flat, Ko, Co, ce, n, nu, m, mu, cg, alpha_g, lambda_p, lambda_c, tau_c, coulomb_clamp, rho_g_g, dt, lengths, stack, rec)

Dispatch the mode-B SFR erosion step on law_code.

siim._core.eroders._erode_modeb_dinf(law_code, z, zo, water_flux, ice_flux, H, Ko, Co, ce, n, nu, m, mu, cg, alpha_g, lambda_p, lambda_c, tau_c, coulomb_clamp, rho_g_g, dt, stack, nb_receivers, receivers, weights, lengths)

Dispatch the mode-B D-inf erosion step on law_code.

siim._core.eroders._erode_modeb_sfr_levels(law_code, z, zo, water_flux, ice_flux, H_flat, Ko, Co, ce, n, nu, m, mu, cg, alpha_g, lambda_p, lambda_c, tau_c, coulomb_clamp, rho_g_g, dt, lengths, rec, order, offsets, nlev)

Level-parallel twin of _erode_modeb_sfr() (same dispatch, same per-node arithmetic; order/offsets/nlev from siim._core.routing._levels_sfr()).

siim._core.eroders._erode_modeb_dinf_levels(law_code, z, zo, water_flux, ice_flux, H, Ko, Co, ce, n, nu, m, mu, cg, alpha_g, lambda_p, lambda_c, tau_c, coulomb_clamp, rho_g_g, dt, nb_receivers, receivers, weights, lengths, order, offsets, nlev)

Level-parallel twin of _erode_modeb_dinf() (same dispatch, same per-node arithmetic; order/offsets/nlev from siim._core.routing._levels_dinf()).

Sub-grid glacier-width carving — siim._core.carve

Sub-grid glacier-width carving (mode B only).

Consolidated into the numerical core in the pre-v1.0 rewrite. See docs/guides/concepts.md for the public description of sub-grid width carving. numpy/numba only – no model/fastscape imports.

siim._core.carve._carve_offsets(H_flat, rec, alpha_g, offsets_flat, seed_mask)

Seed array for the width-carve power transform: -R^2 with R = alpha_g*H/2 at icy, non-border, seed-allowed cells; PDT_NO_SOURCE elsewhere. Border (self-receiving) cells are left out of the seeds so interior sources inherit their attribution (see _carve_subgrid_width).

seed_mask (per-cell, nonzero = allowed) can de-seed cells that must not anchor a disc; the Mode-C carve passes an all-ones mask (every icy interior cell seeds a disc). Returns the seed count (0 = nothing this step).

siim._core.carve._power_dt_1d(f, h2, d_out, i_out)

1D lower envelope of parabolas: d_out[p] = min_q((p-q)^2*h2 + f[q]), i_out[p] = argmin q. Algorithm 1 of Felzenszwalb and Huttenlocher [FH12] with a physical spacing-squared factor h2 (anisotropic grids) and +inf-aware seeding.

siim._core.carve._power_dt_2d(offsets, dy, dx, D, SRC)

Minimum power distance D[y,x] = min_src(dist^2 - R_src^2) over the grid (physical units, anisotropic spacing), and SRC = the flat index of the argmin source cell (-1 outside any line of sources). offsets holds -R^2 at source cells and >=1e30 elsewhere.

siim._core.carve._power_dt_2d_periodic(offsets, dy, dx, D, SRC, wrap_y, wrap_x)

Periodic-aware power transform: wrap-pad the seed array along the looped axes, run the ordinary FH transform on the padded grid, crop the central window, and remap the argmin labels to original flat indices (numpy wrapper; the njit kernel is untouched).

EXACT for the carve, by two arguments. (1) Only each seed’s NEAREST periodic image matters: farther images of the same seed have the same R and larger d, hence strictly worse power. (2) Footprint membership and in-footprint attribution only involve seeds with d < R <= R_max (a footprint cell’s winner satisfies it, and any competitor that could beat the winner there satisfies it too), so padding by ceil(R_max/spacing) cells covers every image that can win; the half-circumference cap covers the R_max > L/2 regime, where the nearest image always lies within half a wrap. Cells outside every footprint (D >= 0) may in principle miss a remote image, but the carve skips them — and D < 0 vs >= 0 itself is exact by (2).

siim._core.carve._carve_subgrid_width(zb_flat, zb_kern, zb_pre, H_flat, surface_out, rec, D, SRC, offsets, widening_factor, hc_over_H)

Apply the sub-grid width carve (see block comment above) to every footprint cell whose bed stands above its target — bare valley walls AND thin-ice cells inside a bigger glacier’s footprint (“nodes that were in a glacier and didn’t know it”: under broad ice cover nearly the whole footprint is icy, and skipping icy cells starves the carve).

Parameters:

• zb_flat (ndarray) – Bed, MODIFIED IN PLACE (post-kernel on entry).

• zb_kern (ndarray) – Snapshot of zb_flat at carve entry (post-kernel, pre-carve) — source anchors are read here, so the result is independent of cell visit order even when sources are themselves carved.

• zb_pre (ndarray) – Bed before this step’s erosion kernel — descent caps are measured from here, so the kernel’s own erosion at a cell and the carve never add (whichever is lower wins).

• offsets (ndarray) – The -R^2 seed array the power transform consumed — R^2 is read back from it (R2 = -offsets[s]) so footprint definition and carve targets cannot disagree.

• no-op (A self-attributed cell's target is its own bed — a structural)

• as ((the skip is a shortcut); border (self-receiving) cells are skipped)

• sources; (targets (border budget's business) and must not be seeded as)

• carve (sources with no kernel erosion this step (floating / decoupled))

• (bare (nothing. Updates surface_out in place for carved cells)

• icy (bed + hc_over_H*H — the presented surface drops with the bed).)

Explicit H diffusion — siim._core.diffusion

Explicit H diffusion for the siim 2D model (CFL-substepped 5-point FD).

Consolidated into the numerical core in the pre-v1.0 rewrite. numpy/numba only – no model/fastscape imports.

siim._core.diffusion._diffuse_H_2d(H_flat, ny, nx, dx, dy, D, dt, wrap_y=False, wrap_x=False)

Explicit 5-point FD diffusion on H, sub-stepped for CFL: per sub-step ax + ay <= 0.4 with ax = D*dt_sub/dx², ay = D*dt_sub/dy² (the 2D FTCS stability limit is ax + ay <= 0.5). Reduces to the previous square-grid behaviour exactly when dx == dy.

A looped axis (wrap_y / wrap_x) wraps the stencil across the seam so the seam cells diffuse as the interior cells they physically are (matching the fill / facet-scan / carve seam-awareness; audit m16); a non-looped axis holds its outer ring fixed. Bit-for-bit with the old kernel when both are False. No-op when D <= 0. Clamps H >= 0 after each sub-step.

Hillslope diffusion (ADI) — siim._core.hillslope

In-house alternating-direction-implicit hillslope diffuser (topography), the standalone replacement for fastscapelib-fortran’s fs.diffusion. Distinct from siim._core.diffusion (the ice-thickness FD diffuser above).

In-house hillslope diffusion (ADI), the standalone replacement for fastscapelib-fortran’s fs.diffusion (stock LinearDiffusion).

Linear hillslope diffusion \(\partial z/\partial t = \nabla\cdot(k_d\nabla z)\) advanced one step by the alternating-direction-implicit (ADI) scheme, matching fastscapelib-fortran Diffusion.f90 (v2.8.4): two dt/2 half-steps (x-implicit then y-implicit), arithmetic-mean face diffusivities, a Thomas tridiagonal solve per row/column, and the exact Dirichlet (fixed) / one-sided no-flux (free) edge branches keyed to the same ibc digit map. The Thomas solver is independently expressed from the standard row recurrence. The scheme is unconditionally stable and, for spatially-uniform k_d (siim’s scalar D), reproduces fs.diffusion bit-for-bit (twin-gated, rtol ~1e-12).

numpy/numba only – no fastscape/xsimlab imports (framework-free core). This is the HILLSLOPE diffuser (topography); distinct from siim._core.diffusion, the ice-thickness FD diffuser, which is not a migration target.

Edge/ibc map (single source: fs.processes.boundary + Diffusion.f90). cbc = f"{ibc:04d}" and digit '1' = fixed_value (Dirichlet), '0' = free (no-flux):

row j=0      -> cbc[0]      (top)
row j=ny-1   -> cbc[2]      (bottom)
col i=0      -> cbc[3]      (left)
col i=nx-1   -> cbc[1]      (right)

The outer boundary COLUMNS (i=0, i=nx-1) are left unchanged by the scheme (the y-pass writes only interior columns) – a structural quirk of the fortran ADI that is reproduced exactly for the bit-for-bit gate.

siim._core.hillslope._solve_tridiagonal(lower, diagonal, upper, rhs, solution, size)

Solve a tridiagonal system with the standard Thomas recurrence.

lower[i], diagonal[i], and upper[i] multiply the unknowns at indices i-1, i, and i+1. The result is written into the preallocated solution array. This implementation was independently derived from that row equation; it performs no pivoting.

siim._core.hillslope._adi_step(z, kd, dt, dx, dy, fix_r0, fix_rN, fix_c0, fix_cN)

One ADI diffusion step (two dt/2 half-steps) on z (ny, nx), face diffusivities from the arithmetic mean of kd (ny, nx). fix_* are the fixed-value (Dirichlet) flags for the four edges (row 0 / row ny-1 / col 0 / col nx-1); a False flag is the one-sided no-flux branch. Serial (no prange) – grids are small and the tridiagonal solves are sequential; order-independent by construction. Returns the diffused (ny, nx) array. The x-boundary columns are frozen (never written by the y-pass), matching the fortran ADI.

siim._core.hillslope._fixed_flags(ibc)

Decode ibc -> (fix_row0, fix_rowN, fix_col0, fix_colN) booleans.

siim._core.hillslope.diffuse(elevation, diffusivity, dt, nx, ny, xl, yl, ibc)

Diffuse elevation one step (dt) by the ADI scheme; return the diffused surface, same 2D shape (ny, nx). elevation / diffusivity may be flat or 2D; diffusivity is a scalar D or an (ny, nx) field. ibc is the fastscapelib boundary code. Mirrors the stock LinearDiffusion step (fs.diffusion on fs_context['h']), the caller derives erosion = elevation - diffuse(...). dx = xl/(nx-1), dy = yl/(ny-1).

Spectral flexure — siim._core.flexure

In-house scipy.fft thin-plate flexure solve on the native grid, the standalone replacement for fastscapelib-fortran’s fs.flexure.

In-house spectral flexure, the standalone replacement for fastscapelib-fortran’s fs.flexure (fastscape Flexure / siim GlacialFlexure).

Thin elastic plate on an inviscid asthenosphere,

\[D\,\nabla^4 w + \rho_a\,g\,w = q, \qquad D = \frac{E\,T_e^3}{12\,(1-\nu^2)},\]

solved spectrally on the native (ny, nx) grid (no quarter-grid resample) in a type-2 sine basis (DST-II) on the padded box: transfer \(\hat w = \hat q / (\rho_a g + D|k|^4)\) with \(|k|^2 = k_x^2 + k_y^2\), \(k_x = \pi j/(N_x\,d_x)\) and \(k_y = \pi i/(N_y\,d_y)\) for modes \(i, j = 1 \ldots N\), evaluated with scipy.fft dstn/idstn (type=2, workers=-1). Type 2 rather than type 1 because scipy computes DST-I through a real FFT of length 2(N+1), near-prime for a box sized as a fast FFT length, while DST-II runs at N itself: same fidelity, 2-3.5x faster (measurements in the 2026-09-16 decision record). Hardcoded E = 1e11 Pa, nu = 0.25, g = 9.81 – the 9.81 MATCHES fastscape’s flexure (flexure2D.f90), deliberately NOT siim’s constants.GRAVITY = 9.8 (which the shared step’s ice-column term uses); this cross-constant seam is intentional and documented.

The sine basis clamps w = 0 at the padded-box edge (half a sample outside it in type 2, one sample in the fortran’s type 1; invisible at the gates’ resolution) – the flexure2D.f90 basis, here on the native grid – and needs no regime switch (decision, Eric 2026-09-11, superseding the 2026-07-13 domain-mean zeroing):

• domains L << alpha (siim’s valleys; flexural parameter alpha = (4D/(rho_a g))**(1/4), ~55 km at Te = 20 km vs ~20 km domains): a uniform load is RIGIDITY-SUPPORTED by the plate held flat at the padded-box edge, not Airy-compensated – 0.004 m per metre of uniform unloading on a 20 km all-free box, matching the fortran’s far-field-neutral behaviour (full Airy would be 0.875).

• domains >> alpha: a localized load compensates LOCALLY at Airy and its deflection decays to zero well inside the box – the correct large-domain limit. The retired periodic (rfft2) solve zeroed the domain-mean [0, 0] bin, which removed the Airy mean of a localized load and re-emitted it as a uniform uplift of the WHOLE box (+3.5 m per 100 kyr wave step on a 2500 x 250 km domain, an untouched far plateau rising 290 -> 1460 m).

Key differences from the fortran (both by design):

• native grid. The fortran resamples the domain onto the central quarter of a power-of-two grid (~4x coarser), a Nunn & Aires (1988) anti-wraparound device; the native-grid solve is finer and closer to the analytic Kelvin solution. So the two agree only to a tolerance, not bit-for-bit (twin-gated on square grids; the anisotropic gate is vs the Kelvin oracle).

• anisotropy fixed (OQ-6). The fortran builds the y-wavenumber with the x-spacing (pihy = pi/hx, flexure2D.f90:88,145), correct only for square cells. Here k_y uses the true y-spacing dy, so on dx != dy grids the in-house diverges from fortran BY DESIGN and instead tracks the closed-form point-load solution.

Same call signature as the fortran seam flexure(elev_post, elev_eq, nx, ny, xl, yl, rhos, rhoa, Te, ibc) – mutates elev_post in place (adds the deflection w) – so it drops into the S1 injection seam (siim._core.step.glacial_flexure_step()) unchanged.

numpy/scipy only – no fastscape/xsimlab imports (framework-free core).

siim._core.flexure._pad_load(load, free_top, free_bottom, free_left, free_right)

Embed load (ny, nx) in a >=2x, fast-FFT-length grid: mirror-reflect the load across each FREE (no-reflection at fixed) edge – a symmetric extension that gives the zero-gradient (free) plate edge the fortran addw reflection imposes – and zero-pad the fixed edges (clamped far field) plus the anti-wraparound buffer. Returns (padded, off_y, off_x) where the native block sits at padded[off_y:off_y+ny, off_x:off_x+nx].

WARNING – since the sine basis the pad RATIO is physics, not just a numerical buffer: the padded box’s outer edge is where the sine basis clamps w = 0, so it sets how far from the load the plate is held flat, which IS the small-domain rigidity support. Measured on the 31x31 / 20 km all-free config (uniform 1 m unload), mean|w| / Airy runs 0.005 at the current >=2x, 0.028 at 3x, 0.081 at 4x, 0.35 at 6x. Do not retune the RATIO for speed without re-pinning test_flexure_mean_load_rigidity_supported().

siim._core.flexure._transfer(Ny, Nx, dy, dx, D, rhoa)

rho_a g + D |k|^4 on the padded (Ny, Nx) box: DST-II mode j has wavenumber pi*j/(N*d), j = 1..N (module docstring). Cached because it is about a quarter of a solve and identical every step of a run; read-only so a cached array can never be corrupted in place.

siim._core.flexure.flexure(elev_post, elev_eq, nx, ny, xl, yl, rhos, rhoa, Te, ibc)

Flexural deflection of the elastic plate under the load implied by elev_post - elev_eq (the per-step loading increment the caller stacked into elev_post); the deflection w is ADDED to elev_post in place (so rebound = elev_post - elev_eq_pre on return). Signature matches the fortran fs.flexure seam.

elev_post / elev_eq / rhos are flat length-nx*ny arrays (row- major, x fastest, as fastscape passes them); rhoa asthenospheric density; Te effective elastic thickness; ibc the fastscapelib boundary code. dx = xl/(nx-1), dy = yl/(ny-1).

Flow routing & accumulation primitives — siim._core.routing

Flow routing & accumulation primitives for the siim numerical core.

The fastscape-free numba kernels – SFR + D-inf accumulators, the eps-fill priority flood, and the Tarboton (1997) D-inf routing primitives. Pulling them out of the fastscape-importing model is what lets the solver kernels import without the fastscape stack. The DinfFlowRouter xsimlab process that consumes them lives in siim.fastscape.processes. numpy/numba only.

siim._core.routing._flow_accumulate_dinf(field, stack, nb_receivers, receivers, weights)

Single-field D-inf accumulator. Donor-first stack.

siim._core.routing._heap_push(hz, hi, hn, z, i)

Binary min-heap push on parallel (z, idx) arrays; returns new size.

siim._core.routing._heap_pop(hz, hi, hn)

Binary min-heap pop; returns (z, idx, new_size).

siim._core.routing._priority_flood_eps(z_flat, ny, nx, interior_flat, eps, wrap_y, wrap_x, z_fill)

Priority-flood depression filling with an epsilon drainage gradient (Barnes, Lehman & Mulla 2014): flood inward from the outlet cells (interior_flat == 0, the SFR self-receiving set), each cell filled to max(z, spill_path + eps), so on the FILLED surface every interior cell has a strictly lower 8-neighbour and lakes drain toward their spill on eps-gradients. eps accumulates to ~eps * lake diameter — millimetres for any realistic basin. 8-connectivity matches the D-inf facet neighbourhood, so routing the filled surface leaves no interior pits. Looped axes wrap (wrap_y rows, wrap_x columns).

siim._core.routing._dinf_route(z_flat, ny, nx, dx, dy, interior_flat, rec1, rec2, w1, w2, len1, len2, slope_out, e1_dj, e1_di, e2_dj, e2_di, wrap_y, wrap_x)

D-infinity routing. Fills rec1, rec2, w1, w2, len1, len2, slope_out in place.

The Tarboton [Tar97] algorithm on a (possibly anisotropic) rectangular grid: per facet, the cardinal step d1 is dx or dy depending on whether e1 is an E/W or N/S neighbour, and the transverse step d2 is the other spacing:

s1 = (z_c  - z_e1) / d1
s2 = (z_e1 - z_e2) / d2
r  = atan2(s2, s1)                   # flow angle within facet
if r < 0:            r = 0,    s = s1            (along the cardinal)
if r > atan2(d2,d1): r = that, s = (z_c-z_e2)/sqrt(dx²+dy²)  (diagonal)
else:                s = sqrt(s1² + s2²)

The facet with the largest s is selected; flow splits between e1 and e2 in proportion to the angle within the facet (w_e2 = r / facet_angle — the π/4 of the square-grid formula generalises to atan2(d2, d1)). Boundary cells self-receive. Looped axes wrap. Routed on the eps-filled surface (see _priority_flood_eps), no interior cell pits.

Row-parallel (prange over j): every cell writes only its own index and reads only z_flat, so the result is independent of thread count and bit-for-bit identical to the serial scan.

siim._core.routing._dinf_topo_stack(rec1, rec2, w1, w2, n, stack)

Topological sort (Kahn’s algorithm) for D-inf graph. Emits cells in receivers-first order (base cells first, donors last). Caller must reverse for fastscape’s donor-first convention.

siim._core.routing._dinf_pack(rec1, rec2, w1, w2, len1, len2, n, receivers, weights, lengths, nb_receivers)

Pack D-inf routing into (n, 2) receiver/weight arrays.

Single-receiver cases (pit/boundary, both-receivers-same-cell, only one weight nonzero) get nb_receivers[i]=1 with the receiver at index 0. Two-receiver cells get nb_receivers[i]=2.

siim._core.routing._d8_receivers(z_flat, ny, nx, dx, dy, interior_flat, rec, lengths, wrap_y, wrap_x)

Steepest-descent D8 receiver, replicating fortran find_receiver (FlowRouting.f90:302-362) on the (eps-filled) surface: outer jj in {-1,0,1} x inner ii in {-1,0,1} scan (skipping (0,0)), l = sqrt((dx*ii)**2 + (dy*jj)**2), keep on STRICT slope > smax with smax initialised to tiny — so steepest descent wins and, on a tie, the first neighbour in scan order does (byte-identical to fortran on tie-free surfaces). Non-interior cells self-receive (rec[i]=i, lengths[i]=0, the bounds_bc border contract). On the eps-filled surface every interior cell has a strictly-lower 8-neighbour, so it always finds one (no interior self-receivers => acyclic graph). Looped axes wrap.

Serial: each cell writes only its own rec/lengths and reads only z_flat, so the result is independent of iteration order (a prange over the outer loop would be bit-for-bit identical; kept serial per Map 4 §4).

siim._core.routing._d8_stack(rec, n, stack)

Outlet-first topological stack for the single-receiver graph: base cells (rec[i]==i) first, then their donors (receiver-before-donor). Donor-CSR BFS from the base cells — a valid topological order, which is all the accumulators/eroders require (fortran’s exact DFS visitation order is NOT contracted, Map 3 §3(iii)). On the eps-filled surface rec is a forest of trees rooted at the border/base cells, so the BFS emits every cell exactly once.

siim._core.routing._d8_basin(rec, stack, n, basin)

Basin id of each cell = the index of the base (outlet) cell its receiver chain terminates at. One forward pass over the outlet-first stack (the receiver is always processed before its donors), so basin propagates from each outlet down its tributaries. Labeling by outlet INDEX makes the ids reproducible run-to-run (unlike fortran’s unseeded random_number catch labels, Map 4 §4) — but basin is a diagnostic and stays out of every equality gate regardless.

siim._core.routing.d8_interior_mask(border_status, ny, nx)

The interior mask (1 = routable interior cell, 0 = self-receiving boundary) directly from border_status — the in-house replacement for the fortran-derived sfr_rec != arange mask (Map 3 §4). Provably identical: after fortran LocalMinima the only self-receiving cells are bounds_bc = exactly the ‘fixed_value’ border rings (FastScape_ctx.f90:589-592).

‘core’ and ‘looped’ are both plain interior for EVERY router (OQ-1(b), ratified): ‘core’ is non-periodic interior (periodicity is the separate wrap axis, keyed on ‘looped’ only), aligning the SFR mask with the D-inf convention siim already shipped. border_status = [left, right, top, bottom]; node index = j*nx + i.

siim._core.routing._levels_sfr(stack, rec)

Level index for the SFR graph. Returns (order, offsets, nlev): order[offsets[l]:offsets[l+1]] are the nodes of level l.

siim._core.routing._levels_dinf(stack, nb_receivers, receivers)

Level index for the D-inf graph (level = 1 + max over real receivers). Same return contract as _levels_sfr(); nodes are bucketed in receivers-first (reversed donor-first stack) order.

Composition-chain step functions — siim._core.step

Framework-free extractions of the 2D model’s per-step composition chain (one self-free function per siim.fastscape process run_step body). The xsimlab adapter shells and the standalone driver call the same functions.

Framework-free composition-chain step functions for the siim 2D model.

Each run_step body of siim’s own @xs.process classes (siim.fastscape.processes) is extracted here as a self-free module-level function taking state + params explicitly. The xs.process classes become thin shells that unpack self.*, call the function, and assign the outputs; the standalone in-house driver (added later in the migration) calls the same functions. One implementation, two front ends — no divergence.

This module stays numpy/numba/scipy-only (like the rest of siim._core): it imports NOTHING from xsimlab / fastscape / siim.fastscape. The single fortran seam that survives S1 — the FFT flexure solve — is injected as a callable (glacial_flexure_step()’s flexure_solve argument), so this module never imports fastscapelib_fortran either.

The step-order and state-separation subtleties that must remain bit-for-bit stable are summarized below:

• the ice_thickness one-step lag (router/accumulator see H(t-1) — an ordering artifact reproduced by reading H into the routing surface before the kernel overwrites it);

• exactly ONE routing_relax EMA update per step (ema_thickness(); the EMA carry _H_eff is cross-step state AND reused within-step by the trunk subclass — the shell owns that single update);

• size-1 clock-sliced inputs de-squeezed inside the extracted fns (bl = float(asarray(bl).ravel()[-1]), bbu[0] on a (1, ny, nx) slice);

• the FIREWALL (measured twice): the mode-B kernel reconstructs its OWN zs = zb + hc*H from raw H; the relaxed / fabricated surfaces (ema_thickness() / _fabricate_trunk_surface()) feed ONLY the router graph + the mass-balance surface, never a flux closure, the carve, the flexure load, or the outputs.

siim._core.step.build_glacial_params(*, sliding_law, Ko, ce, n, nu, m, mu, Ac, alpha_g, lambda_p, lambda_c, tau_c, coulomb_clamp, hc_over_H, H_diffusivity)

(law_code, GlacialParams) for the given sliding law — the frozen per-run scalars the law_code step skeletons consume. Body of GlacialLaw.initialize + _glacial_params_and_code. Validates sliding_law and hc_over_H; m defaults to n/2 and mu to the per-law constants.derive_*.

siim._core.step.initial_topography(elevation_init, shape, border_status, seed, noise_amplitude)

Initial surface = elevation_init + uniform tie-breaking noise, zeroed on ‘fixed_value’ edges. RNG seeded here (init-only; no per-step RNG).

siim._core.step.plateau_profile(x, plateau_zo, plateau_dz, plateau_frac, plateau_w)

1-D arctan-smoothed plateau profile on the coordinate x: exactly 0 at x[0], rising across an escarpment centred at x_esc = (1 - plateau_frac)*Lx with transition width plateau_w, then dropping plateau_dz linearly to exactly plateau_zo - plateau_dz at x[-1]. The raw arctan ramp is rescaled so the domain ends land on 0 / 1 (the raw ramp only reaches them asymptotically — with the default frac/w on a 50 km domain the low edge sat at 0.25*zo, a permanent sill above the border’s water datum).

siim._core.step.plateau_surface(x, y, shape, border_status, seed, noise_amplitude, plateau_zo, plateau_dz, plateau_frac, plateau_w)

Arctan-smoothed plateau initial topography (plateau_profile() tiled along y) + tie-breaking noise (zeroed on ‘fixed_value’ edges, so a fixed x-border starts exactly on 0 / plateau_zo - plateau_dz). Body of PlateauSurface.initialize.

siim._core.step.wave_uplift(x, y, shape, mask, dt, t, delta_h, wave_width, wave_velocity, x_escarpment, wave_calibration, U_inf)

Moving-Gaussian uplift wave, midpoint-sampled (t + dt/2) with calibration 1.0 = exact delta_h deposition over the passage. Body of WaveUplift.run_step; the border mask is uplift_mask().

siim._core.step.SIDE_SLICES = ((slice(None, None, None), 0), (slice(None, None, None), -1), (0, slice(None, None, None)), (-1, slice(None, None, None)))

Border ring of each domain edge, indexed in border_status order (siim._core.outputs.SIDES = left, right, bottom, top): index 2 is row 0 (y = 0), index 3 is row ny-1.

siim._core.step.uplift_mask(border_status, shape)

Binary uplift mask: 0 on ‘fixed_value’ border rings, 1 elsewhere. Body of fastscape BlockUplift.initialize (the mask half).

siim._core.step.block_uplift(rate, dt, mask, shape)

uplift = rate * dt, zeroed on fixed borders by mask. Body of GlacialBlockUplift.run_step; drops the leading size-1 tstep dim left by xsimlab on a (nt, y, x) slice before broadcasting.

siim._core.step.ema_thickness(H_lag, H_eff_prev, r)

The (optionally EMA-relaxed) lagged thickness feeding this step’s routing + mass-balance surface. r == 0 returns the raw lagged H unchanged (bit-for-bit); otherwise H_eff = r*H_eff_prev + (1-r)*H_lag (seeded at H_lag on the first step, H_eff_prev is None).

The caller owns the single per-step update: it stores the return as the EMA carry (cross-step) AND reuses the SAME value within-step (the trunk subclass). One update per model step — never call twice.

siim._core.step.routing_surface(post_uplift_surface, hc_over_H, H_eff)

The mode-B routing/mass-balance surface zs = post_uplift_bed + hc*H_eff (the reconstructed ice column on the post-uplift bed). H_eff is the relaxed lagged thickness (raw H when routing_relax == 0). FIREWALL: this surface reaches ONLY the router graph + the accumulator’s mass-balance surface — never a physics closure, the carve, the flexure load, or outputs.

siim._core.step._fabricate_trunk_surface(zs_dyn, zb, H_lag, border, alpha_g, dx, dy, k_dip, floor, offsets, D, SRC, wrap_y, wrap_x)

Build the fabricated trunk routing surface described in docs/guides/concepts.md. Pure function (no process state) so the process and the channel-persistence test share one implementation.

zs_dyn (ny, nx) is the dynamic ice surface zb + hc*H_lag; zb the matching bed; H_lag the lagged thickness; border (ny, nx bool) the base-level edges (excluded as seeds + never fabricated). offsets/D/SRC are (ny, nx) scratch buffers. Returns a fresh (ny, nx) elevation: the V-dipped trunk surface at footprint cells (max(zs_geo, zb)), zs_dyn elsewhere.

siim._core.step.accumulate_glacial_flow(surface, surface_upward, zELA, beta, runoff, cell_area, width_hack_k, width_hack_p, shape, stack, receivers, nb_receivers, weights, lengths, basin)

Accumulate water flux, ice flux, drainage area + the routing-topology outputs (basin_ids, receivers_2d, stack_2d). Body of GlacialFlowAccumulator.run_step. surface is the post-uplift routing surface (zs_route); the ELA-relative mass balance b(z) is evaluated on the PRE-uplift climate surface z_clim = surface - surface_upward (no O(U*dt) bias). Returns (ice_flux, water_flux, area, basin_ids, receivers_2d, stack_2d) (caller also sets flowacc = water_flux).

siim._core.step._H_from_QS_modeA(law_code, gp, Qg, S)

Per-law point closure H(Qg, S): the from_slope branch of the shared _modeb_closure() (single-sourced; audit m34).

siim._core.step._solve_border_H_modeA(z_flat, H_flat, ice_flux, receivers, nb_receivers, lengths, law_code, gp)

Mode A: a self-receiving border node with through-flowing ice gets its thickness from the per-law H(Q, S) closure with S the steepest upwind (donor-side) surface slope — matching siim1d’s outlet treatment. Mutates H_flat in place at border cells; the interior is untouched. nb_receivers is read only on the D-inf branch (may be None under SFR).

siim._core.step.run_modeA_step(surface, H, ice_flux, water_flux, law_code, gp, dt, stack, receivers, nb_receivers, weights, lengths, hc_over_H, shape)

Mode A erosion step (ice-surface state). Body of GlacialSPLModeA: dispatch to the routing-specific mode-A skeleton (mutating a fresh copy of the surface + H), solve the border-H closure, then commit. Returns (z_eroded, H_new, bedrock_surface, erosion, denudation) — erosion == denudation (mode A erodes the ice surface as its single hc-invariant state).

siim._core.step.run_modeB_kernel(zb_flat, H_flat, ice_flux, water_flux, law_code, gp, dt, stack, receivers, nb_receivers, weights, lengths, shape, dx, dy, border_bed_uplift, bl, gate, ramp, parallel, wrap_y, wrap_x)

No-carve mode-B kernel dispatch (SFR or D-inf). Mutates zb_flat and H_flat in place; returns the kernel’s surface_out (= zb + hc*H).

FIREWALL: the kernel reconstructs its OWN zs = zb + hc*H from the raw H for every closure and erosion slope — it never sees the relaxed/fabricated routing surface. De-squeezes the size-1 clock slices of border_bed_uplift and bl inside, so the shell and the driver share the treatment. bl may also be the driver’s flat (ny*nx,) per-side border datum (see siim._core.driver._bl_field()).

siim._core.step.carve_bed(zb_flat, H_flat, surface_out, zb_pre, receivers, alpha_g, hc_over_H, widening_factor, shape, dx, dy, wrap_y, wrap_x, offsets=None, D=None, SRC=None, zb_kern=None)

Apply the sub-grid width carve (see siim._core.carve) to the post-kernel bed zb_flat IN PLACE. zb_pre is the pre-kernel bed (the denudation datum + descent-cap origin); surface_out the kernel’s reconstructed ice surface (updated for carved cells). Routing-agnostic: receivers enter only through the border marker rec[i] == i (self-receiving = base-level border). offsets/D/SRC/zb_kern are optional reusable scratch buffers (an optimization, not semantics — allocated fresh if omitted).

siim._core.step._filled_surface_and_mask(elevation, shape, border_status)

The fill-then-route preamble shared by route_d8() / route_dinf(): the border_status interior mask (Map 3 §4 — the in-house replacement for the fortran sfr_rec != i mask), the looped-axis wrap flags, and the eps-filled surface (depression floors filled to spill + eps, so every interior cell has a strictly-lower 8-neighbour and lakes drain toward their spills). Returns (z_route, interior, wrap_y, wrap_x, ny, nx, nn). The kernels/physics keep consuming the TRUE elevation; only routing sees the fill.

siim._core.step.route_d8(elevation, shape, dx, dy, border_status)

In-house D8 single-flow router (fill-then-route) — the framework-free replacement for fastscape’s fortran SingleFlowRouter (fs.flowroutingsingleflowdirection). Returns the SFR router bundle (receivers, weights, lengths, nb_receivers, stack, basin): receivers / lengths 1D (n,), weights / nb_receivers all-ones (SFR is single-receiver), stack outlet-first, basin (ny, nx) labeled by outlet index. The receiver scan replicates fortran find_receiver on the eps-filled surface (Map 3 §3); the routing delta vs fortran is confined to depression/tie cells (behavioral gate).

siim._core.step.route_dinf(elevation, shape, dx, dy, border_status)

D-infinity flow directions on the eps-filled surface. Body of DinfFlowRouter.run_step — fully fortran-free (S4): the interior/boundary mask comes from border_status directly (Map 3 §4, provably identical to the old fortran sfr_rec != i mask), and basin from the in-house outlet labeling. Returns (receivers, weights, lengths, nb_receivers, stack, basin) — receivers/weights/lengths (n, 2), basin (ny, nx).

siim._core.step.accumulate_sediment(denudation, cell_area, stack, receivers, nb_receivers, weights, shape)

Route this step’s denuded rock volume max(denudation, 0) * cell_area down the flow graph in one accumulation pass. Body of SedimentTracker.run_step (the per-step flux; the caller owns the cross-step running integral _cum).

siim._core.step.edge_sediment(flux, border_status)

Sum the routed per-node flux (the (ny, nx) field accumulate_sediment() returns) over each domain-edge outlet ring — the volume (m^3) leaving the domain across that edge this step.

Returns a length-4 array in border_status order (siim.siim2d.siim._BL_SIDES = left, right, bottom, top), NaN on any side that is not 'fixed_value': only fixed-value rings are self-receiving outlets (siim._core.routing.d8_interior_mask()), so a ‘core’ or ‘looped’ edge has no sediment to deliver, which is not the same as delivering none. Each corner node is counted exactly ONCE, by the same rule siim._core.driver._bl_field uses for the water datum: an outlet x-side (left/right) owns it, else the y-side.

Shared by the in-house driver and SedimentTracker.run_step so the two front ends report the identical sums.

siim._core.step.glacial_flexure_step(elevation, denudation, surface_upward, ice_thickness, alpha_g, cell_area, lithos_density, asthen_density, e_thickness, ibc, shape, length, col_prev, ice_load, flexure_solve)

Incremental flexural isostasy: rock unloading (surface_upward - denudation) plus, when ice_load, the per-step glacial ice load (rho_ice/lithos)*d(col) with col = alpha_g*H**2/L (hc-free, mass-conserving). Body of GlacialFlexure.run_step; the biharmonic plate solve is INJECTED as flexure_solve(elev_post, elev_eq, nx, ny, xl, yl, lithos, asthen, Te, ibc) (fortran fs.flexure at S1, in-house FFT later — keeps this module framework-free). Returns (rebound, col_new); the caller keeps col_new as the cross-step _col_prev.

The returned rebound is exactly 0 on every 'fixed_value' border row and column (read off ibc, so the in-house driver and the xsimlab adapter share it) — those nodes already get no block uplift and no erosion, so letting them subside was the one place their value was not actually fixed.

siim._core.step.sum_erosion(*erosion_terms)

Reproduce fastscape TotalErosion.height = sum(erosion group) = glacial_spl.erosion + diffusion.erosion. Sums from 0 like the stock builtin sum, so a single term returns a copy (0 + a == a).

siim._core.step.compose_vertical_motion(surface_forcing, bedrock_forcing, rebound, erosion_total)

Reproduce fastscape’s TotalVerticalMotion group-sum composition:

surface_up = sum(surface_upward group) - sum(surface_downward group)
           = (surface_forcing + rebound) - erosion_total
bedrock_up = sum(bedrock_upward group)  = bedrock_forcing + rebound

surface_forcing / bedrock_forcing are the tectonic forcings (TectonicForcing.{surface,bedrock}_upward = block uplift here); rebound the flexural rebound in both upward groups (None when flexure is off); erosion_total the TotalErosion.height (sum_erosion()). The finalize commit topo += surface_up lives in the driver.

In-house time loop — siim._core.driver

The standalone driver: siim’s own merged step loop + two-cadence snapshot, calling the same step functions the adapter shells call. Selected via run(driver='inhouse') (default constants.DRIVER_DEFAULT).

The in-house time loop for the standalone siim 2D model.

Owns the merged step chain and the two-cadence snapshot that xsimlab’s driver performed for the siim.fastscape adapter path, calling the SAME framework-free step functions (siim._core.step) the @xs.process shells call — one implementation, two front ends. The public model and output contracts are documented in docs/guides/concepts.md and docs/guides/outputs_and_io.md.

Merged step order (a valid topological sort of the composition graph):

uplift -> tectonics -> surf2erode -> route -> accumulate -> kernel ->
diffusion -> flexure -> (erosion sum + vertical-motion compose) -> sediment
-> [SNAPSHOT] -> finalize (topo += surface_up)

Snapshot-timing invariant (load-bearing, Map 2 §2): intermediate frames are captured AFTER the step but BEFORE the uplift-committing finalize; the LAST frame (out_idx[-1] == nt-1) is captured AFTER the loop (post the final finalize). Frame 0’s topography__elevation is thus the initial topo while H / flux / area / erosion already reflect one solved step.

Cross-step state OWNED by the driver: the ice_thickness one-step lag (the router/accumulator read H(t-1) — an ordering artifact reproduced by reading H into the routing surface before the kernel overwrites it), the _H_eff EMA carry (routing_relax), the flexure _col_prev, the sediment _cum (and its per-domain-edge twin), and the topography state (committed only at finalize).

FIREWALL (Map 1 §2): the relaxed / fabricated routing surface (zs_route) reaches ONLY the router graph, the mass-balance surface, and the ice-surface hillslope diffusion — the mode-B kernel reconstructs its OWN zs = zb + hc*H from raw H for every closure/erosion slope.

Routing is INJECTED as a callable (cfg.route) — the fortran backend at S3, the in-house D8 producer at S4 — so this module stays framework-free (no xsimlab / fastscape / fastscapelib import). The flexure plate solve is injected too (cfg.flexure_solve); the in-house hillslope diffuser (siim._core.hillslope.diffuse()) is called directly.

siim._core.driver._edge_border_mask(border_status, shape)

Non-looped domain edges (base-level borders), as the trunk-surface fabrication excludes them (TrunkSurfaceToErode.initialize).

siim._core.driver._slice_forcing(series, scalar, k)

Per-step forcing value: series[k] when a clock series is present, the static scalar/field otherwise (the driver indexes arr[k] directly, replacing xsimlab’s (('tstep',), arr) groupby slice, Map 2 §4).

siim._core.driver._bl_field(bl_sides, shape, k)

Per-node water datum for step k (PER-SIDE bl): each fixed_value border node carries its own side’s datum; interior entries are placeholders (0) that the mode-B kernel overwrites with the datum of the outlet its basin drains to. Corner nodes touched by two fixed sides take the x-side (left/right) value — the x-sides are written last. bl_sides is the 4-list of (scalar, series) entries parsed by siim.siim2d.siim.set_and_check_parameters (None off fixed sides).

siim._core.driver.run_loop(cfg)

Run the merged in-house time loop; return the packed ds_out step buffers (a {name: ndarray} dict, one array per output-spec row at that row’s dims). cfg is the resolved-parameter bundle assembled by siim.siim2d.siim._run_inhouse() (see that method for the field contract).

siim._core.driver._snapshot(buffers, jf, cur)

Copy this frame’s step values into the preallocated output buffers (dtype cast happens on assignment: int32 buffers, float64 fields).

Output packing — siim._core.outputs

Packs the driver’s step buffers into the exact ds_out contract (variable names, (time, y, x) dims, dtypes, coords) that siim.siim2d.siim._unpack_outputs reads.

Pack the in-house driver’s step buffers into the exact ds_out contract.

The standalone driver (siim._core.driver) fills per-step numpy buffers; this module allocates them to the exact dtypes and wraps them in the xarray.Dataset siim.siim2d.siim._unpack_outputs() reads — same variable names, dims ((time, y, x) for every per-node field, (time, side) for the per-domain-edge sediment totals), coords (time/x/y/tstep, plus side), dtypes (float64 fields; int32 for basin_ids/receivers_2d/stack_2d) and the mode/flag-conditional rows. _unpack_outputs / save / load consume the exact keys constructed here; see docs/guides/outputs_and_io.md for the public output contract.

xarray is the only heavy import; it stays function-local so import siim._core.outputs costs nothing until a dataset is built.

siim._core.outputs.SIDES = ('left', 'right', 'bottom', 'top')

domain edges in border_status order. Index 2 (‘bottom’) is row 0, i.e. y = 0, and index 3 (‘top’) is row ny-1 — the fastscape ibc digit naming calls those same two edges the other way round.

Type:

The side coordinate

siim._core.outputs.output_spec(mode, flexure, sediment, sediment_edge)

The active (name, dtype, dims) rows for this run. bedrock_surface iff mode A (keys the _unpack mode-A branch, siim2d.py:965), sediment__* iff sediment, sediment__edge_* (the per-edge totals, dims (time, side)) iff sediment_edge, flexure__rebound iff flexure.

siim._core.outputs.allocate_buffers(spec, nt_out, ny, nx)

One output buffer per spec row, at the row’s dtype and dims (the raster rows (nt_out, ny, nx), the per-edge rows (nt_out, 4)).

siim._core.outputs.build_dataset(buffers, spec, t_out, x, y, tstep)

Wrap the filled step buffers into the ds_out xarray Dataset: each row at its spec dims; coords time=t_out, x, y, tstep, plus side when a per-edge row is active. Reads only variable values by name downstream (_unpack_outputs), so this in-memory build is interchangeable with the retired zarr path.