fastscape integration — siim.fastscape

The public, composable surface for using siim’s glacial-erosion physics inside your own fastscape/xsimlab model: the assembly helpers (glacial_processes, glacial_model) and the @xs.process classes they wire in. This is the optional adapter (pip install siim-lem[fastscape]; conda for the Fortran backend) — the standalone siim.siim2d model runs its own in-house driver and does not require it, though siim2d(...).run(driver='xsimlab') will drive the model through this adapter when the stack is present.

Public surface — siim.fastscape

siim.fastscape — siim’s glacial processes as composable fastscape citizens.

This is the public surface for using siim’s glacial-erosion physics inside your own fastscape/xsimlab model. The names exported here are @xs.process classes that drop into fastscape’s basic_model once its stock 'spl' and 'drainage' slots are removed. Build the override dict with glacial_processes(), or get a ready-to-run model from glacial_model():

from fastscape.models import basic_model
from siim.fastscape import glacial_processes

model = basic_model.drop_processes(['spl', 'drainage']).update_processes(
    glacial_processes(mode='B'))

For the full coupled model with forcing, I/O and plotting, use siim.siim2d.siim instead. Its default in-house driver calls the shared framework-free numerical core directly; its optional driver='xsimlab' path assembles these process classes. The two paths share the same law parameters and step kernels and are kept aligned by driver-conformance tests.

Implementation lives in submodules: the glacial @xs.process classes in siim.fastscape.processes, and two reusable generic forcing processes (WaveUplift, PlateauSurface) in siim.fastscape.forcing.

API stability

siim is pre-1.0 (0.x). This surface is intended to stay stable, but may still evolve before 1.0 — it is not a frozen contract yet. Modes A, B and C are settled: the sub-grid width-carving path (carve=True → GlacialSPLModeC) is a mode-B citizen too (bed + ice thickness + carve), and is siim’s flagship carved mode.

siim.fastscape.glacial_processes(*, mode='B', carve=None, routing='single', router_backend=None, flexure=False, sediment=False, trunk_surface=None, numerics_backend='inhouse')

Build the slot→process-class override dict for fastscape’s basic_model.

Apply the returned dict after dropping the stock fluvial slots:

basic_model.drop_processes(['spl', 'drainage']).update_processes(
    glacial_processes(...))

or use glacial_model(), which does exactly that. The optional driver='xsimlab' path of siim.siim2d.siim builds its model from this function; the default in-house driver consumes the same shared law records and numerical kernels without xsimlab.

Every flag defaults off, except that mode C (mode='C', or mode='B' with carve=True) resolves trunk_surface to the mode-C standard exactly as siim.siim2d.siim does (siim.constants.mode_c_standard()), so glacial_model(mode='C') with no extra inputs is native mode C.

Parameters:

• mode ({'B', 'A', 'C'}) – One-breath taxonomy: 'A' = ice-surface state (carved troughs heal instantly); 'B' = bedrock zb + ice-thickness H state (bed memory; siim’s native regime); 'C' = B + sub-grid width carving (equivalent to mode='B', carve=True — implies carve, and an explicit carve=False contradicts it). Case-insensitive.

• carve (bool, optional) – Mode-B sub-grid glacier-width carving — wires GlacialSPLModeC (the citizen bed-and-thickness class plus the carve). Ignored for mode='A'. Defaults off here so the advertised default is the plain Mode-B citizen. siim.siim2d.siim defaults to Mode C, which enables carving; explicit Mode B leaves it off unless requested. mode='C' is the alias for turning it on.

• trunk_surface (bool or None, optional) – Fabricated trunk-surface routing (mode B/C) — swaps in TrunkSurfaceToErode so routing + accumulation converge trunk flow onto the centerline (the centerline’s raw flux is then the full cross-section). None (default) resolves to the mode-C standard: on for mode C, off for plain mode B. An explicit True / False wins. routing_relax follows the provider class here (0.5 on TrunkSurfaceToErode, 0 on GlacialSurfaceToErode) where siim.siim2d.siim follows the mode, so two explicit combinations differ unless surf2erode__routing_relax is set: mode B with trunk_surface=True (0.5 here, 0 in siim2d) and mode C with trunk_surface=False (0 here, 0.5 in siim2d). Ignored for mode='A'.

• routing ({'single', 'dinf'}, optional) – D8 single-flow ('single', default) or D-infinity ('dinf').

• router_backend ({'inhouse_d8'}, optional) – Single-flow routing backend (constants.ROUTER_DEFAULT). 'inhouse_d8' — siim’s numba D8FlowRouter — is the only accepted value since the 0.9.1 standalone flip (the fortran SingleFlowRouter wiring was retired); the parameter is the router-contract plug point for future backends. Only affects routing='single' — the D-inf directions + mask + basin ride the same in-house contract.

• flexure (bool, optional) – Opt in to glacial-isostatic flexure / sediment tracking.

• sediment (bool, optional) – Opt in to glacial-isostatic flexure / sediment tracking.

• numerics_backend ({'inhouse'}, optional) – Flexure + hillslope-diffusion backend (constants.NUMERICS_BACKEND). 'inhouse' — HillslopeDiffusion (siim’s numba ADI) in the stock LinearDiffusion slot and GlacialFlexure on siim’s scipy.fft plate solve — is the only accepted value since the 0.9.1 standalone flip. Does not affect routing.

Returns:

Maps xsimlab slot name → process class, ready for update_processes.

Return type:

dict

siim.fastscape.glacial_model(*, mode='B', carve=None, routing='single', router_backend=None, flexure=False, sediment=False, trunk_surface=None, numerics_backend='inhouse')

Return a ready-to-run xsimlab Model: fastscape’s basic_model with its 'spl'/'drainage' slots replaced by siim’s glacial processes.

Convenience wrapper over glacial_processes() (same parameters). For the full coupled model with forcing, I/O and plotting, use siim.siim2d.siim instead.

Glacial processes — siim.fastscape.processes

xsimlab process classes for siim2d.

Three processes plug into fastscape’s basic_model:

• InitialTopography replaces FlatSurface for the initial elevation field (precomputed elevation + zeroed-on-fixed-edges noise).

• GlacialFlowAccumulator extends FlowAccumulator to track water flux, ice flux, drainage area, basin IDs, and the per-step routing topology (receivers_2d, stack_2d) needed by downstream extract_channel / extract_basin / strahler_order code in siim2d.py.

• GlacialSPLModeA / GlacialSPLModeB / GlacialSPLModeC (all GlacialSPLBase subclasses) plug into the model’s erosion slot, one per surface-evolution mode (build-time selection): A tracks the ice surface, B the bed + ice thickness (citizen), C = B plus the sub-grid width carve. They dispatch on (sliding_law, routing) and call the corresponding law_code step skeleton in siim._core.skeleton.

• DinfFlowRouter overrides fastscape’s FlowRouter with Tarboton (1997) D-infinity flow directions (defined at the bottom of this module).

The numba flow-accumulation + D-infinity routing primitives live in siim._core.routing; the ice-thickness + erosion solvers in siim._core (solvers / eroders / skeleton).

class siim.fastscape.processes.InitialTopography

Bases: object

Initialize surface topography from a pre-computed elevation array plus uniform tie-breaking noise. Noise is zeroed on ‘fixed_value’ edges so the boundary z stays exactly at the elevation_init value across the run (BlockUplift and the glacial erosion process both skip fixed cells; without this any IC noise there gets locked in and drifts the base level).

class siim.fastscape.processes.GlacialBlockUplift(*args: Any, **kwargs: Any)

Bases: BlockUplift

BlockUplift that tolerates a time-varying (nt, ny, nx) uplift rate.

xsimlab slices the clock axis of a (('tstep', 'y', 'x'), array) input per step, but xarray no longer squeezes groupby, so the per-step rate keeps a leading size-1 dim and stock BlockUplift’s np.broadcast_to((1, ny, nx), (ny, nx)) raises. Drop that leading dim before broadcasting; scalar / (y, x) rates are unaffected.

class siim.fastscape.processes.GlacialLaw

Bases: object

Raw sliding-law inputs + the derived-constant algebra, hoisted out of the glacial erosion process (GlacialSPLBase).

Owns the per-run physical/law scalars (sliding law, erosion coefficients, ice-rheology parameters, the channel-floor ratio hc_over_H, the channel aspect ratio alpha_g) and turns them into the frozen (law_code, GlacialParams) record the law_code step skeletons consume (siim._core.skeleton). Exposed as params so the erosion process (and any sibling) reads one record by foreign instead of re-deriving the constants. hc_over_H is the centerline/mean channel-depth ratio (default constants.HC_OVER_H); the surfaces built from (zb, H) state are zs = zb + hc_over_H * H (the tracked bed is the sub-grid channel floor, H the width-mean depth).

class siim.fastscape.processes.GlacialFlowAccumulator(*args: Any, **kwargs: Any)

Bases: FlowAccumulator

Accumulates water flux, ice flux, and drainage area through the flow graph.

Sub-grid mass balance: ablation below the ELA scales with the glacier’s plan-view area (width × flow length) wherever that exceeds cell_area. Width is set by Hack-style drainage-area scaling (Hack [Hac57]):

glacier_width = width_hack_k * upstream_area ** width_hack_p

This decouples width from instantaneous ice thickness, so the kinematic- wave mode at the toe doesn’t feed back through ablation. Matches the strictly-elevation-driven mass balance of the 1D analytical SS.

class siim.fastscape.processes.GlacialSPLBase

Bases: object

Glacial + fluvial erosion — shared base for the per-mode erosion processes (GlacialSPLModeA, GlacialSPLModeB, GlacialSPLModeC). Abstract: it owns the shared declarations + setup but defines NO run_step (the concrete subclass that fills the glacial_spl slot supplies it; siim2d picks the subclass at build time by mode). It is @xs.process-decorated only so its variables are discoverable as foreign targets (siblings + subclasses); it is never assigned to a model slot itself.

The tracked bed is the sub-grid channel floor zb = z - hc_over_H*H (H the width-mean depth). The raw law params + constant algebra live in GlacialLaw; this base reads the derived (law_code, GlacialParams) record via the law foreign. Subclasses dispatch on the record’s law_code and on routing dimensionality (SFR vs D-inf), calling the law_code step skeletons in siim._core.skeleton.

Sibling processes foreign onto THIS base (xsimlab resolves foreign(Base) via the MRO to whichever subclass fills the slot): SedimentTracker and GlacialFlexure read denudation / ice_thickness here.

class siim.fastscape.processes.GlacialSPLModeA(*args: Any, **kwargs: Any)

Bases: process

Mode A erosion (ice-surface state). Fills the glacial_spl slot when mode == 'A'. The ice surface is the single tracked state (pinned at base level by the BCs); H is solved from the local pre-erosion surface slope each step, the bed is the derived z - hc_over_H*H. Historical siim2d behaviour; never carves (the bed is reconstructed, not stored).

class siim.fastscape.processes.GlacialSPLModeB(*args: Any, **kwargs: Any)

Bases: process

Mode B erosion as a fastscape citizen (Fork B), no carve.

The tracked state topography__elevation IS the bed zb (the sub-grid channel floor); ice_thickness carries H. Each step hands the kernel the POST-uplift bed (persisted state + this step’s TectonicForcing bed motion; fastscape composes the same uplift at finalize) and H from ice_thickness, runs the UNCHANGED mode-B kernel (which reconstructs zs = zb + hc_over_H*H internally), and reports a genuine erosion height:

denudation = zb_in - zb_out (the bed lowering this step) erosion    = denudation (the erosion-group height)

fastscape then finalizes zb_new = zb + uplift - erosion = zb + uplift - delta-zb = the kernel’s eroded post-uplift bed, so the bed evolves correctly. No bedrock_surface output — topography IS the bed (the derived display surface zb + hc_over_H*H is reconstructed by the consumer / plotter). The sub-grid width carve rides on the GlacialSPLModeC subclass; routing/erosion slopes come from the reconstructed ice surface (GlacialSurfaceToErode).

class siim.fastscape.processes.GlacialSPLModeC(*args: Any, **kwargs: Any)

Bases: process

Mode C erosion: the citizen mode-B bed-and-thickness class (GlacialSPLModeB) plus the sub-grid glacier-width carve. Fills the glacial_spl slot for mode == 'B' with carve_width on (the default siim2d path).

Everything the citizen does is inherited: topography IS the tracked bed zb (bed memory; no flicker, no bedrock_surface output), routing + erosion slopes come off the reconstructed ice surface (GlacialSurfaceToErode), and erosion == denudation == delta-zb so fastscape finalizes zb_new = zb + uplift - delta-zb. This subclass runs the UNCHANGED mode-B kernel and then carves the sub-grid footprint into the bed (the topography state) — so the carve deepening flows straight into denudation (the sediment / flexural-unloading source), the same rock the retired surface-replace carve class removed but now on citizen semantics. The carve edits the BED after the kernel; surface_out (the kernel’s reconstructed ice surface) is handed to the carve for its per-cell updates but discarded by the citizen.

With carve_width=False this is GlacialSPLModeB bit-for-bit (the carve is gated on self._carve); the class exists so the carve rides the citizen path rather than the surface-replace one.

class siim.fastscape.processes.GlacialSurfaceToErode(*args: Any, **kwargs: Any)

Bases: SurfaceAfterTectonics

Surface used for routing + erosion slopes in citizen Mode B (GlacialSPLModeB).

Subclasses fastscape’s SurfaceAfterTectonics, whose elevation = topo_elevation + forced_motion is the post-uplift surface (here topo IS the bed zb, so this is the post-uplift bed). Adds the reconstructed ice column on top:

elevation = (post-uplift bed zb) + hc_over_H * H_lag = zs

i.e. the ice surface the flow router stacks on and the erosion kernel takes its slopes from. The uplift add is inherited (no explicit +uplift); H is the step-start (lagged) ice_thickness global, consistent with the order surf2erode -> glacial_flow -> glacial_spl. hc_over_H is read from the law record via GlacialLaw. (A pre-uplift climate variant topo + hc_over_H*H is deferred to Phase 3.)

Anti-flicker relaxation (routing_relax = r): the once-per-step H -> zs -> D8-receivers -> flux -> closure-H loop lags by one step and D8 receiver choice is discrete, so near-tied surfaces flip whole subtrees of flux each step and H ~ Q^(1/4.. 1/5) turns those O(1) swaps into a period-2 ice-thickness slosh (a cosmetic planview flicker; integrated area/volume + attractor stats are unaffected). With r > 0 the routing thickness is the EMA H_eff(t) = r*H_eff(t-1) + (1-r)*H_lag(t) (seeded H_eff(0) = H_lag) instead of the raw lagged H. r = 0 (default) uses the raw H, bit-for-bit.

State-separation firewall: the relaxed H reaches ONLY this provider’s elevation, which feeds ONLY the flow router graph and GlacialFlowAccumulator’s mass-balance surface. The mode-B kernel reconstructs its own zs = zb + hc*H from the tracked bed + raw H for every closure and erosion slope (it never reads this elevation), so no relaxed/geometric value ever enters a physics closure, the carve, the flexure load or the outputs.

class siim.fastscape.processes.TrunkSurfaceToErode(*args: Any, **kwargs: Any)

Bases: process

Fabricated trunk-surface routing (mode B/C, trunk_surface; on by default for mode C, as part of the mode-C standard).

routing_relax defaults to constants.MODE_C_ROUTING_RELAX here (the parent’s is constants.ROUTING_RELAX), so this class in the surf2erode slot is the whole mode-C routing standard with no extra inputs.

A trunk glacier of mean thickness H occupies width W = alpha_g*H — routinely several cells — but the flow graph carries its ice down 2-3 parallel chains, each with its own diluted ice_flux and its own closure H, understating the trunk H. This provider replaces the plain reconstructed ice surface (zb + hc*H, GlacialSurfaceToErode) with a fabricated surface that has a LINEAR cross-valley dip toward the trunk centerline. The flow router stacks on it and the accumulator reads it, so (a) flow converges onto the centerline chain — its raw accumulated flux then IS the full cross-section discharge (no accounting correction needed), and (b) mass balance is evaluated at the trunk-surface elevation.

Fabrication (per-step, from LAGGED closure-H — this provider runs before the router). The power transform the carve uses gives footprint membership D < 0 and per-cell source SRC (the thickest disc wins — the medial-axis-transform inversion). For a footprint cell i with source s (whose disc radius R_s = alpha_g*H_s/2 > cell):

zs_geo(i) = zs_dyn(s) + S_c(s) * (d_i - R_s)

with d_i the distance to the source (d_i^2 = D_i + R_s^2, from the transform), zs_dyn(s) = zb_s + hc*H_s the source’s own dynamic surface, and S_c(s) = TRUNK_DIP_K * max(|grad zs_dyn|(s), TRUNK_DIP_FLOOR) the cross-slope (0 at the rim d = R_s, -S_c*R_s at the axis). The routing value inside a footprint is max(zs_geo, zb) — bare cells / nunataks present their rock (routing goes around them); a footprint cell’s own dynamic ice column is deliberately IGNORED (its ice is part of the trunk; taking a max with it would re-erect a spine). Outside footprints the surface is the dynamic zb + hc*H unchanged.

State-separation rule: the fabricated surface reaches ONLY the flow router (graph) and the accumulator’s mass-balance surface. The mode-B kernel is UNTOUCHED — it uses the router’s receiver graph but reconstructs its own zs = zb + hc*H for every closure and erosion slope, so no fabricated (geometric) elevation ever enters a closure. A flank cell routed to a taller centerline hits the well-posed negative-a (from-a) branch, bounded by flux; the carve reads post-kernel closure state. The dip is a routing numerics device (like the priority-flood eps), not physics.

class siim.fastscape.processes.SedimentTracker

Bases: object

Optional sediment-throughput tracker (added only when siim’s track_sediment is on — zero cost otherwise).

Routes each step’s denuded rock volume — max(denudation, 0) * cell_area — down the flow graph in one accumulation pass, giving per node the total upstream-eroded volume passing through it that step (flux, m^3) and its running time-integral (cumulative, m^3). Detachment-limited bookkeeping only: no deposition, every node passes all upstream sediment through. Difference cumulative along the time axis to recover per-interval volumes; the outlet node’s cumulative is the whole-basin yield.

edge_flux / edge_cumulative report the same routed flux summed over each domain-edge outlet ring instead (dims ('side',), the border_status order left/right/bottom/top, NaN off 'fixed_value' sides) — one shared sum (siim._core.step.edge_sediment()) with the in-house driver. Both reports ride the one accumulation pass; siim’s track_sediment value picks which are stored, and clears basin when only the edge totals are wanted so the per-node running integral is skipped.

Reads the erosion process’s denudation (via GlacialSPLBase: the true rock removed — delta-zb incl. sub-grid carve in mode B, delta-zs in mode A; per-step, already includes dt) so it runs after the erosion step and reuses the same receiver/stack flow graph.

class siim.fastscape.processes.GlacialFlexure(*args: Any, **kwargs: Any)

Bases: Flexure

Flexural isostasy with true glacial isostatic adjustment: ice loading + erosional/tectonic unloading + elastic rebound.

fastscape’s Flexure loads the elastic plate from TotalErosion.height (the ice-SURFACE change in siim2d mode B — ice-thickness change leaks in). This override does two things:

1. Re-sources the rock load from the erosion process’s per-step denudation (via GlacialSPLBase; mode B: delta-zb incl. sub-grid carve; mode A: delta-zs), so the plate flexes in response to rock erosional + tectonic unloading, correct under transient ice.

2. Adds the per-step ICE load. The glacial ice load is the channel cross-section carried over the cell, col = alpha_g * H**2 / L (L = sqrt(cell_area)) — the mass-conserving, hc-FREE ice volume per cell. alpha_g*H**2 = Qg/V (the ice flux over the channel velocity), and the centerline/mean hc factor cancels the parabola’s 2/3, so it never appears. Its per-step change folds in as a rock-equivalent column, (rho_ice/lithos_density)*d(col), into the diff that drives fastscape’s incremental plate solve (incremental: the full column would re-apply the whole load every step). Width-aware (∝ H**2): a glacier wider than the cell (alpha_g*H > L) piles proportionally more ice than the mean depth H. Gated by ice_load (default True; False = erosion-only, the pre-GIA behaviour). The mass-conserving Qg/V load is used in preference to a per-cell thickness because siim does not resolve 3D ice flow: sheet accumulation Qg = area*balance is trustworthy, the unresolved per-cell H is not (see DECISIONS 2026-06-16).

The plate solve, densities, boundary handling and rebound feedback are inherited unchanged. The ice LOAD is applied everywhere ice is present, including afloat cells (zs = zb + hc*H < bl): the waterline-flotation gate gates glacial EROSION, not the flexural load. This is an accepted secondary seam — the ice mass is physically present regardless of grounding, and a true floating-shelf load would be reduced by the water it displaces; siim’s overdeepenings sit right at the flotation draft, so the residual is small. Hillslope diffusion of the ice surface is deliberately not in the load (an accepted small approximation: it is a near- conservative redistribution of the ice surface, not rock unloading).

class siim.fastscape.processes.HillslopeDiffusion

Bases: object

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

Same input/output contract as LinearDiffusion — reads SurfaceToErode.elevation + the diffusivity input, writes the per-step erosion into the "erosion" group — so it drops into the 'diffusion' model slot when numerics_backend='inhouse'. The unconditionally-stable ADI solve (siim._core.hillslope.diffuse()) is bit-for-bit with the fortran for the uniform kd siim uses. Boundary handling keys off the same ibc code (BorderBoundary), read directly rather than via the fortran context.

class siim.fastscape.processes.D8FlowRouter(*args: Any, **kwargs: Any)

Bases: process

In-house D8 single-flow router (S4) — the framework-free replacement for fastscape’s fortran SingleFlowRouter (fs.flowroutingsingleflowdirection). Fills the flow slot when routing='single' and router_backend='inhouse_d8'. Produces the SFR bundle (1D receivers / lengths, all-ones weights / nb_receivers, outlet-first stack) via siim._core.step.route_d8() on the eps-filled surface; basin labeled by outlet index. The routing delta vs the fortran SFR is confined to depression/tie cells (behavioral/attractor gate).

class siim.fastscape.processes.DinfFlowRouter(*args: Any, **kwargs: Any)

Bases: process

D-infinity flow router [Tar97].

Subclasses the in-house router shell — fully fortran-free (S4): the interior/boundary mask comes from border_status (Map 3 §4, provably identical to the old fortran sfr_rec != i mask) and basin from the in-house outlet labeling, so the last fs.flowrouting() call is gone. Outputs use the (n_nodes, 2) receiver/weight shape consumed by GlacialFlowAccumulator and the mode-B / mode-A D-inf kernels.

Continuous facet-direction weights replace the discrete neighbor weighting jumps of Quinn-style multi-flow routing, removing per-step receiver-flip kicks (~10x reduction at 100x100 vs Quinn MFR).

Depressions: directions are computed on the eps-filled surface (priority-flood from the outlet cells, _priority_flood_eps), so flux crosses closed basins toward their spills — the same depression semantics as the fortran SFR. The PHYSICS stays on the true surface: the mass balance samples real z (flux crossing a deep trough melts at the drowned elevation and dies) and the erosion/H kernels clip negative true-surface slopes to zero (lake interiors do not erode). Looped boundaries wrap (fill and facet scan both). See docs/guides/concepts.md for the public routing overview.

Forcing processes — siim.fastscape.forcing

Generic fastscape forcing processes.

Two reusable @xs.process forcing classes — not glacial, not siim-specific — that drop into any fastscape model. They came out of the escarpment work (siim.escarpment) but depend only on fastscape’s grid/boundary/surface processes:

• WaveUplift — a moving Gaussian uplift wave (a drop-in replacement for fastscape’s BlockUplift).

• PlateauSurface — an arctan-smoothed plateau initial topography (a drop-in replacement for an initial-elevation process).

class siim.fastscape.forcing.WaveUplift

Bases: object

Moving Gaussian uplift wave (replaces fastscape’s BlockUplift).

rate(x, t) = U_inf + U0 * exp(-(x - wave_center(t))^2 / wave_width^2), wave_center(t) = x[0] + x_escarpment + t * wave_velocity, U0 = wave_calibration * delta_h * wave_velocity / (wave_width * sqrt(pi)).

class siim.fastscape.forcing.PlateauSurface

Bases: object

Arctan-smoothed plateau initial topography (replaces InitialTopography).

Arctan ramp from exactly 0 on the low-x side to plateau_zo on the high-x side, centered at x_escarpment = (1 - plateau_frac) * Lx with transition width plateau_w (rescaled so the domain ends land on 0 / plateau_zo - plateau_dz). A small slope plateau_dz across the plateau seeds the divide. Uniform noise is added for D8 tie-breaking and zeroed on ‘fixed_value’ edges (matching InitialTopography, so the boundary z stays put across the run).