Module Data

The module Data is created in the module where @init_parallel_stencil is called from. It provides the following types:

Data.Number

The type of numbers used by @zeros, @ones, @rand and @fill and in all array types of module Data (selected with argument numbertype of @init_parallel_stencil).

Data.Array{ndims}

Expands to Data.Array{numbertype, ndims}, where numbertype is the datatype selected with @init_parallel_stencil and the datatype Data.Array is chosen to be compatible with the package for parallelization selected with @init_parallel_stencil (Array for Threads, CUDA.CuArray or CUDA.CuDeviceArray for CUDA and AMDGPU.ROCArray or AMDGPU.ROCDeviceArray for AMDGPU; @parallel and @parallel_indices convert CUDA.CuArray and AMDGPU.ROCArray automatically to CUDA.CuDeviceArray and AMDGPU.ROCDeviceArray in kernels when required).

Data.CellArray{ndims}

Expands to Data.CellArray{numbertype, ndims}, where numbertype is the datatype selected with @init_parallel_stencil and the datatype Data.CellArray is chosen to be compatible with the package for parallelization selected with @init_parallel_stencil (CPUCellArray for Threads, CuCellArray or CuDeviceCellArray for CUDA and ROCCellArray or ROCDeviceCellArray for AMDGPU; @parallel and @parallel_indices convert CellArray automatically to DeviceCellArray when required).

Data.Cell{S}

Expands to Union{StaticArrays.SArray{S, numbertype}, StaticArrays.FieldArray{S, numbertype}}, where numbertype is the datatype selected with @init_parallel_stencil.

Data.NumberTuple{N_tuple} | Data.NamedNumberTuple{N_tuple, names} | Data.NumberCollection{N_tuple}

Expands to: NTuple{N_tuple, Data.Number} | NamedTuple{names, NTuple{N_tuple, Data.Number}} | Union{Data.NumberTuple{N_tuple}, Data.NamedNumberTuple{N_tuple}}

Data.ArrayTuple{N_tuple, N} | Data.NamedArrayTuple{N_tuple, N, names} | Data.ArrayCollection{N_tuple, N}

Expands to: NTuple{N_tuple, Data.Array{N}} | NamedTuple{names, NTuple{N_tuple, Data.Array{N}}} | Union{Data.ArrayTuple{N_tuple, N}, Data.NamedArrayTuple{N_tuple, N}}

Data.CellArrayTuple{N_tuple, N, B} | Data.NamedCellArrayTuple{N_tuple, N, B, names} | Data.CellArrayCollection{N_tuple, N, B}

Expands to: NTuple{N_tuple, Data.CellArray{N, B}} | NamedTuple{names, NTuple{N_tuple, Data.CellArray{N, B}}} | Union{Data.CellArrayTuple{N_tuple, N, B}, Data.NamedCellArrayTuple{N_tuple, N, B}}

Data.CellTuple{N_tuple, S} | Data.NamedCellTuple{N_tuple, S, names} | Data.CellCollection{N_tuple, S}

Expands to: NTuple{N_tuple, Data.Cell{S}} | NamedTuple{names, NTuple{N_tuple, Data.Cell{S}}} | Union{Data.CellTuple{N_tuple, S}, Data.NamedCellTuple{N_tuple, S}}

Data.DeviceArray{ndims}

Data.DeviceCellArray{ndims}

ForwardPlan

Allias for in-place precomputed plans from FFTW or CUFFT. Used to compute forward FFT.

InversePlan

Allias for in-place precomputed plans from FFTW or CUFFT. Used to compute inverse FFT.

ComputationDomain
ComputationDomain(W, n)
ComputationDomain(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:

Omega = ComputationDomain(W, n)

If a rectangular domain is needed, run:

Omega = ComputationDomain(Wx, Wy, Nx, Ny)

CurrentState

Return a mutable struct containing the geostate which will be updated over the simulation. The geostate contains all the states of the [FastIsoProblem] to be solved.

DenseOutputs()

Return a struct containing the fields of viscous displacement, viscous displacement rate, elastic displacement, geoid displacement, sea level and the computation time resulting from solving a FastIsoProblem.

FastIsoProblem(Omega, c, p, t_out)
FastIsoProblem(Omega, c, p, t_out, Hice)
FastIsoProblem(Omega, c, p, t_out, t_Hice, Hice)

Return a struct containing all the other structs needed for the forward integration of the model over Omega::ComputationDomain with parameters c::PhysicalConstants and p::LayeredEarth. The outputs are stored at t_out::Vector{<:AbstractFloat}.

FastIsoTools(Omega, c, p)

Return a struct containing pre-computed tools to perform forward-stepping of the model. This includes the Green's functions for the computation of the lithosphere and geoid displacement, plans for FFTs, interpolators of the load and the viscosity over time and preallocated arrays.

InversionConfig

Struct containing configuration parameters for a [InversionProblem].

Need to choose regularization factor α ∈ (0,1], When you have enough observation data α=1: no regularization

update_freq 1 : approximate posterior cov matrix with an uninformative prior 0 : weighted average between posterior cov matrix with an uninformative prior and prior

InversionData

Struct containing data (either observational or output of a golden standard model) for a [InversionProblem].

InversionProblem

Struct containing variables and configs for the inversion of Solid-Earth parameter fields. For now, only viscosity can be inverted but future versions will support lithosphere rigidity. For now, the unscented Kalman inversion is the only method available but ensemble Kalman inversion will be available in future.

LayeredEarth(Omega; layer_boundaries, layer_viscosities)

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

Omega = ComputationDomain(3000e3, 7)
lb = [100e3, 300e3]
lv = [1e19, 1e21]
p = LayeredEarth(Omega, 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 = 200 \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 Omega::ComputationDomain.

OceanSurfaceChange(; z0 = 0.0)

Return a mutable struct OceanSurfaceChange containing:

• z_k: the GMSL at current time step k.
• A_k: the ocean surface at current time step k.
• z: a vector of GMSL values used as knots for interpolation.
• A: a vector of ocean surface values used as knots for interpolation.
• A_itp: an interpolator of ocean surface over depth. Bias-free for present-day.
• A_pd: the present-day ocean surface.
• res: residual of the nonlinear equation solved numerically.

An osc::OceanSurfaceChange can be used as function to update osc.z_k and osc.A_k based on osc.A_itp and an input delta_V by running:

osc(delta_V)

PhysicalConstants

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

ReferenceEarthModel

Return a struct with vectors containing the:

• radius (distance from Earth center),
• depth (distance from Earth surface),
• density,
• P-wave velocities,
• S-wave velocities,

which are typically used to characterize the properties of a spherically symmetrical solid Earth.

ReferenceState

Return a struct containing the reference state.

Options

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

analytic_solution(r, t, c, p, H0, R0, domains; n_quad_support)

Return the analytic solution of the bedrock displacement resulting from a cylindrical ice load with radius R0 and height H0 on a flat Earth represented by an elastic plate overlaying a viscous half space. Parameters are provided in c, p. The points at which the solution is computed are specified by the distance r to the center of the domain. The time at which the solution is computed is specified by t.

apply_bc!(u::M, bcm::M, nbc::T)

A generic function to update u such that the boundary conditions are respected on average, which is the only way to impose them for Fourier collocation.

build_greenintegrand(distance::Vector{T}, 
    greenintegrand_coeffs::Vector{T}) where {T<:AbstractFloat}

Compute the integrands of the Green's function resulting from a load at a given distance and based on provided greenintegrand_coeffs. Reference: Deformation of the Earth by surface Loads, Farell 1972, table A3.

calc_viscous_green(GV, kei2D, L_w, D_lith, dx, dy)

Calculate the viscous Green's function. Note that L_w contains information about the density of the upper mantle.

dist2angulardist(r::Real)

Convert Euclidean to angular distance along great circle.

elra!(dudt, u, fip, t)

Update the displacement rate dudt of the viscous response according to ELRA.

extract_inversion()

Extract results of parameter inversion from the priors and ukiobj that resulted from solve!(paraminv::InversionProblem).

gauss_distr(X::KernelMatrix{T}, Y::KernelMatrix{T},
    mu::Vector{<:Real}, sigma::Matrix{<:Real})

Compute Z = f(X,Y) with f a Gaussian function parametrized by mean mu and covariance sigma.

get_differential_fourier(W, N2)

Compute the matrices representing the differential operators in the fourier space.

get_effective_viscosity(
    layer_viscosities::Vector{KernelMatrix{T}},
    layers_thickness::Vector{T},
    Omega::ComputationDomain{T, M},
) where {T<:AbstractFloat}

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

get_elasticgreen(Omega, quad_support, quad_coeffs)

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

get_flexural_lengthscale(litho_rigidity, rho_uppermantle, g)

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_geoidgreen(fip::FastIsoProblem)

Return the Green's function used to compute the geoid anomaly as in ^[Coulon2021].

get_greenintegrand_coeffs(T)

Return the load response coefficients with type T. Reference: Deformation of the Earth by surface Loads, Farell 1972, table A3.

get_kei(filt, L_w, dx, dy; T = Float64)

Calculate the Kelvin filter in 2D.

Arguments

• filt::Matrix: The filter array to be filled with Kelvin function values.
• L_w::Real: The characteristic length scale.
• dx::Real: The grid spacing in the x-direction.
• dy::Real: The grid spacing in the y-direction.

Returns

• filt::Matrix{T}: The filter array filled with Kelvin function values.

get_kei_value(r, L_w, rn_vals, kei_vals)

Calculate the Kelvin function (kei) value based on the radius from the point load r, the flexural length scale L_w, and the arrays of normalized radii rn_vals and corresponding kei values kei_vals.

This function first normalizes the radius r by the flexural length scale L_w to get the current normalized radius. If this value is greater than the maximum value in rn_vals, the function returns the maximum value in kei_vals. Otherwise, it finds the interval in rn_vals that contains the current normalized radius and performs a linear interpolation to calculate the corresponding kei value.

get_loadgreen(r::T, rm::Vector{T}, greenintegrand_coeffs::Vector{T},     
    interp_greenintegrand_::Interpolations.Extrapolation) where {T<:AbstractFloat}

Compute the integrands of the Green's function resulting from a load at a given distance and based on provided greenintegrand_coeffs. Reference: Deformation of the Earth by surface Loads, Farell 1972, table A3.

get_normalized_lin_transform(x1, x2)

Return parameters of linear function mapping x1, x2 onto -1, 1.

get_quad_coeffs(T, n)

Return support points and associated coefficients with specified Type for Gauss-Legendre quadrature.

get_r(x::T, y::T) where {T<:Real}

Get euclidean distance of point (x, y) to origin.

get_rigidity(t::T, E::T, nu::T) where {T<:AbstractFloat}

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

init(fip)

Initialize an ode::CoupledODEs, aimed to be used in step!.

kernelpromote(X, arraykernel)

Promote X to the kernel (Array or CuArray) specified by arraykernel.

latlon2stereo(lat, lon, lat0, lon0)

Compute stereographic projection (x,y) for a given latitude lat longitude lon, reference latitude lat0 and reference longitude lon0. Optionally one can provide lat::KernelMatrix and lon::KernelMatrix if the projection is to be computed for the whole domain. Note: angles must be provided in degrees! Reference: John P. Snyder (1987), p. 157, eq. (21-2), (21-3), (21-4).

load_dataset(name) → (dims), field, interpolator

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"

load_prem()

Load Preliminary Reference Earth Model (PREM) from Dzewonski and Anderson (1981).

get_greenintegrand_coeffs(T)

Return the load response coefficients with type T. Reference: Deformation of the Earth by surface Loads, Farell 1972, table A3.

lon360tolon180(lon, X)

Convert longitude and field from lon=0:360 to lon=-180:180.

lv_elra!(dudt, u, fip, t)

Update the displacement rate dudt of the viscous response according to LV-ELRA.

lv_elva!(dudt, u, fip, t)

Update the displacement rate dudt of the viscous response according to LV-ELVA.

m_per_sec2mm_per_yr(dudt::Real)

Convert displacement rate dudt from $ m \, s^{-1} $ to $ mm \, \mathrm{yr}^{-1} $.

matrify(x, Nx, Ny)

Generate a vector of constant matrices from a vector of constants.

meshgrid(x, y)

Return a 2D meshgrid spanned by x, y.

normalized_lin_transform(y, m, p)

Apply normalized linear transformation with slope m and bias p on y.

quadrature1D(f, n, x1, x2)

Compute 1D Gauss-Legendre quadrature of f between x1 and x2 based on n support points.

quadrature2D(f, x, w, x1, x2, y1, y2)

Return the integration of f over [x1, x2] x [y1, y2] with x, w some pre-computed support points and coefficients of the Gauss-Legendre quadrature.

reinit_structs_cpu(Omega, p)

Reinitialize Omega::ComputationDomain and p::LayeredEarth on the CPU, mostly for post-processing purposes.

relax_u_eq!(fip)

Get equilibrium displacement u_eq by integrating LV-ELVA equation, which gives the same equilibrium as LV-ELRA.

samesize_conv_indices(N, M)

Get the start and end indices required for a samesize_conv

savefip(filename, fip; T = Float32)

Save the output of fip::FastIsoProblem as NetCDF file under filename.

scalefactor(lat::T, lon::T, lat0::T, lon0::T) where {T<:Real}
scalefactor(lat::M, lon::M, lat0::T, lon0::T) where {T<:Real, M<:KernelMatrix{T}}

Compute scaling factor of stereographic projection for a given (lat, lon) and reference (lat0, lon0). Angles must be provided in radians. Reference: ^[Snyder1987], p. 157, eq. (21-4).

seconds2years(t::Real)

Convert input time t from seconds to years.

simple_euler!()

Update the state u by performing an explicit Euler integration of its derivative dudt over a time step dt.

solve!(fip)

Solve the isostatic adjustment problem defined in fip::FastIsoProblem.

solve!(paraminv::InversionProblem)

Return priors and ukiobj that allow to extract the results of the parameter inversion as initialized in paraminv.

step!(fip)

Step fip::FastIsoProblem over tspan and based on ode::CoupledODEs, typically obtained by init.

stereo2latlon(x, y, lat0, lon0)

Compute the inverse stereographic projection (lat, lon) based on Cartesian coordinates (x,y) and for a given reference latitude lat0 and reference longitude lon0. Optionally one can provide x::KernelMatrix and y::KernelMatrix if the projection is to be computed for the whole domain. Note: angles must be para elloprovided in degrees!

Convert stereographic (x,y)-coordinates to latitude-longitude. Reference: John P. Snyder (1987), p. 159, eq. (20-14), (20-15), (20-18), (21-15).

surfacechange_residual!(osc, z, delta_V)
surfacechange_residual!(Vresidual, z, z_k, A_itp, delta_V)

Update the residual of the piecewise linear approximation, used to solve the sea-level/ocean-surface nonlinearity.

three_layer_scaling(Omega::ComputationDomain, kappa::T, visc_ratio::T,
    channel_thickness::T)

Return the viscosity scaling for a three-layer model and based on a the wave number kappa, the visc_ratio and the channel_thickness. Reference: Bueler et al. 2007, below equation 15.

update_V_af!(fip::FastIsoProblem)

Update the volume contribution from ice above floatation as in goelzer-brief-2020 (eq. 13). Note: we do not use (eq. 1) as it is only a special case of (eq. 13) that does not allow a correct representation of external sea-level forcings.

update_V_den!(fip::FastIsoProblem)

Update the volume contribution associated with the density difference between meltwater and sea water, as in goelzer-brief-2020 (eq. 10).

update_V_pov!(fip::FastIsoProblem)

Update the volume contribution to the ocean (from isostatic adjustement in ocean regions), which corresponds to the "potential ocean volume" in goelzer-brief-2020 (eq. 14). Note: we do not use eq. (8) as it is only a special case of eq. (14) that does not allow a correct representation of external sea-level forcings.

update_bsl!(fip::FastIsoProblem)

Update the sea-level contribution of melting above floatation and density correction. Note that this differs from goelzer-brief-2020 (eq. 12) because the ocean surface is not assumed to be constant. Furthermore, the contribution to ocean volume from the bedrock uplift is not included here since the volume displaced on site is arguably blanaced by the depression of the peripherial forebulge.

update_deformation_rhs!(fip)

Update the right-hand side of the deformation equation.

update_diagnostics!(dudt, u, fip, t)

Update all the diagnotisc variables, i.e. all fields of fip.now apart from the displacement, which requires an integrator.

update_elasticresponse!(fip::FastIsoProblem)

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

update_geoid!(fip::FastIsoProblem)

Update the geoid by convoluting the Green's function with the load anom.

update_loadcolumns!(fip::FastIsoProblem, u::KernelMatrix{T}, H_ice::KernelMatrix{T})

Update the load columns of a ::CurrentState.

update_seasurfaceheight!(fip::FastIsoProblem)

Update the sea-level by adding the various contributions as in coulon-contrasting-2021. Here, the constant term is used to impose a zero geoid perturbation in the far field rather thank for mass conservation and is embedded in convolution operation.

where_significant()

Find points of parameter field that can be inverted. We here assume that

write_out!(fip::FastIsoProblem)

Write results in output vectors if the load is updated internally. If the load is updated externally, the user is responsible for writing results.

years2seconds(t::Real)

Convert input time t from years to seconds.