Analytical steady state

Analytical steady-state machinery for the coupled glacial-fluvial profile.

Lightweight by design: importing this subpackage needs numpy/scipy only (no fastscape/xsimlab, no matplotlib until plot is called), so the theory classes can be used for paper figures without the model stack.

Entry points:

from siim.analytical import GeneralProfile          # general exponents
from siim.analytical import MarginalCoulombProfile  # exact d*phi = 1 case
from siim.analytical import SteadyStateProfile      # physical parameters

GeneralProfile is the general steady state in steepness form (ks, cs, zELA, L plus exponents); MarginalCoulombProfile is the closed-form (arcosh) solution of the marginal-Coulomb special case and doubles as the cross-check oracle for the general machinery; SteadyStateProfile (alias analytical_steady_state_solution) maps physical model parameters onto the same engine. RegimeMap evaluates the same closure over (kappa = Ng/Nf, Y = zELA/zfo) grids for regime diagrams, including the saddle-node bistability boundary. All share the Solution / AARResult result types, the shared closure solver, and the incomplete_beta kernel below.

Math kernel — siim.analytical.core

Math kernel for the analytical steady-state machinery.

Pure numpy/scipy: the analytically continued incomplete beta (including the exact logarithmic limit at b = 0, i.e. the marginal-Coulomb case d*phi = 1), the closure root-solver shared by the profile classes, and the shared result containers. Nothing here imports the model stack.

siim.analytical.core.incomplete_beta(x, a, b)[source]

Non-regularized incomplete beta function, analytically continued in b.

\[B(x;\,a,\,b) = \int_0^x t^{a-1} (1-t)^{b-1}\,dt\]

Valid for a > 0 and any real b, including b = 0 (the logarithmically divergent marginal case d*phi = 1, handled exactly) and non-positive integers. x in [0, 1]; at x = 1 the value is the complete beta for b > 0 and +inf for b <= 0. Vectorizes over x.

When x is close to 1 and you know eps = 1 - x exactly (e.g. (xo/Lt)^k), call incomplete_beta_compl() instead to avoid the 1 - x roundoff.

Parameters:

  • x (float or array_like) – Upper integration limit(s), in [0, 1].

  • a (float) – First shape parameter; must be > 0.

  • b (float) – Second shape parameter; any real value (analytic continuation).

Returns:

B(x; a, b), scalar when x is scalar, else an array of x’s shape.

Return type:

float or ndarray

Raises:

ValueError – If a <= 0, if x is outside [0, 1], or if b is so negative that the b -> b+1 recurrence would exceed the internal step cap.

siim.analytical.core.incomplete_beta_compl(eps, a, b)[source]

B(1 - eps; a, b) with eps supplied directly (exact for tiny eps).

The profile shape function G evaluates the incomplete beta at 1 - (x/Lt)^k; for x << Lt the argument is so close to 1 that forming it in floating point loses (or zeroes) the information. Passing eps = (x/Lt)^k here keeps the evaluation exact all the way into the divergent tail — this replaces the old upper-limit clip, which silently froze G once (x/Lt)^k dropped below 1e-12. Vectorizes over eps.

Parameters:
Returns:

B(1 - eps; a, b), scalar when eps is scalar.

Return type:

float or ndarray

Raises:

ValueError – Same conditions as incomplete_beta().

siim.analytical.core.find_closure_roots(f, lo, hi, *, kink=None, n=4000, hi_max=None, n_per_decade_ext=400, refine=64)[source]

All roots of a vectorized scalar function f on [lo, hi(, hi_max)].

Log-spaced scan + brentq on every strict sign change. Grid nodes where f is exactly zero count as roots once (np.sign(0) double-bracketing is avoided). Local |f| minima without a sign change are re-scanned refine times finer to catch fold pairs narrower than one cell. If hi_max is given, the grid is extended decade-by-decade past hi (coarser, n_per_decade_ext points per decade) so far-out roots are not missed. kink (e.g. Lt = L) is inserted as an exact grid node so no bracket straddles the non-smooth point. brentq failures on a valid bracket propagate — a sign change that cannot be refined is a bug in f, not noise to be swallowed.

Returns:

Sorted roots (duplicates within rtol 1e-9 merged).

Return type:

list of float

siim.analytical.core.closure_slope(f, root, *, kink=None, lo=None, rel_h=0.0001)[source]

One-sided secant estimate of df/dLt at a root.

Sided away from kink when the root is within 2h of it, so the estimate never straddles a point where f is continuous but not differentiable (at the mixed/glacial seam the glacial side can have a divergent derivative, which would otherwise invert the stability flag).

class siim.analytical.core.Solution(regime: str, Lt: float, zt: float, zo: float, stable: bool)[source]

Bases: object

A single steady-state solution.

stable reflects the sign of d(zELA)/d(Lt) at the root: a glaciated state is stable when cooling (lowering zELA) grows the glacier, i.e. d(zELA)/d(Lt) < 0. The warm saddle of a fold has the opposite sign. The fluvial state, when viable, is always stable.

class siim.analytical.core.AARResult(surface: str, aar: float, u_ela: float, fc: float, N: float, Aa_over_Ac: float, eta_bar: float)[source]

Bases: object

Steady-state accumulation-area ratio and its ingredients.

Two surface closures fix the ELA crossing u_ELA = x_ELA / Lt and hence the catchment accumulation fraction fc = u_ELA^d (see GeneralProfile.aar):

  • surface='powerlaw': u_ELA^k = lam -> fc = lam^(d/k)

  • surface='gsurface' : G(x_ELA) = (1-lam) G(xo) solved on the consistent glacial shape -> fc = u_ELA^d

For a fully glacial state the ice footprint is truncated at base level (x = L): the ablation integral stops there, and a crossing at or beyond L means the whole on-orogen glacier accumulates (aar = 1).

eta_bar is the mean glacierized fraction of the below-ELA swath (A_a divided by the swath catchment area). Ice cannot be wider than its valley, so eta_bar <= 1 physically; a value above 1 means the ice-width closure (alpha_g * kH) and the Hack closure (kh) are mutually inconsistent for these parameters and the set should be rejected — aar warns when this happens.

Profile classes — siim.analytical.profiles

Steady-state glacial-fluvial profile classes (steepness-index form).

Driven by ks, cs, zELA and L directly, not by the underlying erosion/sliding parameters — uplift, erodibility, sliding law and mass balance are absorbed into ks and cs.

GeneralProfile solves the full coupled steady state for arbitrary exponents (d, sigma, phi, theta, k). It computes the elevation profile z(x), ice thickness H(x), reports the regime (fluvial / mixed / glacial), and detects bistability by finding all roots of the continuous zELA(Lt) closure.

MarginalCoulombProfile solves the exact phi = theta = 1/2, d = 2 special case in closed (arcosh) form. That marginal case (d*phi = 1) is precisely where the general beta-function machinery is logarithmically divergent, so this class doubles as the closed-form cross-check of GeneralProfile: the two must agree in that limit.

Both classes share one Solution/AARResult vocabulary (siim.analytical .core), the same constructor convention (Class(params_dict) == Class(**params_dict)), the same argument order, and the package-wide constants from siim.constants. matplotlib is imported lazily inside plot; importing this module needs numpy/scipy only.

siim.analytical.profiles.draw_profile(ax, x, surface, bed, *, ice_x_range=None, topo_color='k', ice_color='#b0d6e2', ice_alpha=0.55, bedrock_color='#888888', bedrock_fade='auto', bedrock_alpha=0.6, show_bedrock=True, show_thickness=True, label=None)[source]

Render a steady-state profile in the paper’s canonical visual language: a depth-faded bedrock fill below bed, the ice column shaded between bed and surface, the topographic surface line, and the sub-ice bed line. The single source of truth behind both GeneralProfile.plot and analytical_steady_state_solution.plot.

x, surface, bed are in metres (x strictly increasing; plotted in km). bed is the channel-floor datum (surface - HC_OVER_H*H under ice, surface where ice-free). ice_x_range = (x_lo, x_hi) in metres bounds the sub-ice bed line — H -> 0 at the glacier endpoints, so a geometric extent draws cleaner than an H > 0 mask; None draws it wherever surface > bed. show_thickness is a no-op on a fully fluvial profile (no cells with surface > bed).

class siim.analytical.profiles.GeneralProfile(ks=None, cs=None, zELA=None, L=None, *, xo=None, d=2.0, sigma=0.5, phi=0.5, theta=0.5, k=1.0, lam=None, kH=0.05, gamma=None, alpha_g=5.0, kh=5.0)[source]

Bases: object

General steady-state coupled glacial-fluvial profile (steepness-in).

Parameters:

  • ks (float) – Fluvial and glacial steepness indices. ks has units m^(d*theta); cs has units m^(1-r). They roll up uplift, erodibility and climate.

  • cs (float) – Fluvial and glacial steepness indices. ks has units m^(d*theta); cs has units m^(1-r). They roll up uplift, erodibility and climate.

  • zELA (float) – Equilibrium-line altitude [m].

  • L (float) – Orogen half-width / base-level distance [m].

  • xo (float, optional) – Hillslope / channel-head cutoff [m]. Defaults to L/1000. Enters G(xo) in both closures and the reference relief zfo.

  • d (float, optional) – Hack-area exponent (A = kh * x^d). Default 2.0.

  • sigma (float, optional) – Hack-flux exponent for ice-flux integration. Enters only via lam = d*sigma/(d*sigma + k). Default 0.5.

  • phi (float, optional) – Glacial concavity index (mu/nu). Coulomb -> 1/2, power-law -> 4/15. Default 0.5.

  • theta (float, optional) – Fluvial concavity index (m/n). Default 0.5.

  • k (float, optional) – Ice-accumulation profile shape exponent (Beff ~ 1 - (x/Lt)^k). Default 1.0 (the linear-ansatz baseline).

  • lam (float, optional) – AAR-like ratio. Defaults to d*sigma/(d*sigma + k); override to decouple.

  • kH (float, optional) – Ice-thickness prefactor in H = kH * Lt^(g(1+d)/(1+phi)) * shape^g. Model-agnostic; tune to taste. Default 0.05.

  • gamma (float, optional) – Ice-thickness shape exponent. Defaults to phi (Coulomb closure); power-law uses (2/3)(1/3 + phi).

  • alpha_g (float, optional) – Valley width-to-thickness ratio (W = alpha_g H). Sets the AAR ablation area; only the group alpha_g/kh matters. Default constants.ALPHA_G (currently 5).

  • kh (float, optional) – Hack coefficient (catchment area A = kh x^d). Sets the AAR accumulation area. Default constants.KH (5).

solutions

All steady states found (stable branches + warm saddles).

Type:

list[Solution]

primary

Largest-Lt stable solution (fluvial Lt = NaN ranks below any glacier). Drives profile, thickness, plot and the convenience attrs.

Type:

Solution

bistable

True when more than one stable steady state exists.

Type:

bool

solve()[source]

Re-derive every downstream quantity from the current raw parameter attributes and re-run the steady-state solve. Called by __init__; call it again after mutating a parameter attribute (e.g. self.zELA) to refresh the object (construct-then-mutate-then-solve(); idempotent; audit m27).

property Nf

Fluvial erosion number.

\[N_f = \frac{k_s\, L^{1-d\theta}}{z_{\rm ELA}}\]

+inf at zELA = 0 (ELA at base level).

property Ng

Glacial erosion number.

\[N_g = \frac{c_s\, L^{r}}{z_{\rm ELA}}, \qquad r = \frac{1 - d\phi}{1 + \phi}\]

+inf at zELA = 0.

property kappa

Master control ratio kappa = Ng / Nf = cs L^r / (ks L^(1-d*theta)). Equals cs/ks only in the marginal-Coulomb case (where r = 1-d*theta = 0). zELA cancels in the ratio, so kappa stays finite at zELA = 0.

property kappa_c

Critical steepness ratio \(\kappa_c = 1/(1-\lambda)\).

The critical value of kappa in the zELA(Lt) closure (theory paper, Coupled glacial-fluvial steady state). With the default exponent-derived lam, this is \((k+d\sigma)/k\); an explicit lam override supplies the corresponding closure scale directly.

property zfo

divide elevation of the unglaciated orogen.

Type:

Reference fluvial relief

NH(Lt: float | None = None) → float[source]

Dimensionless ice-thickness scale \(H/z_{\rm ELA} = (k_H/z_{\rm ELA})\, L_t^{\Lambda}\).

Returns:

The H/zELA scale; NaN when there is no glacier (Lt non-finite or Lt <= xo).

Return type:

float

property Lambda

Ice-thickness size exponent Lambda = gamma (1+d)/(1+phi).

(Same exponent that drives NH; named to match the AAR derivation.)

aar(surface='gsurface', *, solution=None, alpha_g=None, kh=None)[source]

Steady-state accumulation-area ratio AAR = A_c / (A_c + A_a).

Accumulation area is the full catchment above the ELA, \(A_c = k_h x_{\rm ELA}^d\); ablation area is the ice footprint below the ELA,

\[A_a = \frac{\alpha_g k_H}{k}\, L_t^{\Lambda+1}\, B(a, b) \left[1 - I_{u_{\rm ELA}^k}(a, b)\right], \qquad a = \frac{d\gamma + 1}{k},\quad b = \gamma + 1.\]

Both reduce to a single number once u_ELA is fixed by the surface closure (see AARResult).

Parameters:

  • surface ({'gsurface', 'powerlaw', 'both'}) – Which ELA-crossing closure to use. Default ‘gsurface’ (the exact in-model G(xo) crossing). ‘powerlaw’ reads the crossing from the power-law accumulation ansatz (u_ELA^k = lam); the ansatz serves only to close the accumulation integral and misplaces the crossing, so treat it as a comparison tool, not a result. ‘both’ returns a dict keyed by both labels.

  • solution (Solution, optional) – Which steady state to use. Defaults to primary.

  • alpha_g (float, optional) – Override the stored width-to-thickness ratio / Hack coefficient.

  • kh (float, optional) – Override the stored width-to-thickness ratio / Hack coefficient.

Returns:

  • AARResult, or {‘powerlaw’ (AARResult, ‘gsurface’: AARResult}.)

  • aar is NaN for a fluvial (unglaciated) state.

profile(n_points: int = 3000, *, solution: Solution | None = None)

Surface elevation profile (x, z) for the chosen solution.

For mixed solutions x is concentrated such that both the glacial and fluvial sections are well resolved.

Parameters:

  • n_points (int, optional) – Number of sample points (mixed profiles keep at least 1000 per section). Default 3000.

  • solution (Solution, optional) – Steady state to sample. Defaults to primary.

Returns:

  • x (ndarray) – Distance from the divide [m], geometrically spaced on [xo, L].

  • z (ndarray) – Surface elevation above base level [m], same shape as x. Both arrays are empty when no steady state exists.

thickness(x, *, solution: Solution | None = None)

Ice thickness profile for the chosen solution.

\[H(x) = z_{\rm ELA}\, N_H \left[\bigl(1 - (x/L_t)^k\bigr)\, (x/L_t)^d\right]^{\gamma}\]
Parameters:

  • x (array_like) – Distance from the divide [m].

  • solution (Solution, optional) – Steady state to sample. Defaults to primary.

Returns:

Width-mean ice thickness [m], same shape as x; zero outside the ice extent and identically zero for fluvial (or absent) solutions.

Return type:

ndarray

plot(ax=None, *, solution: Solution | None = None, n_points: int = 3000, show_ela: bool = True, show_thickness: bool = True, show_bedrock: bool = True, topo_color='k', ice_color='#b0d6e2', ice_alpha: float = 0.55, bedrock_color='#888888', bedrock_fade='auto', bedrock_alpha: float = 0.6, label: str | None = None)

Plot the steady-state surface and (optionally) the ice column.

Parameters:

  • ax (matplotlib.axes.Axes, optional) – Axis to draw into; a new 8x3-inch figure is created when None.

  • solution (Solution, optional) – Steady state to draw. Defaults to primary.

  • n_points (int, optional) – Sample count for the surface profile. Default 3000.

  • show_ela (bool, optional) – Draw the ELA as a dashed horizontal red line. Default True.

  • show_thickness (bool, optional) – Fill the ice column (bed to surface) and draw the sub-ice bed line. Default True.

  • show_bedrock (bool, optional) – Render the bedrock fill below the bed surface. Default True.

  • topo_color (color, optional) – Line color for the topographic surface and sub-ice bed. Default 'k'.

  • ice_color (color, optional) – Fill color of the ice column. Default '#b0d6e2'.

  • ice_alpha (float, optional) – Opacity of the ice fill. Default 0.55.

  • bedrock_color (color, optional) – Fill color of the bedrock. Default '#888888'.

  • bedrock_fade (float | 'auto' | None) – Vertical e-folding scale [m] for the bedrock alpha fade below the bedrock surface. 'auto' picks 0.25 * relief. None falls back to a flat fill from a baseline up to the surface.

  • bedrock_alpha (float) – Maximum opacity of the bedrock right at the surface. Default 0.6.

  • label (str, optional) – Legend label for the surface line.

Returns:

The axis drawn into (unchanged when there is no solution to draw).

Return type:

matplotlib.axes.Axes

Notes

The bed under ice is reconstructed on the channel-floor datum, zb = z - HC_OVER_H * H; H is the width-mean depth and is not rescaled. See docs/guides/concepts.md for the public datum convention.

siim.analytical.profiles.sweep(ks, cs, zELA, L, *, cls=<class 'siim.analytical.profiles.GeneralProfile'>, **fixed)[source]

Solve a grid of steady states (regime diagrams, paper figures).

Broadcasts ks, cs, zELA and L against each other (numpy rules) and constructs cls at every grid point; parameters held fixed across the sweep go in fixed. Returns a dict of arrays in the broadcast shape:

  • 'regime' : int8 code, -1 none / 0 fluvial / 1 mixed / 2 glacial

  • 'Lt', 'zt', 'zo' : primary-solution values (NaN where absent)

  • 'n_stable' : number of stable steady states at that point

  • 'bistable' : bool

plus 'regime_codes' mapping the names to the codes. Example:

out = sweep(ks=120.0, cs=np.linspace(50, 400, 80),
            zELA=np.linspace(0, 2500, 120)[:, None], L=1e5)
plt.pcolormesh(..., out['regime'])
class siim.analytical.profiles.MarginalCoulombProfile(ks=None, cs=None, zELA=None, L=None, *, sigma=0.5, xo=None, k=1.0, Ac=2.5e-24, alpha_g=5.0, tau_c=100000.0, lam_c=1000.0, beta=0.01, kh=5.0)[source]

Bases: object

Steady-state marginal-Coulomb profile (phi = theta = 1/2, d = 2).

The exact closed (arcosh) form of the case d*phi = 1, where the general beta-function machinery is logarithmically divergent — which makes this class the closed-form cross-check of GeneralProfile in that limit.

Parameters:

  • ks (float) – Fluvial and glacial steepness indices [m].

  • cs (float) – Fluvial and glacial steepness indices [m].

  • zELA (float) – Equilibrium-line altitude [m].

  • L (float) – Orogen half-width [m].

  • sigma (float, optional) – Hack-flux exponent for ice-flux integration. Default 0.5.

  • xo (float, optional) – Hillslope cutoff [m]. Defaults to L/1000.

  • k (float, optional) – Ice-flux profile exponent (Beff(x) = beta*(zo-zELA)*[1-(x/Lt)^k]). Default 1.0 (the linear-ansatz baseline, matching GeneralProfile).

  • Ac (float, optional) – Glen’s flow-law coefficient A [Pa^-3 s^-1]; the deformation prefactor 2A/5 is applied internally. Default constants.AC (2.5e-24).

  • alpha_g (float, optional) – Valley width-to-thickness ratio. Default constants.ALPHA_G (currently 5).

  • tau_c (float, optional) – Coulomb yield stress [Pa]. Default constants.TAU_C (1e5).

  • lam_c (float, optional) – Coulomb sliding length [m]. Default constants.LAMBDA_C (1e3).

  • beta (float, optional) – Mass-balance gradient [yr^-1]. Default constants.BETA (1e-2).

  • kh (float, optional) – Hack’s-law prefactor (A = kh * x^d). Default constants.KH (5).

solutions

All steady states found (stable + the warm saddle if it exists).

Type:

list[Solution]

primary

The cold/largest-Lt stable solution. Used by profile, thickness, plot and the Lt, zt, zo, regime convenience attributes.

Type:

Solution

bistable

True when more than one stable steady state exists.

Type:

bool

solve()[source]

Re-run the marginal-Coulomb solve from the current raw parameter attributes. Called by __init__; call it again after mutating a parameter attribute (e.g. self.zELA) to refresh the object (construct-then-mutate-then-solve(); idempotent; audit m27).

property kappa

Master control ratio \(\kappa = c_s/k_s\).

In the marginal-Coulomb case both length exponents vanish (r = 1 - d*theta = 0), so kappa is the bare steepness ratio.

property kappa_c

Critical steepness ratio \(\kappa_c = (k + d\sigma)/k = 1/(1-\lambda)\).

property lam

Flux-partition ratio \(\lambda = d\sigma/(d\sigma + k)\), the fraction of accumulation-area ice flux already committed at the ELA crossing (fixed by the exponents in this class).

property zfo

ks * log(L/xo).

Type:

Reference fluvial relief

property zgo

cs * log(2^(2/k) * L/xo).

Type:

Reference glacial-to-base-level relief

NH(Lt: float | None = None) → float[source]

Dimensionless ice-thickness scale H/zELA for given Lt.

\[N_H = \frac{L_t}{z_{\rm ELA}} \left[ \left(\frac{c_s}{\lambda_\tau}\right)^{3} + \frac{\lambda_c\, c_g\, \lambda_\tau^{2}} {k_h\,\beta\,\lambda\,G_o} \right]^{-1/3}\]

with \(c_g = k_t\,\alpha_g\,(2A_c/5)\,(\rho_i g)^3\), \(\lambda_\tau = \tau_c/(\rho_i g)\), and the marginal-Coulomb shape factor \(G_o = G(x_o; L_t) = B_o/k\) where \(B_o = 2\,\mathrm{arccosh}\bigl((L_t/x_o)^{k/2}\bigr)\) — matching the physical front end’s coulomb NH (quadrature-verified; the old phi factor in place of the 1/k was a structural error, ~phi/k off).

Returns:

The H/zELA scale; NaN when there is no glacier (or when lam = 0, i.e. no ice flux).

Return type:

float

profile(n_points: int = 3000, *, solution: Solution | None = None)

Surface elevation profile (x, z) for the chosen solution.

For mixed solutions x is concentrated such that both the glacial and fluvial sections are well resolved.

Parameters:

  • n_points (int, optional) – Number of sample points (mixed profiles keep at least 1000 per section). Default 3000.

  • solution (Solution, optional) – Steady state to sample. Defaults to primary.

Returns:

  • x (ndarray) – Distance from the divide [m], geometrically spaced on [xo, L].

  • z (ndarray) – Surface elevation above base level [m], same shape as x. Both arrays are empty when no steady state exists.

thickness(x, *, solution: Solution | None = None)

Ice thickness profile for the chosen solution.

\[H(x) = z_{\rm ELA}\, N_H \left[\bigl(1 - (x/L_t)^k\bigr)\, (x/L_t)^d\right]^{\gamma}\]
Parameters:

  • x (array_like) – Distance from the divide [m].

  • solution (Solution, optional) – Steady state to sample. Defaults to primary.

Returns:

Width-mean ice thickness [m], same shape as x; zero outside the ice extent and identically zero for fluvial (or absent) solutions.

Return type:

ndarray

plot(ax=None, *, solution: Solution | None = None, n_points: int = 3000, show_ela: bool = True, show_thickness: bool = True, show_bedrock: bool = True, topo_color='k', ice_color='#b0d6e2', ice_alpha: float = 0.55, bedrock_color='#888888', bedrock_fade='auto', bedrock_alpha: float = 0.6, label: str | None = None)

Plot the steady-state surface and (optionally) the ice column.

Parameters:

  • ax (matplotlib.axes.Axes, optional) – Axis to draw into; a new 8x3-inch figure is created when None.

  • solution (Solution, optional) – Steady state to draw. Defaults to primary.

  • n_points (int, optional) – Sample count for the surface profile. Default 3000.

  • show_ela (bool, optional) – Draw the ELA as a dashed horizontal red line. Default True.

  • show_thickness (bool, optional) – Fill the ice column (bed to surface) and draw the sub-ice bed line. Default True.

  • show_bedrock (bool, optional) – Render the bedrock fill below the bed surface. Default True.

  • topo_color (color, optional) – Line color for the topographic surface and sub-ice bed. Default 'k'.

  • ice_color (color, optional) – Fill color of the ice column. Default '#b0d6e2'.

  • ice_alpha (float, optional) – Opacity of the ice fill. Default 0.55.

  • bedrock_color (color, optional) – Fill color of the bedrock. Default '#888888'.

  • bedrock_fade (float | 'auto' | None) – Vertical e-folding scale [m] for the bedrock alpha fade below the bedrock surface. 'auto' picks 0.25 * relief. None falls back to a flat fill from a baseline up to the surface.

  • bedrock_alpha (float) – Maximum opacity of the bedrock right at the surface. Default 0.6.

  • label (str, optional) – Legend label for the surface line.

Returns:

The axis drawn into (unchanged when there is no solution to draw).

Return type:

matplotlib.axes.Axes

Notes

The bed under ice is reconstructed on the channel-floor datum, zb = z - HC_OVER_H * H; H is the width-mean depth and is not rescaled. See docs/guides/concepts.md for the public datum convention.

Physical front end — siim.analytical.steady_state

Physical-parameter front end for the analytical steady state.

SteadyStateProfile is the preferred public name; the class object is analytical_steady_state_solution (kept as the historical/primary name, with SteadyStateProfile a module-level alias). It maps the model’s physical parameters — uplift, erodibility, climate, sliding law — onto the steepness-form theory and solves the coupled glacial-fluvial steady state.

Unlike the steepness-form classes in siim.analytical.profiles, the glacial steepness index here is not a constant: cs depends on Lt through the shape factor Go (C = Co (kh beta lam Go)^mu), so the closure zELA(Lt) is assembled with cs(Lt) inside and handed to the same shared root solver (siim.analytical.core.find_closure_roots) that GeneralProfile and MarginalCoulombProfile use. One engine, three front ends.

The k = 1 linear mass-balance ansatz recovers the sliding-ice incision model of Deal and Prasicek [DP21]; the fluvial limb is the stream-power model of Whipple and Tucker [WT99].

Imports stay numpy/scipy-light: usable without the model stack.

class siim.analytical.steady_state.analytical_steady_state_solution(user_params=None)[source]

Bases: object

Standalone steady-state analytical solution for the coupled glacial-fluvial profile.

Constructed from a dict of user parameters; computes the full SS profile eagerly. Same parameter conventions as siim1d, restricted to the subset the analytical needs (no simulation-loop params). U and zELA must be scalars.

Use:

a = analytical_steady_state_solution({"zELA": 1000, "Ko": 1e-6, ...})
a.surface, a.bed   # arrays on a.x  (None if no glacier-or-fluvial regime)
a.zo, a.Lt         # divide elevation, terminus position
a.analytical_z(xp), a.analytical_zb(xp), a.analytical_ice_thickness(xp)
a.solutions        # all steady states found (incl. warm saddles)
solve()[source]

Re-run the full SS solve (idempotent). Updates self.surface/.bed.

plot(ax=None, *, show_ela=True, show_thickness=True, show_bedrock=True, topo_color='k', ice_color='#b0d6e2', ice_alpha=0.55, bedrock_color='#888888', bedrock_fade='auto', bedrock_alpha=0.6, label=None)[source]

Plot the steady-state profile in the paper’s canonical profile style — the same renderer as GeneralProfile.plot (analytical.profiles.draw_profile), so a basin-fitted analytical reference overlays in the identical visual language as the standalone profiles. bed is the channel-floor datum (z - HC_OVER_H*H).

Plotted against distance from the divide, whichever end holds base level. Returns the axis (a new 8x3-inch figure when ax is None); a no-op axis when there is no steady-state glacier to draw.

plot_profile(fields='elevation', *, axes=None)[source]

Stacked analytical-profile panels (filled-band notebook style), one per field. Mirrors the five fields of siim1d’s plot.profile – ‘elevation’, ‘ice_thickness’, ‘erosion_rate’, ‘ice_flux’, ‘water_flux’ – but draws the analytical steady state only (no extracted 2D channel).

Unlike plot (the canonical single-axis draw_profile renderer), this reproduces the two-panel bed/ice/thickness layout. bed is the channel-floor datum (z - HC_OVER_H*H); at steady state the erosion rate is uniformly U. Plotted against x (distance, km).

Parameters:

  • fields (str or sequence of str) – One field name, or an ordered list -> one stacked panel each.

  • axes (matplotlib.axes.Axes or sequence, optional) – Draw into existing axes (length must match fields); otherwise a new figure is created.

Return type:

(fig, axes), or (None, None) when there is no steady-state glacier.

analytical_z_fluvial(xp)[source]

Fluvial profile (x > Lt branch).

\[z(x) = k_s\, L^{1-d\theta}\, B\!\left(1 - x/L;\ 1,\ 1-d\theta\right)\]

with xp = x/L the dimensionless along-flow coordinate. xp is clamped to >= xo/Ld (the channel-head cutoff) so the q=1 log singularity at xp=0 doesn’t appear and x < xo reads as ‘at the channel head’.

Returns:

Surface elevation above base level [m], shape of xp.

Return type:

ndarray or float

analytical_z_glacial(xp)[source]

Glacial profile (x <= Lt branch).

\[z(x) - z(L_t) = c_s\, L_t^{r}\, G(x, L_t)\]

with xp = x/Ld the dimensionless along-domain coordinate. xp is clamped to >= xo/Ld so x < xo reads as ‘at the channel head’; the clamp guarantees upper = 1 - (x/Lt)^k <= 1 - (xo/Lt)^k, so no separate physical cap is needed inside _G_at.

Returns:

Surface elevation above base level [m], shape of xp.

Return type:

ndarray or float

analytical_ice_thickness(xp)[source]

Ice thickness profile (dimensional, despite the internal H’).

\[H(x) = z_{\rm ELA}\, N_H \left[\bigl(1-(x/L_t)^k\bigr)\, (x/L_t)^d\right]^{\gamma}\]

Dispatches on self.sliding_law:

  • coulomb: gamma = phi with the coulomb NH

  • eff-exp / power: gamma = (2/3)(1/3 + phi) with the eff-exp NH

The Lt-exponent is gamma*(1+d)/(1+phi) for both. Pre-factors (cs, Lt, Go, lam, r) carry the law-specific phi via the per-law dispatch.

Returns:

Width-mean ice thickness [m] (the dimensionless H’ = H/zELA is multiplied back by zELA before returning), shape of xp. Zero outside the ice extent (x > L_t) and identically zero for a fluvial or absent solution (L_t is NaN) — matching the profiles._thickness family (audit m26); no NaN/RuntimeWarning past the terminus.

Return type:

ndarray or float

analytical_z(xp)[source]

Surface elevation [m] at xp = x/Ld: glacial branch where xp <= Lt/Ld, fluvial elsewhere. Scalar in, scalar out.

analytical_zb(xp)[source]

Bed elevation [m] at xp = x/Ld on the channel-floor datum: zb = z - HC_OVER_H * H under ice, zb = z elsewhere.

siim.analytical.steady_state.SteadyStateProfile

alias of analytical_steady_state_solution

Regime map — siim.analytical.regime

Nondimensional regime map for the coupled glacial-fluvial steady state.

The steady state of the coupled profile is controlled by two dimensionless numbers: the fluvial and glacial erosion numbers \(N_f\) and \(N_g\). This module works in the derived coordinates used for regime diagrams,

\[\kappa = N_g/N_f, \qquad Y = z_{\rm ELA}/z_{fo},\]

with all lengths in units of the orogen length L (so xo here is \(x_o/L\) and Lt is \(L_t/L\)), and all elevations in units of the reference fluvial relief \(z_{fo} = k_s L^{1-d\theta} F(x_o/L)\).

RegimeMap evaluates, over arrays of \((\kappa, Y)\), the same continuous closure that siim.analytical.profiles.GeneralProfile solves per point in physical units:

\[\begin{split}Y(L_t';\,\kappa)\,F(x_o') = \begin{cases} F(L_t') + (\kappa/\kappa_c)\, L_t'^{\,r}\, G(x_o'/L_t') & L_t' \le 1 \text{ (mixed)} \\[2pt] \kappa\, L_t'^{\,r}\left[(1/\kappa_c)\, G(x_o'/L_t') - G(1/L_t')\right] & L_t' > 1 \text{ (glacial)} \end{cases}\end{split}\]

with \(r = (1-d\phi)/(1+\phi)\) and \(\kappa_c = 1/(1-\lambda)\). For the default exponent-derived \(\lambda\), \(\kappa_c = (k+d\sigma)/k\). The two branches join continuously at \(L_t' = 1\). The shape functions are evaluated through the analytically continued incomplete beta of siim.analytical.core, so the marginal-Coulomb case \(d\phi = 1\) (and \(d\theta = 1\)) needs no special casing — the kernel’s exact b = 0 logarithmic branch reproduces the \(\operatorname{arcosh}\) closed forms.

Because the closure is linear in \(\kappa\) at fixed \(L_t'\), the shape functions are evaluated once on a shared \(L_t'\) grid; per- \(\kappa\) closure traces are then cheap vector operations. Array solves locate roots by monotone-segment interpolation on that trace (accuracy set by n_trace); size-1 inputs are solved trace-free by damped Newton from a closed-form initial guess (machine precision at the root, and cheap enough that consumers may construct a fresh map per evaluation).

Everything here is numpy/scipy only (no model stack, no matplotlib).

class siim.analytical.regime.RegimeMap(*, d=2.0, sigma=0.5, phi=0.5, theta=0.5, k=1.0, xo=0.001, lam=None, Lt_max=50.0, n_trace=4096)[source]

Bases: object

Steady-state regime map in \((\kappa, Y)\) coordinates.

Parameters:

  • d (float, optional) – Hack-area exponent. Default 2.0.

  • sigma (float, optional) – Hack-flux exponent for ice-flux integration; enters via \(\lambda = d\sigma/(d\sigma+k)\). Default 0.5.

  • phi (float, optional) – Glacial concavity index, in [0, 1). Default 0.5 (Coulomb).

  • theta (float, optional) – Fluvial concavity index. Default 0.5.

  • k (float, optional) – Ice-accumulation shape exponent. Default 1.0.

  • xo (float, optional) – Channel-head cutoff \(x_o/L\), in (0, 1). Default 1e-3.

  • lam (float, optional) – AAR-like ratio override. Defaults to \(d\sigma/(d\sigma+k)\).

  • Lt_max (float, optional) – Far end of the fully-glaciated closure trace, in units of L. Default 50.

  • n_trace (int, optional) – Shared-grid resolution of the closure trace, built lazily on the first array solve (size-1 solves are trace-free Newton). Default 4096.

kappa_c

Critical steepness ratio \(1/(1-\lambda)\). With the default exponent-derived \(\lambda\), this is \((k + d\sigma)/k\).

Type:

float

alpha

Geometric cutoff factor \(G(x_o')/F(x_o')\) — the slope of the cold (fully-glaciated) regime boundary \(Y = \alpha\kappa/\kappa_c\).

Type:

float

r

Glacial length exponent \((1-d\phi)/(1+\phi)\).

Type:

float

lam

AAR-like elevation-interval ratio \(\lambda\).

Type:

float

Notes

The array methods (Lt(), zo(), zt(), masks()) return the cold branch: the largest-\(L_t'\) root of the continuous closure — the most-glaciated stable state, the one a cooling history lands on. branch='warm' returns the smallest root instead (the small-glacier branch followed when warming out of a bistable state). Where no glaciated root exists the state is fluvial: Lt is NaN and zo is 1.

F(u)[source]

Dimensionless fluvial profile shape \(F(u)\), u = x/L.

\[F(u) = \mathcal{B}(1-u;\,1,\,1-d\theta) = \frac{1 - u^{1-d\theta}}{1-d\theta} \;\xrightarrow{\;d\theta\to 1\;}\; \ln(1/u)\]

Evaluated through expm1 so the \(d\theta = 1\) logarithmic limit is reached continuously, with no branch tolerance. Vectorizes over u.

G(u)[source]

Dimensionless glacial profile shape \(G(u)\), u = x/Lt.

\[G(u) = \tfrac{1}{k}\, \mathcal{B}\!\left(1-u^k;\,1-\phi,\,\tfrac{1-d\phi}{k}\right)\]

evaluated through the eps-form of the analytically continued incomplete beta (exact for u << 1 and at the marginal case \(d\phi = 1\), where b = 0). Vectorizes over u.

closure_Y(Lt, kappa)[source]

The continuous closure \(Y(L_t'; \kappa)\).

The climate \(Y = z_{\rm ELA}/z_{fo}\) for which Lt is a steady-state terminus at steepness ratio kappa. This is the dimensionless form of siim.analytical.profiles.GeneralProfile._zELA_of_Lt(). Broadcasts over Lt and kappa.

Lt(kappa, Y, *, branch='cold')[source]

Steady-state terminus position \(L_t'/L\) over (kappa, Y).

NaN where no glaciated steady state exists (fluvial), or where the cold root lies beyond Lt_max. Size-1 inputs are solved to machine precision; arrays by trace interpolation (accuracy set by n_trace).

zt(kappa, Y, *, branch='cold', Lt=None)[source]

Terminus elevation \(z_t/z_{fo}\).

\(F(L_t')/F(x_o')\) on the mixed branch; \(-(\kappa/F(x_o'))\,L_t'^{\,r}\,G(1/L_t')\) (below base level) on the glacial branch; NaN where fluvial. Pass Lt to reuse an already-computed terminus array.

zo(kappa, Y, *, branch='cold', Lt=None)[source]

Divide elevation \(z_o/z_{fo}\).

\(z_o = z_t + (\kappa/F(x_o'))\,L_t'^{\,r}\,G(x_o'/L_t')\) on a glaciated root, exactly as _make_solution() builds it; 1 where the state is fluvial.

masks(kappa, Y, *, branch='cold')[source]

(fluvial, mixed, glacial) boolean masks from the cold-branch root.

mixed where the root sits in (xo, 1]; glacial where it exceeds 1 (fully glaciated to base level); fluvial where no glaciated root exists. Note the bistable strip (kappa < kappa_c, Y > 1 below the saddle-node) reads as mixed under this convention — the diagram shows the glaciated attractor there; overlay saddle_node() to mark the strip.

closure_extrema(kappa)[source]

Interior (Y_max, Y_min) of the mixed-branch closure trace.

The fold pair of the partial-glacial closure on (xo, 1): an interior maximum (fluvial/mixed bistability, the warm saddle) and/or an interior minimum (small/large-glacier bistability, the cooling spinodal). Either is NaN when the trace is monotone there.

saddle_node(n_samples=400)[source]

Parametric saddle-node (L1) boundary of the mixed closure.

Setting \(\partial Y/\partial L_t' = 0\) on the mixed branch gives, with \(u = x_o'/L_t'\) and \(D(u) = r\,G(u) + u^{1-d\phi}(1-u^k)^{-\phi}\),

\[\frac{\kappa}{\kappa_c} = \frac{L_t'^{\,1-r-d\theta}}{D(u)},\]

and the fold’s \((Y, z_o/z_{fo})\) follow from the closure and the profile geometry at that \((L_t', \kappa)\). The trace is sampled with log concentration toward both endpoints and prepended with the exact corner \((0, 1, 1)\) at \(L_t' = x_o'\); in the marginal-Coulomb case it reduces to the \(\kappa/\kappa_c = \sqrt{1 - v^{-2}}\), \(v = (L_t'/x_o')^{k/2}\) arcosh trace.

Returns:

Arrays (kappa, Y, zo, Lt) along the fold, ordered by Lt.

Return type:

SaddleNodeBoundary

class siim.analytical.regime.SaddleNodeBoundary(kappa: ndarray, Y: ndarray, zo: ndarray, Lt: ndarray)[source]

Bases: object

Parametric saddle-node (fold) trace of the partial-glacial closure.

The locus in \((\kappa, Y)\) where the warm and cold mixed roots merge — the warm boundary of the bistable strip on a regime diagram. Sampled parametrically in the terminus position Lt (in units of L); the first entry is the exact \(L_t' \to x_o'\) corner \((\kappa, Y, z_o/z_{fo}) = (0, 1, 1)\).