2D landscape model

The constructor accepts a parameter dictionary; see Numerical-model parameter reference for every key and mode-dependent default.

Wrapper / dispatch — siim.siim2d

class siim.siim2d.siim(user_params=None)

Bases: object

Standalone two-dimensional coupled glacial-fluvial landscape model.

Parameters:

user_params (dict, optional) – Parameter overrides. Omitted keys use the documented model defaults; unknown keys raise ValueError. See the public parameter reference for units, accepted values, and mode-specific behavior.

Notes

The default driver is the in-house NumPy/Numba implementation and does not require Fastscape or xsimlab. Call run() to integrate, use the *_out arrays or ds_out for scientific output, and use plot for visualization. save() and load() use versioned pickle files; load only files from trusted sources.

run(hooks=None, driver=None)

Run the model and unpack outputs. driver selects the time-loop backend: 'xsimlab' (the fastscape/xsimlab adapter orchestration) or 'inhouse' (siim’s own framework-free loop, siim._core.driver); default constants.DRIVER_DEFAULT. Both produce the identical ds_out (bit-for-bit on the same backend). hooks (an xsimlab RuntimeHook) is honored ONLY by the xsimlab driver — the standalone driver drops it (OQ-2).

save(filename=None)

Atomically pickle the current model state (user_params + ds_out) to ./model_outputs/saved_models/ under the current working directory. With no filename, an auto-name combining a UTC timestamp and run_id is generated. Returns the path written.

classmethod load(filename)

Rebuild a siim instance from a previously save()-ed pickle. filename may be a name under ./model_outputs/saved_models/ (with or without the .pkl suffix; a subfolder such as 'batch/run003' round-trips with save()) or an absolute / relative path.

Pickle files can execute code while loading. Only load files from a trusted source. Unversioned legacy payloads are rejected rather than interpreted under current model conventions.

extract_channel(i=-1, basin_rank=0)

Extract the main channel from a basin at output time step i, then follow those same nodes through all output time steps.

Parameters:

• i (int) – Reference time step for channel extraction (default -1, last step).

• basin_rank (int) – Which basin to extract from, ranked by node count. 0 = largest (default), 1 = second largest, etc.

• with (Returns a SimpleNamespace) – nodes — flat node indices along the channel (headwater → outlet) x_coord — x coordinates (m) along the channel y_coord — y coordinates (m) along the channel distance — cumulative distance from headwater (m) z, zb, H, Qg, Qf, area, erosion_rate — (nt_out, n_nodes) arrays k_h, d, xo, L — Hack’s law fit from time step i

strahler_order(i=-1, channel_threshold=1000.0)

Compute Strahler stream order over the full grid, threshold-gated.

Parameters:

• i (int) – Output time step index (default -1, last step).

• channel_threshold (float) – Minimum upstream area in cell units to be considered a channel node. Nodes below this threshold get order 0.

Returns:

orderndarray, shape (grid_ny, grid_nx)

Strahler order at each node (0 for sub-threshold non-channel nodes).

areandarray, shape (grid_ny, grid_nx)

Upstream catchment area (m²) at each node.

area_by_orderdict

Mean-area member arrays keyed by Strahler order (1..max_order).

max_orderint

Highest Strahler order present.

sigmandarray

Per-order area ratios log(2)/log(A_{k+1}/A_k) (Horton scaling).

Return type:

SimpleNamespace with

Notes

Works correctly with flow_routing=’single’ and ‘dinf’ (SFR-equivalent donor topology in both cases — D-inf’s primary-receiver collapse is meaningful because facet-weighted routing picks a clear dominant direction).

siim.siim2d.load(filename)

Module-level shortcut for siim.load().

filename must be a trusted SIIM pickle; pickle can execute code before the saved-state envelope is validated.

The standalone driver uses the fastscape-free numerical core for routing, ice-thickness closure, erosion, sub-grid carving, diffusion, and flexure; see Numerical core — siim._core. Optional xsimlab process shells around the same step functions are documented in fastscape integration — siim.fastscape.

Plotting — siim.plotting

siim plotting package — siim_plotter assembled from per-area mixins.

Split out of the old siim2d_plotting.py for file-size hygiene and a trimmed v1.0 keeper set; the public API (m.plot.<method>) is unchanged. See docs/guides/outputs_and_io.md for the public plotting contract.

siim_plotter is assembled from the per-area mixins below; siim2d.py imports it from here, and the old siim2d_plotting.py is gone.

class siim.plotting.siim_plotter(model)

Bases: MapMixin, ProfileMixin, LandscapeMixin, BasinsMixin, DiagnosticsMixin

Plotting facade bound to a run model instance (m.plot).

The 1D facade remains in siim1d.py. Both inherit the same public profile methods and share field metadata, visual styles and serial movie handling. The 1D limit-cycle and phase-portrait helpers are private research tools, outside the public plotting API.

Raw 2D fields — siim.plotting.maps

map / view_map / animate_map — a straight 2D raster of a stored model field (bedrock, ice, …).

Registry-driven: adding a field is one entry in FIELD_REGISTRY (the data is always stored). No smoothing / hillshade / sub-grid ice — that is landscape. See docs/guides/outputs_and_io.md for the public plotting contract.

(“field” here is which stored variable to draw — the field= argument and FIELD_REGISTRY — distinct from the map method that draws it.)

class siim.plotting.maps.MapMixin

Bases: object

Straight raster of a stored 2D field: map / view_map / animate_map (frame / interactive / mp4 over the same render).

map(field='bedrock', i=-1, field_min=None, field_max=None, cmap=None, ax=None, fig_width=6)

Stored node-valued raster at snapshot i; returns the Axes.

Ice is width-mean thickness, with ice-free ground masked. Color limits cover all outputs unless field_min/field_max are supplied.

view_map(field='bedrock', field_min=None, field_max=None, cmap=None, fig_width=6)

Notebook slider; update the same raster and colorbar in place.

animate_map(field='bedrock', path=None, run_id=None, field_min=None, field_max=None, cmap=None, fps=20, interval=42, fig_width=6, frames=None)

MP4 of map(). Explicit fps wins; fps=None uses interval (ms/frame).

frames encodes a subset of the saved frames: a slice, a range or a sequence of indices (negatives count from the end).

Channel profiles & cross-sections — siim.plotting.profiles

profile / view_profile / animate_profile — the extracted-channel profile vs the analytical steady state, plus view_xsection (raster cross-section) and the analytical overlay helper.

Unified call shape (identical in 1D and 2D):

profile(fields=’elevation’, i=-1, field_min=None, field_max=None,

basin_rank=0, ref=-1, analytical=True, bistable=True, ax=None)

fields is a str, a list (stacked panels), or a {field: (min, max)|None} dict; None (or a None tuple side) means auto = the global, all-frame stable limit (so the axis is constant across a view/animation). field_min / field_max are the single-field shorthand. 5 fields: elevation, ice_thickness, erosion_rate, ice_flux, water_flux. See docs/guides/outputs_and_io.md for the public plotting and output contract.

class siim.plotting.profiles.ProfileMixin

Bases: ProfileMethods

2D channel extraction and references for the common profile renderer.

analytical(ax=None, bistable=True)

Plot the analytical SS reference profile (distance axis). Assumes a prior profile/channel extraction set the analytical grid.

view_xsection(fields='elevation', figsize=None, y_km=None, z_min=None, z_max=None, H_max=None, fig_width=8, aspect=0.38)

Interactive x-direction cross-section viewer: two sliders (time, y-row). The cross-section is the raster row at the chosen y — no flow graph / analytical. fields ⊂ {elevation, ice_thickness, erosion_rate}. Uses ipympl for a smooth live canvas (managed locally, so your other plots stay on the default backend).

Landscape rendering — siim.plotting.landscape

landscape / animate_landscape — the deluxe atlas-style cartographic render: supersampled + Gaussian-smoothed terrain, sub-grid across-channel-width ice (footprint fill — the display dual of the width carve), priority-flooded lakes, hillshade, contours, and an optional cross-section + hypsometry panel. Frame + animate only (no viewer, by design). See docs/guides/outputs_and_io.md for the public plotting contract.

field ∈ {'bedrock', 'bedrock+ice', 'bedrock+lakes'}.

The numba render backend (footprint ice, smoothing, hillshade, path tracing) lives in _render.py.

class siim.plotting.landscape.LandscapeMixin

Bases: object

Deluxe cartographic renderer: landscape + animate_landscape.

landscape(field=None, i=-1, style='smooth', oversample=None, channel_threshold=100000.0, H_threshold=None, area_threshold=None, z_min=0.0, z_max=None, cmap_bed='gist_earth', ice_extent=None, ice_color='#e5e8ed', ice_shading=None, trunk_display=None, trunk_width_cells=1.0, trunk_alpha=0.85, ice_smoothing='mask', min_ice_cells=0, ice_time_avg=None, ice_cmap=None, H_min=0.0, H_max=None, lake_color='#a8c8e8', lake_min_depth=0.5, lake_min_area=0.0, sigma_cells=None, ice_sigma_cells=None, hillshade=None, azdeg=315, altdeg=45, ve=2.0, contour_interval=None, contour_color='k', contour_lw=0.2, contour_alpha=0.4, show_margin=None, margin_color='#1b4f72', margin_lw=1.0, show_trimline=None, trimline_color=None, trimline_lw=None, cross_section=None, cross_section_color='#716477', hyp_max=15.0, show_smoothed_paths=False, smoothed_path_color='red', smoothed_path_lw=0.4, fig_width=8, fig=None, ax=None, ax_cs=None, ax_hyp=None, colorbar=True, save=None)

Cartographic view of the stored bed, ice or lake display layers.

Returns (fig, ax). i selects the output snapshot (default last). Display processing never changes stored model arrays.

Presets and layers

style='smooth' draws bedrock + ice with oversample=4, sigma_cells=oversample, hillshade, 100 m contours, footprint ice, depth-graded translucent shading, resolved trunk ribbons and margins. style='raw' uses oversample=1, unsmoothed channel cells, flat ice, no contours/margins/hillshade/ribbons and area_threshold=1e6 m². Both use H_threshold=0 and ice_time_avg=1. An explicit argument always overrides the preset; None selects the preset’s value.

field selects ‘bedrock’, ‘bedrock+ice’ (default), or ‘bedrock+lakes’. Bedrock is the stored true bed zb_out. Ice is a display reconstruction; lakes are priority-flooded on the true bed. These are separate recipes, not names of stored quantities.

Ice extent and filtering

ice_extent='footprint' fills the claimed width W=alpha_g*H; ‘cells’ shows the ice-bearing channel cells. trunk_display='ribbons' traces ice downstream of cells spanning trunk_width_cells native cells; ‘none’ disables this extra layer. Ribbons carry their source ice surface into hillshade and the cross-section.

H_threshold is a gate on WIDTH-MEAN thickness H (m), and area_threshold gates upstream area (m²). Both govern the source cells for footprints and ribbons; they do not clip the width grown from a surviving source. Zero disables each gate. The raw area’s default can be disabled to inspect every ice-bearing cell.

ice_sigma_cells smooths the ice mask in SUBGRID pixels (smooth default 2, raw 0). Its native-grid strength is divided by oversample. ice_smoothing='mask' smooths a binary indicator; ‘field’ smooths clipped thickness min(H, 2*H_threshold) before thresholding. The field option requires extent=’cells’ and H_threshold>0. min_ice_cells removes small connected components and fills small enclosed holes in the final mask, including ribbons. Its unit is native cells; it respects looped seams and defaults to 0 (off). ice_time_avg averages this many trailing OUTPUT snapshots for the ice display only. Terrain, sections and hypsometry retain snapshot i. Averaging and component removal are opt-in because they can obscure real changes or small glaciers.

Colors and terrain

ice_shading='veil' blends the depth ramp over terrain; ‘flat’ uses ice_color. ice_cmap overrides either with opaque column-depth coloring. H_min/H_max bound COLUMN depth, not H. The veil’s automatic maximum is run-global hc_over_H*max(H_out); an explicit ice_cmap otherwise scales its maximum per snapshot. trunk_alpha is the minimum ribbon opacity under veil shading. lake_color, lake_min_depth (m), lake_min_area (m²) govern lakes.

z_min/z_max and cmap_bed set the terrain color scale. oversample sets subgrid resolution. sigma_cells is terrain Gaussian smoothing in subgrid pixels; 0 disables smoothing. hillshade, azdeg/altdeg (degrees) and ve set lighting. contour_interval is in meters, anchored at zero; 0 disables it. contour_color/contour_lw/contour_alpha style those contours.

Overlays and composition

show_margin, margin_color and margin_lw control the current ice outline (ribbons alone when present). The deprecated aliases show_trimline, trimline_color and trimline_lw still work. show_smoothed_paths adds a diagnostic centerline overlay, gated by channel_threshold (upstream area in m²), with smoothed_path_color/smoothed_path_lw styling. cross_section is y in km, or 'mean' to average over the rendered rows of the section grid (an interpolated oversample grid, so the two end rows carry less ground than the interior ones) — an averaged profile whose ice band is the average ice column per unit x, with the ice thickness’s 25-75% band across those rows shaded on the mean bed and the bed’s quartiles dashed, both named in a legend, and with no map locator line and no lake layer (a reduced water table is not a water table). Either way it adds a section and bed hypsometry. cross_section_color styles its map locator. hyp_max fixes the percent-area scale (default 15); None auto-scales it.

fig_width is in inches. fig and supplied ax/ax_cs/ ax_hyp must belong to one figure. With a supplied map ax, supply ax_cs (and optionally ax_hyp) for section panels. colorbar=False permits caller-managed shared bars. save names the image output. See docs/guides/outputs_and_io.md for examples and output paths.

animate_landscape(path='landscape_animate', run_id=None, *, interval=42, fps=None, frames=None, fig_width=8, n_workers=None, **landscape_kwargs)

MP4 of landscape (its kwargs passed through). Returns the path.

interval is milliseconds per frame; an explicit fps overrides it on both render paths. frames encodes a subset of the saved frames — a slice, a range or a sequence of indices, negatives counting from the end — in the order given.

Frames are independent, so they render IN PARALLEL by default (audit N35): the arrays landscape reads are dumped once and loaded by n_workers spawned processes — ONLY those five arrays plus a small metadata pickle, never the full run state, and the worker count is additionally capped so the workers’ combined array memory stays within ~25% of system RAM. Each worker renders its share of frames to PNG via the unchanged landscape(); ffmpeg assembles the MP4. n_workers=None auto-sizes (up to 8, capped by cores, frame count and the RAM budget, and only when a timed probe frame projects the serial render past ~30 s); n_workers=1 is the original in-process serial render (also the automatic fallback when fig/ax are passed or the model lacks the arrays the workers need). Frame content is identical either way — only the wall clock changes.

The auto colour scales are FROZEN over the movie: z_max and the ice norm’s H_max are resolved once across the ENCODED frames instead of per frame, so a frame’s colours mean the same thing throughout the movie; pass explicit values to override. (Over the whole run, H_max freezes to the run-global hc_over_H * max H_out, which is already the still render’s default for the depth-graded veil; a frames window scales to that window, not to relief it never shows.)

Inherits landscape’s style='smooth' default (the cartographic view — supersampled + hillshaded terrain, footprint-width ice drawn as a depth-graded translucent veil, with nothing hidden by thickness — H_threshold=0). Pass style='raw' for the naked model output (one pixel per cell, no hillshade or smoothing, flat ice cells shown directly), or override individual knobs (explicit values always win over the preset). For DENSELY-SAMPLED output (frames close in model time), the classic anti-flicker movie recipe is ice_sigma_cells=3, ice_time_avg=2 — deliberately NOT defaulted: ice_time_avg blends output frames (a double exposure of two glacial epochs when saves are far apart) and larger ice_sigma globs and erases thin ice. min_ice_cells=6 drops specks but HIDES real small glacierets — deliberately in no preset. Keep ice_extent='footprint' (the default) for the width-honest view (ice across the claimed width W = alpha_g*H); 'cells' shows the raw channel-cell state. See landscape() for the full knob reference.

Basin & channel diagnostics — siim.plotting.basins

Public basin plots, with computation delegated to _basin_data.

class siim.plotting.basins.BasinsMixin

Bases: object

Channel/basin analysis methods.

flux(which='ice', i=-1, ax=None)

Single log-scale map of a flux field at step i.

which ∈ {'ice', 'water', 'area'}.

hacks_law(i=-1, *, ref=-1, basin_rank=0, plot=True, verbose=False, ax=None, fig_width=8)

Main-channel area and its Hack’s-law fit.

i selects the area snapshot; ref selects the fitted reference channel and defaults to the final output. Returns an Axes, or a data namespace with plot=False. Scatter/curve coordinates in that namespace are named distance_km/area_km2 and fit_distance_km/fit_area_km2; fit coefficients retain SI units. Only verbose=True prints the fit summary.

largest_basins(n_basins=4, i=-1, channel_threshold=100000.0, plot=True, z_max=None, *, ax=None, fig_width=10)

Extract the largest basins and optionally draw a map/profile pair.

i selects both the basin ranking and the displayed snapshot. Returns a namespace (channels, k_h, d, L, zo, xt, sigma), also when plot=False. Distances and elevations are meters. ax accepts a pair of axes; fig_width controls newly created figures in inches.

largest_basins_history(n_basins=4, ref=-1, t_start=0, t_end=None, n_samples=None, channel_threshold=100000.0, plot=True, *, i_ref=None, verbose=False, ax=None, fig_width=10)

Trace reference outlets through an inclusive output-step window.

ref fixes basin identity, matching profile(ref=...). i_ref is a deprecated keyword alias. t_start/t_end are output-step indices, not years; negative indices count from the final output. Returns a namespace of (basin, sample) arrays and stats. Missing channels/fits are NaN; stats omit them and include valid counts in stats[metric]['n']. plot=False creates no figure. Printing is opt-in with verbose=True. ax accepts nine axes in row order.

sediment_history(n_basins=4, ref=-1, quantity='volume', ax=None, *, i_ref=None, plot=True, fig_width=10)

Outlet sediment output for basins selected at ref.

Needs track_sediment=True. quantity is 'volume' (cumulative km³) or 'flux' (km³/yr). Returns a namespace with t (years), outlets, series and quantity, including with plot=False. The default figure has a map and time series; ax accepts either a single time-series axis or a map/series pair. i_ref is deprecated.

Steady-state diagnostics — siim.plotting.diagnostics

Steady-state measurements and their public plotting wrapper.

class siim.plotting.diagnostics.DiagnosticsMixin

Bases: object

Convergence diagnostics.

steady_state(tol_E=0.05, tol_z=0.1, window=5, z_max=None, *, plot=True, verbose=False, ax=None, fig_width=10)

First output where erosion and elevation trend satisfy the tolerances.

Spatial-mean erosion must match mean uplift within tol_E (relative), and the rolling mean-elevation slope must be below tol_z * abs(U). Fixed-value boundaries are excluded from spatial means. Returns the first matching output-step index, or -1 (including single-output runs). plot=False creates no figure; printing requires verbose=True. ax accepts a pair of axes for the elevation and erosion panels.

Escarpment variant — siim.escarpment

Escarpment-specialized siim2d (thin convenience wrapper).

Rides on top of the standard siim2d.siim model, wiring two generic fastscape forcing processes — WaveUplift and PlateauSurface, which now live in siim.fastscape.forcing (reusable by any fastscape model) — in through the base class extension seams (_default_params / _process_overrides / _forcing_input_vars):

• WaveUplift — a moving Gaussian uplift wave (replaces fastscape’s BlockUplift), for a retreating-escarpment / passing uplift-pulse setting.

• PlateauSurface — an arctan-smoothed plateau initial topography (replaces InitialTopography), for starting from a high plateau with a sharp escarpment edge.

Everything else (solvers, routing, plotting, analytical, save/load, channel extraction, …) is inherited unchanged from siim2d.siim, so improvements to the core model flow through to the escarpment model automatically.

Usage:

from siim.escarpment import siim_escarpment
m = siim_escarpment({
    "uplift_type": "wave", "delta_h": 1500, "wave_velocity": 15e-3,
    "init_type": "plateau", "plateau_zo": 1500,
    "sliding_law": "coulomb", ...
})
m.run()
m.plot.landscape()

siim.escarpment.plateau_topography(nx, ny, Lx, zo=1000.0, frac=0.8, w=10000.0, dz=1.0)

Plateau initial surface for plain siim2d runs, with its own base level.

Returns (topo, bl): the (ny, nx) surface — exactly 0 at the left edge, rising to zo across an escarpment centred at x = (1 - frac)*Lx with transition width w (m), then dropping dz linearly to exactly zo - dz at the right edge — and the per-side base level {'right': zo - dz} to pass as the bl parameter. The right outlet’s water line is then the plateau’s own edge elevation, so every basin draining onto the plateau is graded to plateau height and base level does not incise the plateau from behind; the left outlet keeps the default datum. Same profile as siim_escarpment(init_type='plateau') (siim._core.step.plateau_profile()); siim2d adds its own noise_amplitude (zeroed on fixed borders) on top. Usage:

params['initial_topography'], params['bl'] = plateau_topography(
    nx, ny, Lx, zo=2000, frac=0.85, w=25e3, dz=1000)

class siim.escarpment.siim_escarpment(user_params=None)

Bases: siim

siim2d specialized for escarpment problems.

Accepts the full standard siim2d parameter set plus:

uplift (uplift_type='wave'):
  uplift_type      'block' (default, standard BlockUplift) or 'wave'
  delta_h          target integrated uplift as the wave passes (m); required for 'wave'
  wave_width       1/e half-width of the wave (m)
  wave_velocity    wave propagation velocity (m/yr)
  x_escarpment     initial wave-center position relative to the left edge (m)
  wave_calibration calibration factor on the peak rate (1.0 = exact:
                   the wave deposits delta_h as it passes)
  U_inf            steady background uplift (m/yr)

initial topography (init_type='plateau'):
  init_type        'sloped' (default, standard InitialTopography) or 'plateau'
  plateau_zo       plateau elevation (m); required for 'plateau'
  plateau_dz       slope across the plateau (m); seeds the divide
  plateau_frac     fraction of x occupied by the plateau
  plateau_w        escarpment transition width (m)

The two switches are independent: you can mix e.g. a sloped start with a wave uplift, or a plateau start with block uplift.