Constants and support modules

Package entry — siim

siim — Sliding Ice Incision Model.

Coupled glacial-fluvial landscape evolution: a 1D profile model (siim.siim1d), a standalone 2D raster model (siim.siim2d, with an escarpment variant in siim.escarpment), and the analytical steady-state solution (siim.analytical) both numerical models embed as their reference.

Entry points:

from siim.siim1d import siim                       # 1D profile model
from siim.siim2d import siim                       # 2D landscape model
from siim.fastscape import glacial_processes       # optional adapter for your own fastscape model
from siim import analytical_steady_state_solution  # analytical reference
from siim.analytical import GeneralProfile         # steepness-form theory
from siim.analytical import MarginalCoulombProfile # exact d*phi=1 case
from siim.forcing import ela_sawtooth, uplift_step, interp_forcing # time-varying forcing builders

The 1D and 2D models run directly on the package’s NumPy/Numba numerical core; the 2D driver uses in-house routing, flexure, and diffusion. xsimlab, fastscape, and fastscapelib-fortran are needed only for the optional siim.fastscape adapter. siim.analytical itself needs only NumPy/SciPy. This module imports nothing eagerly so the optional adapter never affects a standard import.

Constants — siim.constants

Shared physical constants and per-sliding-law parameter dispatch.

Single source of truth for the numerical models (siim1d, siim2d, the standalone fastscape process module) and the analytical machinery (siim.analytical). Imports nothing heavier than numpy, so the lightweight analytical layer can depend on it without dragging in the model stack.

siim.constants.cg_prefactor(alpha_g=5.0, Ac=2.5e-24, rho_ice=920.0, g=9.8)[source]

Glacial rheology group cg = alpha_g * kt * (2*Ac/5) * (rho_ice*g)^3 [m^-3 yr^-1] (the s/yr conversion kt is absorbed). The 2/5 is the depth-integrated Glen deformation prefactor 2A/(n_c+2) with n_c = 3 (model-paper convention; Ac is the physical Glen coefficient).

class siim.constants.LawConstants(ell, nu, mu, phi, Co)
Co

Alias for field number 4

ell

Alias for field number 0

mu

Alias for field number 2

nu

Alias for field number 1

phi

Alias for field number 3

siim.constants.derive_coulomb(ce, alpha_g, tau_c, rho_g, g, *, nu=None, ell=None)[source]

Coulomb law: nu = 2*ell, mu = ell. Co = ce * ((rho_g g)^2 / (alpha_g tau_c^2))^ell.

siim.constants.Co_power(ce, cg, lambda_p, alpha_g, mu)[source]

Power/eff-exp erosion prefactor Co = ce * (cg^(5/4) * lambda_p^3 / alpha_g^(9/4))^mu. Single-sourced so an explicit mu override can recompute Co without re-literalling the formula (siim1d / analytical.steady_state, audit B5).

siim.constants.derive_power(ce, cg, lambda_p, alpha_g, *, nu=None, ell=None)[source]

Power and eff-exp laws: nu = 5*ell/3, mu = 4*ell/9. Co = ce * (cg^(5/4) * lambda_p^3 / alpha_g^(9/4))^mu.

siim.constants.mode_c_standard(is_mode_c)[source]

The mode-C standard defaults (trunk_surface, routing_relax).

(MODE_C_TRUNK_SURFACE, MODE_C_ROUTING_RELAX) for mode C (mode B with sub-grid carving), (False, ROUTING_RELAX) otherwise. Callers apply each value only to a sentinel left at None; an explicit user value wins.

siim.constants.normalize_mode(value)[source]

Map a user-facing surface-evolution mode to its internal code.

Accepts the legible names 'ice_surface' (one tracked state, the ice surface) and 'bedrock+ice_thickness' (two tracked states, bedrock + ice thickness), plus the permanent short aliases 'A' / 'B' (exact match) and 'C' / 'c' (the mode-B + sub-grid-carve alias, resolved to 'B' + carving by siim2d). Returns 'A', 'B' or 'C'.

siim.constants.widening_factor_from_rate(value)[source]

User-facing widening_rate (eta >= 0) -> internal E_widening/E_c factor.

Returns 1 + eta, the multiple of the centerline erosion rate the carve kernel applies to a footprint cell. None / inf / 'inf' / 'infinity' (any case) request instant U-imposition (factor = inf). Negative eta is rejected – disable carving with carve_width=False.

Time-varying forcing builders — siim.forcing

Time-series builders for time-varying model forcing.

Each function builds the run’s time vector internally from T (run time, years) and nt (number of steps) and returns (t, series): the time vector (handy for plotting) and the length-nt forcing array that drops straight into a model parameter dict, e.g.:

from siim.forcing import ela_sawtooth
t, zELA = ela_sawtooth(params['T'], params['nt'], ela_high=2400, ela_low=1400)
params['zELA'] = zELA

Both siim.siim1d and siim.siim2d accept a scalar or a length-nt array for zELA, U, and P. numpy-only (no model import), so this stays a cheap import.

siim.forcing.ela_sawtooth(T, nt, ela_high=1500, ela_low=300, period=100000.0, buildup_frac=0.88)[source]

Asymmetric sawtooth ELA(t): slow glacial buildup, fast termination.

Builds t = linspace(0, T, nt) internally. Linear ramp down from ela_high to ela_low over buildup_frac of each period (the slow drop into a glacial), then linear ramp back up over the remaining 1 - buildup_frac (the fast termination).

Parameters:

  • T (float) – Total run time in years (the time vector spans 0 to T).

  • nt (int) – Number of time steps (length of the returned array).

  • ela_high (float) – ELA (m) at the start of buildup and at termination onset.

  • ela_low (float) – ELA (m) at the start of buildup and at termination onset.

  • period (float) – Cycle length in years.

  • buildup_frac (float) – Fraction of the period spent in the slow buildup (0 < frac < 1).

Returns:

(t, ela) – The time vector and the ELA(t) series, each length nt.

Return type:

tuple of ndarray

siim.forcing.uplift_step(T, nt, U_init, U_final, step_frac=0.5)[source]

Step change in uplift rate partway through the run.

Builds t = linspace(0, T, nt) internally and returns U_init before the step and U_final at and after it. The step falls at step_frac * T (so step_frac is a fraction of the run time).

Parameters:

  • T (float) – Total run time in years (the time vector spans 0 to T).

  • nt (int) – Number of time steps (length of the returned array).

  • U_init (float) – Uplift rate (m/yr) before and after the step.

  • U_final (float) – Uplift rate (m/yr) before and after the step.

  • step_frac (float) – Where the step lands, as a fraction of the run time (0 to 1).

Returns:

(t, U) – The time vector and the uplift-rate series, each length nt.

Return type:

tuple of ndarray

siim.forcing.interp_forcing(T, nt, times, values, left=None, right=None)[source]

Piecewise-linear forcing series from coarse (times, values) nodes.

Interpolates values (defined at times, in run-time years from 0) onto t = linspace(0, T, nt). A generic table-to-series builder, reusable for any scalar forcing — P, zELA, or U:

from siim.forcing import interp_forcing
# falling precipitation over the run (m/yr)
t, P = interp_forcing(params['T'], params['nt'],
                      times=[0, 1.5e6, 3e6], values=[2.0, 1.2, 0.65])
params['P'] = P

times are model-time years (0 → T); map geological time (e.g. Ma) onto that axis before calling, the same way you would for zELA.

Parameters:

  • T (float) – Total run time in years (the time vector spans 0 to T).

  • nt (int) – Number of time steps (length of the returned array).

  • times (array_like) – Node positions (run-time years) and the forcing value at each node. times must be increasing (np.interp requirement).

  • values (array_like) – Node positions (run-time years) and the forcing value at each node. times must be increasing (np.interp requirement).

  • left (float, optional) – Value returned below times[0] / above times[-1]. Defaults to the nearest endpoint (flat extrapolation, np.interp default).

  • right (float, optional) – Value returned below times[0] / above times[-1]. Defaults to the nearest endpoint (flat extrapolation, np.interp default).

Returns:

(t, series) – The time vector and the interpolated forcing series, each length nt.

Return type:

tuple of ndarray