Public API

A superstruct needed for the forward integration of the model.

Fields

• domain: the AbstractDomain defining the spatial discretization.
• c: the PhysicalConstants defining the physical constants of the model.
• bcs: the BoundaryConditions defining the boundary conditions of the model.
• sealevel: the RegionalSeaLevel defining the sea level evolution.
• solidearth: the SolidEarth defining the solid earth properties.
• opts: the SolverOptions controlling the solver options.
• tools: the GIATools providing tools for GIA computations.
• ref: the ReferenceState defining the reference state of the model.
• now: the CurrentState defining the current state of the model.
• ncout: the NetcdfOutput controlling the NetCDF output.
• nout: the NativeOutput controlling the native output.
• timer: the Timer controlling and recording timing information.
• simobs: a vector of SimulatedObservable defining simulated observables to be computed during integration.

Control options relative to solving a Simulation.

Fields

• integ: the AbstractIntegrator used to integrate the ODE forward in time, one of BS3Integrator (adaptive, default), Tsit5Integrator (adaptive), RKCIntegrator (adaptive, stabilised for stiff problems) or EulerIntegrator (fixed step). Each integrator carries its own settings — tolerances, step-size bounds — as fields of its own struct.
• dt_sparse_diagnostics: the time interval between updates of the diagnostics variables (elastic displacement, sea-surface elevation, etc.).
• show_progress: whether to report the simulation progress. When true, run! displays a live progress bar (ForwardProgress).
• dt_walltime: minimum wall time in seconds between two refreshes of that progress bar. Refreshing reduces over the whole grid, so this bounds the reporting cost by wall time rather than by step count.
• fft: the AbstractFFTBackend used for the spectral step. A purely numerical choice, orthogonal to the mantle rheology.
• transition: the AbstractTransition used to smooth the transition between grounded and floating ice, and between ocean and land.

integ is a concrete type parameter rather than an abstract field, so sim.opts.integ infers to the integrator's own type. Code that branches on the integrator — the AD extensions in particular — then resolves that branch at compile time instead of paying to compile every arm of it.

run!(sim::Simulation)

Solve the isostatic adjustment problem defined in sim::Simulation, integrating it forward over sim.timer.t_span with the integrator in sim.opts.integ::AbstractIntegrator and writing output at the requested times.

step!(integ::FastIsostasy.AbstractIntegratorState, Δt)
step!(
    integ::FastIsostasy.AbstractIntegratorState,
    Δt,
    force_dt::Bool
)

Advance integrator by the interval Δt, using internal adaptive substeps and writing any output that falls within the interval. force_dt is accepted for backward compatibility; the step always stops exactly at t + Δt.

init_integrator(f!, u0, tspan, alg::AbstractIntegrator, p = nothing)

Build (and prime) the integrator state for the in-place ODE f!(du, u, p, t), ready to be advanced by solve_to! / step!. Returns a TableauIntegratorState, or an RKCIntegratorState when alg is an RKCIntegrator.

Tolerances and step-size bounds are read off alg — set them there, e.g. Tsit5Integrator(reltol = 1e-8) or EulerIntegrator(dt = 100.0).

See init_integrator(sim::Simulation) for the Simulation-level entry point, which takes alg from sim.opts.integ.

init_integrator(
    sim::Simulation
) -> Union{FastIsostasy.RKCIntegratorState{_A, _B, typeof(update_diagnostics!), P} where {_A, _B, P<:Simulation}, FastIsostasy.TableauIntegratorState}

Initialise the integrator of sim::Simulation, which can subsequently be advanced manually with step!(integrator, Δt, force_dt) (e.g. when coupling to an external ice-sheet model).

Carry physical constants, with default values that can be changed by the user. For instance:

c = PhysicalConstants(rho_ice = 0.93)   # (kg/m^3)

All constants are given in SI units (kilogram, meter, second).

Fields

• mE: Earth's mass (kg)
• r_equator: Earth radius at equator (m)
• r_pole: Earth radius at pole (m)
• A_ocean_pd: Ocean surface (m^2) as in Goelzer (2020) before Eq. (9)
• g: Mean Earth acceleration at surface (m/s^2)
• G: Gravity constant (m^3 kg^-1 s^-2)
• seconds_per_year: (s)
• rho_ice: (kg/m^3)
• rho_water: (kg/m^3)
• rho_seawater: (kg/m^3)
• rho_sw_ice: (dimensionless)

Define a self-contained, dependency-free time integrator for in-place ODEs.

Available subtypes:

• EulerIntegrator
• BS3Integrator
• Tsit5Integrator
• RKCIntegrator

EulerIntegrator(dt)
EulerIntegrator(; dt)

Fixed-step explicit Euler. Non-adaptive, so dt is the step size and is required — there is no error estimate to adapt it with.

Fields

• dt::AbstractFloat: the fixed step size

BS3Integrator(; kwargs...)
BS3Integrator{T}(; kwargs...)

Bogacki-Shampine 3(2) embedded pair (FSAL, adaptive). Third-order accurate solution with a second-order embedded error estimate.

Fields

• reltol::AbstractFloat: relative error tolerance of the adaptive controller

• abstol::AbstractFloat: absolute error tolerance of the adaptive controller

• dt_min::AbstractFloat: lower bound on the adaptive step size

• dt_max::AbstractFloat: upper bound on the adaptive step size

• dt0::AbstractFloat: initial step size; the controller grows it from here, at most 10x per step

Tsit5Integrator(; kwargs...)
Tsit5Integrator{T}(; kwargs...)

Tsitouras 5(4) embedded pair (FSAL, adaptive). Fifth-order accurate solution with a fourth-order embedded error estimate.

Fields

• reltol::AbstractFloat: relative error tolerance of the adaptive controller

• abstol::AbstractFloat: absolute error tolerance of the adaptive controller

• dt_min::AbstractFloat: lower bound on the adaptive step size

• dt_max::AbstractFloat: upper bound on the adaptive step size

• dt0::AbstractFloat: initial step size; the controller grows it from here, at most 10x per step

RKCIntegrator(; kwargs...)
RKCIntegrator{T}(; kwargs...)

Stabilised explicit Runge-Kutta-Chebyshev method (RKC2, Sommeijer-Shampine- Verwer 1997), second order, damped. Real-axis stability interval grows with the square of the stage count. For stiff systems with ~1.75x speedup.

Fields

• damping::AbstractFloat: SSV damping parameter ε; must be > 0

• safety::AbstractFloat: safety factor for dt * λ_max when choosing the stage count

• smax::Int64: hard cap on the stage count per step

• reestimate_every::Int64: re-run the power iteration every N accepted steps (0 = never)

• reltol::AbstractFloat: relative error tolerance of the adaptive controller

• abstol::AbstractFloat: absolute error tolerance of the adaptive controller

• dt_min::AbstractFloat: lower bound on the adaptive step size

• dt_max::AbstractFloat: upper bound on the adaptive step size

• dt0::AbstractFloat: initial step size; the controller grows it from here, at most 10x per step

integrate(f!, u0, tspan, alg::AbstractIntegrator, p = nothing; kwargs...) -> (ts, us)

Integrate the in-place ODE f!(du, u, p, t) from tspan[1] to tspan[2].

Returns the sorted save times ts and a vector us of solution snapshots (copies) at those times.

Error tolerances and step-size bounds are fields of alg — e.g. integrate(f!, u0, tspan, Tsit5Integrator(reltol = 1e-8)) or integrate(f!, u0, tspan, EulerIntegrator(dt = 0.05)).

Keyword arguments:

• saveat: times at which to store the solution (defaults to tspan[2]). Must lie within tspan. The endpoint is always included.
• maxiters: safety cap on the number of steps per save interval.

AbstractTransition

Supertype of the transition traits, which decide how the model's non-smooth switches — grounded/ocean masks, the max(0, ·) clamps in the water column — are evaluated. Stored on SolverOptions.

SharpTransition (the default) keeps the exact max/Heaviside behaviour at zero cost. SmoothTransition replaces them by ε-smoothed counterparts, which is what makes the model differentiable across the grounding line; the masks then become eltype-T arrays instead of Bool.

SharpTransition()

Exact, non-differentiable grounding-line switches (max, min, >). Default; produces Boolean masks and leaves the forward model unchanged.

SmoothTransition(ε)

Differentiable grounding-line switches with transition length scale ε (in metres of height-above-flotation). Uses the branch-free approximations

max(x, 0) ≈ (x + √(x² + ε²)) / 2
min(x, 0) ≈ (x − √(x² + ε²)) / 2
(x > 0)   ≈ (1 + x / √(x² + ε²)) / 2

which are symmetric, GPU-friendly and converge to the exact operators as ε → 0 (error confined to |x| ≲ ε). Masks become floating-point fields in [0, 1].

Abstract type for domain representation in the model.

Available subtypes:

• RegionalDomain
• GlobalDomain: Not implemented yet, but hopefully in v3.0!

Define a regional domain, including its geometry and the architecture it should be running on.

Initialization

For a square domain with half-width W and 2^n grid points in each dimension:

domain = RegionalDomain(W, n)

For a rectangular domain with half-widths Wx and Wy and nx and ny grid points in each dimension:

domain = RegionalDomain(Wx, Wy, nx, ny)

For a rectangular domain with spanning vectors x and y:

domain = RegionalDomain(x, y)           # rectangular domain: spanning vectors x, y

Hardware

The backend keyword selects where the arrays live and where the kernels run. It takes a KernelAbstractions.jl backend, so FastIsostasy is not tied to any single GPU vendor:

using CUDA
domain = RegionalDomain(3000e3, 7, backend = CUDABackend())

Everything downstream — SolidEarth, BoundaryConditions, Simulation — picks the backend up from the domain, so this is the only place hardware is named.

Not implemented yet!

Version 3.0 will allow for global domains.

Define the boundary conditions of the problem.

Fields

• ice_thickness: an instance of AbstractIceThickness that defines how the ice thickness is updated.
• viscous_displacement: a boundary condition for the viscous displacement, defined as an OffsetBC.
• elastic_displacement: a boundary condition for the elastic displacement, defined as an OffsetBC.
• sea_surface_perturbation: a boundary condition for the sea surface perturbation, defined as an OffsetBC.

apply_bc!(H, t, it::TimeInterpolatedIceThickness)

Update the ice thickness H at time t using the method defined in it.

apply_bc!(X, bc::OffsetBC)

Apply the boundary condition bc to the matrix X in-place.

Determine how the ice thickness is updated in the model by implementing update_ice! for different subtypes.

Available subtypes

• TimeInterpolatedIceThickness
• ExternallyUpdatedIceThickness

Update the ice thickness based on a time interpolation.

Fields

• t_vec: a vector of time points at which the ice thickness is defined.
• H_vec: a vector of ice thickness values corresponding to t_vec.
• H_itp: a function that interpolates the ice thickness based on time.

Update the ice thickness externally, without any internal update.

An abstract type representing the space in which boundary conditions are defined. This typically needs to be defined when initializing an AbstractBCSpace.

Available subtypes

• RegularBCSpace
• ExtendedBCSpace

Singleton struct to impose boundary conditions at the edges of the computation domain.

Singleton struct to impose boundary conditions at the edges of the extended computation domain, which naturally arises from convolutions.

An abstract type representing a boundary condition in the context of a computational domain.

Available subtypes

• OffsetBC
• NoBC

Apply an offset to the values at the boundaries of a computational domain.

Fields

• space: the AbstractBCSpace in which the boundary condition is defined.
• x_border: the offset value to be applied at the boundaries.
• W: a weight matrix to apply the boundary condition according to some AbstractBC.

A singleton struct representing the absence of a boundary condition.

Impose a Dirichlet-like boundary condition at the corners of the computational domain.

Impose a Dirichlet-like boundary condition at the borders of the computational domain.

Impose a Dirichlet-like boundary condition at the borders of the computational domain, weighted by the distance from the center of the domain.

Impose a mean value for the field.

A struct that gathers the modelling choices for the sea-level component of the simulation. It contains:

• surface: an instance of AbstractSeaSurface to represent the sea surface.
• load: an instance of AbstractSealevelLoad to represent the sea-level load.
• bsl: an instance of AbstractBSL to represent the barystatic sea level.
• update_bsl: an instance of AbstractBSLUpdate to represent the update mechanism for the barystatic sea level.

Abstract type to compute the evolution of the barystatic sea level.

Available subtypes

• ConstantBSL
• ConstantOceanSurfaceBSL
• PiecewiseConstantBSL
• PiecewiseLinearOceanSurfaceBSL (requires using NLsolve)
• ImposedBSL
• CombinedBSL

All subtypes implement the update_bsl! function:

T, delta_V, t = Float64, 1.0e9, 0.0             # Example values
bsl = PiecewiseConstantBSL()        # or any other subtype!
update_bsl!(bsl, delta_V, t)

Define a reference of barystatic sea level and ocean surface area. Used in all subtypes of AbstractBSL to compute the BSL evolution.

Fields

• z: the reference BSL (m), which defaults to 0 (reference year 2020).
• A: the reference ocean surface area (m^2) computed based on z.
• z_vec: a vector of BSL values (m) used for interpolation.
• A_vec: a vector of ocean surface area values (m^2) used for interpolation.
• A_itp: an interpolator function for ocean surface area over BSL.

In the constructor, T determines the floating point arithmetic used in all computations, and itp_kwargs allows customization of the interpolation.

Example usage:

ref = ReferenceBSL()              # assume BSL = 0

Custom:

ref = ReferenceBSL(z = 0.1)       # assume BSL = 0.1 m and compute A accordingly

A mutable struct containing:

• ref: an instance of ReferenceBSL.
• z: the BSL, considered constant in time.
• A: the ocean surface area, considered constant in time.

Assume that the BSL is constant in time.

A mutable struct containing:

• ref: an instance of ReferenceBSL.
• z: the BSL at current time step.
• A: the ocean surface area, considered constant in time.

Assume that the ocean surface is constant in time and that the BSL evolves only according to the changes in ice volume covered by the RegionalDomain.

A mutable struct containing:

• ref: an instance of ReferenceBSL.
• z: the BSL at current time step.
• A: the ocean surface at current time step.

Assume that the ocean surface evolves in time according to a piecewise constant function of the BSL, which evolves in time according to the changes in ice volume covered by the RegionalDomain.

A mutable struct containing:

• ref: an instance of ReferenceBSL.
• z: the BSL at current time step.
• t_vec: the time vector.
• z_vec: the BSL values corresponding to the time vector.
• z_itp: an interpolation of z_vec over t_vec.

Impose an externally computed BSL, which is internally computed via a time interpolation.

A mutable structcontaining:

• bsl1: an ImposedBSL.
• bsl2: an AbstractBSL.

This imposes a mixture of BSL. For instance, if you simulate Antarctica over the LGM, you can impose an offline BSL contribution from the other ice sheets via bsl1. The contribution of Antarctica will be intercatively added to this via bsl2.

An abstract type to determine how the barystatic sea level (BSL) is updated over time.

Available subtypes

• InternalBSLUpdate
• ExternalBSLUpdate

Update BSL internally by the model based on the change in ice volume.

Update the BSL is externally, without any internal update.

update_bsl!(bsl::ConstantBSL, delta_V, t)

Update the BSL and ocean surface based on the input delta_V (in m^3) and on a subtype of AbstractBSL.

Abstract type for sea surface representation. Available subtypes are:

• LaterallyConstantSeaSurface
• LaterallyVariableSeaSurface

Assume a laterally constant sea surface across the domain. This means that the gravitatiional response is ignored.

Assume a laterally variable sea surface across the domain. This means that the gravitational response is included in the sea surface perturbation.

update_dz_ss!(
    sim::Simulation,
    sl::LaterallyVariableSeaSurface
)

Update the SSH perturbation dz_ss by convoluting the Green's function with the load anom.

Abstract type for sea-level load representation. Available subtypes are:

• NoSealevelLoad
• InteractiveSealevelLoad

Assume no sea-level load, i.e. the bedrock deformation is not influenced by the sea-level.

Assume an interactive sea-level load, i.e. the bedrock deformation is influenced by the sea-level.

columnanom_water!(
    sim::Simulation,
    ol::InteractiveSealevelLoad
)

Update the water column based on the instance of AbstractSealevelLoad.

Return a struct containing all information related to the lateral variability of solid-Earth parameters. To initialize with values other than default, run:

domain = RegionalDomain(3000e3, 7)
lb = [100e3, 300e3]
lv = [1e19, 1e21]
solidearth = SolidEarth(domain, layer_boundaries = lb, layer_viscosities = lv)

which initializes a lithosphere of thickness $T_1 = 100 \mathrm{km}$, a viscous channel between $T_1$and $T_2 = 300 \mathrm{km}$and a viscous halfspace starting at $T_2$. This represents a homogenous case. For heterogeneous ones, simply make lb::Vector{Matrix}, lv::Vector{Matrix} such that the vector elements represent the lateral variability of each layer on the grid of domain::RegionalDomain.

Available subtypes are:

• RigidLithosphere
• LaterallyConstantLithosphere
• LaterallyVariableLithosphere

Assume a rigid lithosphere, i.e. the elastic deformation is neglected.

Assume a laterally constant lithospheric thickness (and rigidity) across the domain. This generally improves the performance of the solver, but is less realistic.

Assume a laterally variable lithospheric thickness (and rigidity) across the domain. This generally improves the realism of the model, but is more computationally expensive.

update_elasticresponse!(
    sim::Simulation,
    lithosphere::AbstractLithosphere
)

Update the elastic response by convoluting the Green's function with the load anom. To use coefficients differing from ^{[Farrell1972]}, see GIATools.

The rheology of the mantle. This axis describes what is modelled; how the spectral step is computed is the orthogonal AbstractFFTBackend axis.

Note that none of these subtypes carries the unrelaxed elastic response: that is computed separately by the Farrell Green's-function convolution, which dispatches on AbstractLithosphere.

Available subtypes are:

• RigidMantle
• RelaxedMantle
• ViscousMantle
• TransientCreepMantle

Assume a rigid mantle that does not deform.

Assume a relaxed mantle that deforms according to a relaxation time. This is generally less realistic and offers worse performance than a viscous mantle. It is only included for legacy purpose (e.g. comparison among solvers).

Assume a viscous mantle that deforms according to a viscosity, i.e. steady-state (secondary) creep only.

The name is deliberately not MaxwellMantle: a Maxwell body is a spring in series with a dashpot, but the spring — the unrelaxed elastic response — is not handled here. It is computed independently by the Farrell Green's-function convolution in update_elasticresponse!, which dispatches on AbstractLithosphere and can be switched off entirely (RigidLithosphere). This type supplies the viscous half-space only.

Whether the spectral step uses complex or half-spectrum FFTs is a separate, orthogonal choice — see AbstractFFTBackend and the fft field of SolverOptions.

This is the most realistic mantle model with a steady creep law and generally offers the best performance. It is the default mantle model used in the solver.

TransientCreepMantle(; shearmodulus, relaxation_strength, kelvin_time)

Transient (primary) creep on top of the steady creep of ViscousMantle: the steady Maxwell dashpot (viscosity η₁, taken from SolidEarth) in series with N Kelvin-Voigt branches forming a Prony series, so that the total viscous displacement splits as u = u_M + Σⱼ u_K[j]. N = 1 is the classic Burgers body; N ≈ 3-5 log-spaced branches approximate the extended Burgers / Faul-Jackson relaxation spectrum of Ivins & Caron (2021).

Named for the phenomenon rather than for one rheological body, because the N-branch Prony form is the general case and Burgers is only its N = 1 member.

Note that, unlike ViscousMantle, this model does carry shear moduli: each branch has μ₂ⱼ = shearmodulus / Δⱼ and η₂ⱼ = τⱼ μ₂ⱼ. These are internal to the transient band and remain distinct from the unrelaxed elastic response, which the Farrell convolution keeps handling as for every other mantle. It extends ViscousMantle rather than replacing it: as Δⱼ → 0 the Kelvin branches lock (u_K → 0) and the model reduces to it exactly.

Scalar (laterally constant) parameters only, and supported solely on the semi-implicit path — LaterallyConstantLithosphere or RigidLithosphere with ComplexFFTBackend and a fixed step (SolverOptions(integ = EulerIntegrator(dt = ...))).

Fields

• shearmodulus::AbstractFloat: unrelaxed shear modulus μ₁ of the mantle [Pa]

• relaxation_strength::NTuple{N, T} where {T<:AbstractFloat, N}: relaxation strength Δⱼ = μ₁/μ₂ⱼ of each Kelvin branch

• kelvin_time::NTuple{N, T} where {T<:AbstractFloat, N}: retardation time τⱼ = η₂ⱼ/μ₂ⱼ of each Kelvin branch [yr]

Example

Classic Burgers body with the Ivins & Caron (2021) Fig. 8 parameters:

julia> using FastIsostasy

julia> m = TransientCreepMantle(shearmodulus = 67e9, relaxation_strength = 1.2,
           kelvin_time = 7.14);

julia> FastIsostasy.nbranches(m)
1

See roadmaps/burgers.md for the design and derivation.

BurgersMantle(
;
    shearmodulus,
    relaxation_strength,
    kelvin_time
)

Classic Burgers body: the N = 1 member of TransientCreepMantle, a Maxwell dashpot in series with a single Kelvin-Voigt element. Thin naming convenience over TransientCreepMantle(; shearmodulus, relaxation_strength, kelvin_time) for when (Δ, τ) are already known (e.g. from a published fit), as opposed to ExtendedBurgersMantle, which fits them from a continuous relaxation-time spectrum.

ExtendedBurgersMantle(
;
    shearmodulus,
    relaxation_strength,
    alpha,
    tau_L,
    tau_H,
    nbranches
)

Extended Burgers mantle: a TransientCreepMantle with nbranches Kelvin branches fit, via fit_prony_series, to the continuous Faul-Jackson absorption-band spectrum of the extended Burgers model (EBM) of Ivins & Caron (2021). See fit_prony_series for the fitting procedure and its fit_error; this constructor discards fit_error — call fit_prony_series directly first to check it before committing to nbranches.

update_dudt!(dudt, u, sim, t, earth::SolidEarth) -> Any

Update the time derivative of the viscous displacement based on a dispatch is along three mainly orthogonal axes: the mantle rheology, the lithosphere, and the FFT backend.

Main supported combinations are:

• RigidMantle: no deformation, dudt is zero.
• RelaxedMantle with LaterallyConstantLithosphere: uses ELRA Le Meur and Huybrechts (1996) to compute the viscous response. This also works with laterally-variable relaxation time, as proposed by Coulon et al. (2021) and by van Calcar et al. (2026).
• RelaxedMantle with LaterallyVariableLithosphere: not implemented. This corresponds to what is described by Coulon et al. (2021) but is not yet implemented.
• ViscousMantle with LaterallyConstantLithosphere or RigidLithosphere: not implemented. This corresponds to what is described by Bueler et al. (2007) but is not yet implemented.
• ViscousMantle with LaterallyVariableLithosphere: This corresponds to the approach of Swierczek-Jereczek et al. (2024).

Which FFT plans the spectral solver builds. This is a numerical choice and is deliberately orthogonal to the rheology: it selects how a transform is computed, never what is being modelled. Set it through the fft field of SolverOptions.

Available subtypes:

• ComplexFFTBackend
• RealFFTBackend

Complex-to-complex plans (plan_fft / plan_ifft) over the full (nx, ny) spectrum. The default, and the FFT backend to use for production runs.

Real-to-complex plans (plan_rfft / plan_irfft). The frequency-domain arrays are (nx÷2+1, ny) rather than (nx, ny), roughly halving the memory and arithmetic cost of the spectral step.

Abstract type for layering models. Subtypes should implement the layer_boundaries function. Available subtypes are:

• UniformLayering
• ParallelLayering
• EqualizedLayering
• FoldedLayering

Struct to enforce uniform layering when passed to get_layer_boundaries. Contains:

• n_layers: the number of layers in the model.
• boundaries: the layer boundaries, which are constant across the domain.

Struct to enforce parallel layering when passed to get_layer_boundaries. Contains:

• n_layers: the number of layers in the model.
• thickness: the thickness of each layer.
• tol: a tolerance value to add to the layer boundaries.

Struct to enforce equalized layering when passed to get_layer_boundaries. Contains:

• n_layers: the number of layers in the model.
• boundaries: the layer boundaries.
• tol: a tolerance value to add to the layer boundaries.

Struct to enforce folded layering when passed to get_layer_boundaries. Contains:

• n_layers: the number of layers in the model.
• max_depth: the maximum depth of the layers.
• tol: a tolerance value to add to the layer boundaries.

get_layer_boundaries(
    domain::RegionalDomain,
    litho_thickness,
    layering
) -> Any

Compute the layer boundaries for a given AbstractLayering. Output is typically passed to SolidEarth to create a layered Earth model.

apply_calibration!(eta, calibraton::NoCalibration) -> Any

Apply the calibration to the viscosity eta.

get_effective_viscosity_and_scaling(
    domain,
    layer_viscosities,
    layer_boundaries,
    maskactive,
    lumping::TimeDomainViscosityLumping
) -> Tuple{Any, Any}

Compute equivalent viscosity for multilayer model by recursively applying the formula for a halfspace and a channel from Lingle and Clark (1975).

get_rigidity(t, E, nu) -> Any

Compute rigidity D based on thickness t, Young modulus E and Poisson ration nu.

get_shearmodulus(youngmodulus, poissonratio) -> Any

Compute shear modulus G based on Young modulus E and Poisson ratio nu, or based on seismic velocities.

get_elastic_green(
    domain::RegionalDomain,
    greenintegrand_function,
    quad_support,
    quad_coeffs
) -> Any

Integrate load response over field by using 2D quadrature with specified support points and associated coefficients.

get_flexural_lengthscale(
    litho_rigidity,
    rho_uppermantle,
    g
) -> Any

Compute the flexural length scale, based on Coulon et al. (2021), Eq. in text after Eq. 3. The flexural length scale will be on the order of 100km.

Arguments

• litho_rigidity: Lithospheric rigidity
• rho_uppermantle: Density of the upper mantle
• g: Gravitational acceleration

Returns

• L_w: The calculated flexural length scale

get_relaxation_time(eta, m, p) -> Any

Convert the viscosity to relaxation times following Van Calcar et al. (in rev.).

get_relaxation_time_weaker(eta) -> Any

Compute the relaxation time for a weaker mantle, following Van Calcar et al. (in rev.).

get_relaxation_time_stronger(eta) -> Any

Compute the relaxation time for a stronger mantle, following Van Calcar et al. (in rev.).

absorption_band_density(τ, alpha, tau_L, tau_H) -> Any

Faul-Jackson power-law absorption-band density underlying the extended Burgers model (EBM) of Ivins & Caron (2021, §2): a probability density (∫[τ_L,τ_H] = 1) of retardation times over [tau_L, tau_H] with exponent alpha.

fit_prony_series(
;
    relaxation_strength,
    alpha,
    tau_L,
    tau_H,
    nbranches,
    ntest,
    nquad
)

Fit an nbranches-branch Prony series (Δⱼ, τⱼ) to the continuous Faul-Jackson absorption-band spectrum underlying the extended Burgers model (EBM) of Ivins & Caron (2021), so that the transient term of their creep function

relaxation_strength * ∫[tau_L, tau_H] F(τ) (1 − exp(−t/τ)) dτ

(F the density returned by absorption_band_density) is approximated by the nbranches-branch Kelvin sum Σⱼ Δⱼ (1 − exp(−t/τⱼ)) used by TransientCreepMantle.

nbranches log-spaced bins partition [tau_L, tau_H]. Each Δⱼ is the exact probability mass of its bin (so sum(Δⱼ) == relaxation_strength to machine precision for any nbranches, since the bins exactly partition the band) and τⱼ is the F-weighted mean retardation time within it — both closed-form (the density is a power law, so its CDF and first moment are elementary), no optimizer needed.

Returns (Δⱼ, τⱼ, fit_error), where fit_error is the maximum relative difference (against relaxation_strength) between the continuous and nbranches-branch transient term, sampled log-spaced over t ∈ [tau_L, tau_H] (ntest points; the continuous term is evaluated by Gauss-Legendre quadrature via quadrature1D with nquad points). Use fit_error to pick nbranches: I&C 2021 report that N ≈ 3-5 typically fits their spectrum to a few percent (roadmaps/burgers.md §6).

Example

julia> Δ, τ, err = fit_prony_series(relaxation_strength = 1.2, alpha = 0.5,
           tau_L = 1.0, tau_H = 100.0, nbranches = 4);

julia> sum(Δ)
1.2

julia> round.(τ, sigdigits = 3)
(1.98, 6.26, 19.8, 62.6)

julia> round(err, sigdigits = 2)
0.014

load_dataset(
    name::String;
    kwargs...
) -> Union{Nothing, Tuple{Union{Tuple{Any, Any}, Tuple{Any, Any, Any}, Vector{Float64}}, Any, Any}}

Return the dims::Tuple{Vararg{Vector}}, the field<:Array and the interpolator corresponding to a data set defined by a unique name::String. For instance:

(lon180, lat, t), Hice, Hice_itp = load_dataset("ICE6G_D")

Following options are available for parameter fields:

• "ICE6GD": ice loading history from ICE6GD.
• "Wiens2022": viscosity field from (Wiens et al. 2022)
• "Lithothickness_Pan2022": lithospheric thickness field from (Pan et al. 2022)
• "Viscosity_Pan2022": viscosity field from (Pan et al. 2022)

Following options are available for model results:

• "Spada2011"
• "LatychevGaussian"
• "LatychevICE6G"

Define the output NetCDF file for the simulation.

Fields

• t: a vector of time points at which the variable is defined.
• filename: the name of the NetCDF file.
• buffer: a buffer to store the output data before writing to the NetCDF file.
• vars3D: a vector of symbols representing the 3D variables to be output.
• vars1D: a vector of symbols representing the 1D variables to be output.
• params2D: a vector of symbols representing the 2D parameters to be output.
• oc: an instance of AbstractOutputCrop that defines the cropping strategy for the output.
• k: the current time step index for writing to the NetCDF file.

Define the native output which will be updated over the simulation.

Initialization example:

nout = NativeOutput(vars = [:u, :ue, :b, :dz_ss, :H_ice, :H_water, :u_x, :u_y],
    t = collect(0:1f3:10f3))

Fields

• t: a vector of time points at which the variable is defined.
• vars: a vector of symbols representing the variables to be output.
• vals: a dictionary mapping each variable symbol to a vector of matrices containing the output data.
• computation_time: the total computation time for the simulation.
• k: the current time step index for writing to the output.