Public API

Simulation(domain, c, solidearth, t_out)
Simulation(domain, c, solidearth, t_out, Hice)
Simulation(domain, c, solidearth, t_out, t_Hice, Hice)

Return a struct containing all the other structs needed for the forward integration of the model over domain::RegionalDomain with parameters c::PhysicalConstants and solidearth::SolidEarth. The outputs are stored at t_out::Vector{<:AbstractFloat}. To perform the whole simulation, run!(sim::Simulation) and to perform a single step:

tspan = (t_start, t_end)
integrator = init_integrator(sim)
step!(integrator, tspan)

Return a struct containing the options relative to solving a Simulation.

Contains:

• alg::ODEsolvers: the algorithm to integrate the ODE forward in time.
• reltol: the relative error tolerance of the integrator.

run!(sim::Simulation)

Solve the isostatic adjustment problem defined in sim::Simulation.

init_integrator(sim::Simulation) -> Any

Initialise the integrator of sim::Simulation, which can be subsequently integrated forward in time by using step!.

Return a struct containing important physical constants. Comes with default values that can however be changed by the user, for instance by running:

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

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

Abstract type for domain representation in the model.

RegionalDomain
RegionalDomain(W, n)
RegionalDomain(Wx, Wy, nx, ny)

Return a struct containing all information related to geometry of the domain and potentially used parallelism. To initialize one with 2*W and 2^n grid cells:

domain = RegionalDomain(W, n)

If a rectangular domain is needed, run:

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

Not implemented yet!

Version 3.0 will allow for global domains.

A struct containing the boundary conditions of the problem:

• 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.

An abstract type that determines how the ice thickness is updated in the model. This is done by implementing the update_ice! function for different subtypes. Available subtypes are:

• TimeInterpolatedIceThickness
• ExternallyUpdatedIceThickness

A struct to update the ice thickness based on time interpolation. Contains:

• 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.

ExternallyUpdatedIceThickness

A struct to indicate that the ice thickness is updated 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 are:

• 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 are:

• OffsetBC

A boundary condition that applies an offset to the values at the boundaries of a computational domain. Contains:

• 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 AbstractBCRule.

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 are:

• 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)

ReferenceBSL(z = 0; T = Float32, itp_kwargs = (extrapolation_bc = Flat()))

A struct used in all subtypes of AbstractBSL to define a reference of barystatic sea level and ocean surface area. Contains:

• 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 are:

• InternalBSLUpdate
• ExternalBSLUpdate

A struct to indicate that the BSL is updated internally by the model based on the change in ice volume.

A struct to indicate that the BSL is updated 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.

Available subtypes are:

• RigidMantle
• RelaxedMantle
• MaxwellMantle

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. This is the most realistic mantle model and generally offers the best performance. It is the default mantle model used in the solver.

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

Update the time derivative of the viscous displacement dudt based on an AbstractMantle:

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

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
) -> Matrix{T} where T<:AbstractFloat

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.).

load_dataset(
    name::String;
    kwargs...
) -> Union{Nothing, Tuple{Any, 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"

A struct that contains all the necessary information to store the output in a NetCDF file.

Can be initilized as:

Return a mutable struct containing 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))