Examples

Multi-layer Earth

FastIsostasy relies on a (polar) stereographic projection that allows to treat the radially-layered, onion-like structure of the solid Earth as a superposition of horizontal layers. Furthermore, FastIsostasy reduces this 3D problem into a 2D problem by collapsing the depth dimension, mainly through the computation of an effective viscosity field that accounts for the superposition of layers with different viscosities. The user is required to provide the 3D information, which will then be used under the hood to compute the effective viscosity. This tutorial shows such an example.

We want to render a situation similar to the one depicted below:

Initializing a LateralVariability with parameters corresponding to this situation automatically computes the conversion from a 3D to a 2D problem. This can be simply executed by running:

using FastIsostasy

W = 3000e3      # (m) half-width of the domain Wx = Wy
n = 7           # implies an Nx x Ny grid with Nx = Ny = 2^n = 64.
Omega = ComputationDomain(W, n)
c = PhysicalConstants(rho_litho = 0.0)

lv = [1e19, 1e21]       # (Pa s)
lb = [88e3, 400e3]      # (m)
p = LateralVariability(Omega, layer_viscosities = lv, layer_boundaries = lb)
extrema(p.effective_viscosity)

(5.939747413003433e20, 5.939747413003433e20)

As expected, the effective viscosity is a homogeneous field. It corresponds to a nonlinear mean of the layered values provided by the user. Note that we have set $\ŗho_\mathrm{litho} = 0$ to prevent the lithosphere from contributing to the hydrostatic, upward force. This is made to comply with the later computed analytical solution, which assumes a purely elastic lithosphere. In reality, this is however arguably wrong and the default choice c = PhysicalConstants() therefore uses $ \ŗho_\mathrm{litho} = 2600 \, \mathrm{kg \, m^{-3}} $.

The next section shows how to use the now obtained p::LateralVariability for actual GIA computation.

Simple load and geometry

We now apply a constant load, here a cylinder of ice with radius $ R = 1000 \, \mathrm{km} $ and thickness $H = 1 \, \mathrm{km}$, over Omega::ComputationDomain introduced in Multi-layer Earth. To formulate the problem conviniently, we use FastIsoProblem, a struct containing the variables and options that are necessary to perform the integration over time. We can then simply apply solve!(fip::FastIsoProblem) to perform the integration of the ODE. Under the hood, the ODE is obtained from the PDE by applying a Fourier collocation scheme contained in dudt_isostasy!. The integration is performed according to FastIsoProblem.diffeq::NamedTuple, which contains the algorithm and optionally tolerances, maximum iteration number... etc.

using CairoMakie

R = 1000e3                  # ice disc radius (m)
H = 1e3                     # ice disc thickness (m)
Hice = uniform_ice_cylinder(Omega, R, H)
t_out = years2seconds.([0.0, 200.0, 600.0, 2000.0, 5000.0, 10_000.0, 50_000.0])
interactive_sealevel = false
fip = FastIsoProblem(Omega, c, p, t_out, interactive_sealevel, Hice)
solve!(fip)

function plot3D(fip, k_idx)
    X, Y, out = Array(fip.Omega.X), Array(fip.Omega.Y), fip.out
    zl = extrema(out.ue[end] + out.u[end])
    fig = Figure(fontsize = 10)
    for j in eachindex(k_idx)
        ax = Axis3(fig[1, j])
        u_tot = out.ue[k_idx[j]] + out.u[k_idx[j]]
        surface!(ax, X, Y, u_tot, colormap = :cool)
        wireframe!(ax, X, Y, u_tot, color = :black, linewidth = 0.1)
        zlims!(ax, zl)
    end
    return fig
end
plot3D(fip, [lastindex(t_out) ÷ 2, lastindex(t_out)])

... and here goes the total displacement at $t = 50 \, \mathrm{kyr}$. You can now access the elastic and viscous displacement at time t_out[k] by respectively calling fip.out.ue[k] and fip.out.u[k]. For the present case, the latter can be compared to an analytic solution that is known for this particular case. Let's look at the accuracy of our numerical scheme over time by running following plotting commands:

fig = Figure()
ax = Axis(fig[1, 1])
cmap = cgrad(:jet, length(t_out), categorical = true)
ii, jj = Omega.Mx:Omega.Nx, Omega.My
x = Omega.X[ii, jj]
r = Omega.R[ii, jj]

# A support vector for computing the analytical solution
vsupport = vcat(1.0e-14, 10 .^ (-10:0.05:-3), 1.0)

for k in eachindex(t_out)
    analytic_solution_r(r) = analytic_solution(r, t_out[k], c, p, H, R, vsupport)
    u_analytic = analytic_solution_r.(r)
    u_numeric = fip.out.u[k][ii, jj]
    lines!(ax, x, u_analytic, color = cmap[k], linewidth = 5,
        label = L"$u_{ana}(t = %$(round(seconds2years(t_out[k]))) \, \mathrm{yr})$")
    lines!(ax, x, u_numeric, color = cmap[k], linewidth = 5, linestyle = :dash,
        label = L"$u_{num}(t = %$(round(seconds2years(t_out[k]))) \, \mathrm{yr})$")
end
axislegend(ax, position = :rb, nbanks = 2, patchsize = (50.0f0, 20.0f0))
fig

GPU support

For about $n > 7$, the previous example can be computed even faster by using GPU parallelism. It could not represent less work from the user's perspective, as it boils down to calling ComputationDomain with an extra keyword argument and passing it to a ::LateralVariability with the viscosity and depth values defined earlier:

n = 8
Omega = ComputationDomain(W, n, use_cuda = true)
p = LateralVariability(Omega, layer_viscosities = lv, layer_boundaries = lb)
Hice = uniform_ice_cylinder(Omega, R, H)
fip = FastIsoProblem(Omega, c, p, t_out, interactive_sealevel, Hice)
solve!(fip)
plot3D(fip, [lastindex(t_out) ÷ 2, lastindex(t_out)])

That's it, nothing more! For postprocessing, consider using reinit_structs_cpu.

Make your own time loop

As any high-level function, solve! has some limitations. An ice-sheet modeller typically wants to embed FastIsostasy within a time-stepping loop. This can be easily done by getting familiar with some intermediate-level functions like init, step! and write_out!:

Omega = ComputationDomain(3000e3, 6)
p = LateralVariability(Omega)
Hice = uniform_ice_cylinder(Omega, R, H)
fip = FastIsoProblem(Omega, c, p, t_out, interactive_sealevel, Hice)

update_diagnostics!(fip.geostate.dudt, fip.geostate.u, fip, 0.0)
write_out!(fip, 1)
ode = init(fip)
@inbounds for k in eachindex(fip.out.t)[2:end]
    step!(fip, ode, (fip.out.t[k-1], fip.out.t[k]))
    write_out!(fip, k)
end
plot3D(fip, [lastindex(t_out) ÷ 2, lastindex(t_out)])

Antarctic deglaciation

We now want to provide an example that presents:

• a heterogeneous lithosphere thickness
• a heterogeneous upper-mantle viscosity
• various viscous channels
• a more elaborate load that evolves over time
• changes in the sea-level

For this we run a deglaciation of Antarctica with lithospheric thickness and upper-mantle viscosity from ^[Wiens2021] and the ice thickness history from ^[Briggs2014]. Since the load is known and the isostatic response does not influence it (one-way coupling), we can provide snapshots of the ice thickness and their associated time to FastIsoProblem. Under the hood, an interpolator is created and called within the time integration.

# Code is coming soon!

Inversion of viscosity field

FastIsostasy relies on simplifications of the full GIA problem and might therefore need a calibration step to match the data, typically obtained from observations or from a "golden-standard" 3D GIA model. By means of an unscented Kalman inversion, one can, for instance, infer the appropriate field of effective mantle viscosity that matches the data best. Whereas this is known to be a tedious step, FastIsostasy is developped to ease the procedure by providing a convenience struct InversionProblem. We demonstrate this on a low-resolution grid since (1) the underlying unscented Kalman filter requires many simulations and (2) estimating high-resolution viscosity field might lead to overfit the problem. The effective viscosity field we estimate^[Wiens2021] can be loaded by using load_wiens2021 with appropriate depths of the layer boundaries:

Omega = ComputationDomain(3000e3, 5)
c = PhysicalConstants()
lb = [88e3, 180e3, 280e3, 400e3]
lv = load_wiens2021(Omega)
p = LateralVariability(Omega, layer_boundaries = lb, layer_viscosities = lv)

LateralVariability{Float64, Matrix{Float64}}([3.670395852558643e20 7.002434790938157e20 … 1.052015201260212e21 1.0425684675203643e21; 5.3281592663754165e20 4.51145048174098e20 … 1.0348145886822903e21 1.0112398581332767e21; … ; 7.503705697052097e20 8.676125521740755e20 … 1.0828624808347536e21 1.0624202481101908e21; 5.744571315636416e20 8.545566694178897e20 … 1.0789331175288494e21 1.0748586847609106e21], [88000.0 88000.0 … 88000.0 88000.0; 88000.0 88000.0 … 88000.0 88000.0; … ; 88000.0 88000.0 … 88000.0 88000.0; 88000.0 88000.0 … 88000.0 88000.0], [4.0669444444444445e24 4.0669444444444445e24 … 4.0669444444444445e24 4.0669444444444445e24; 4.0669444444444445e24 4.0669444444444445e24 … 4.0669444444444445e24 4.0669444444444445e24; … ; 4.0669444444444445e24 4.0669444444444445e24 … 4.0669444444444445e24 4.0669444444444445e24; 4.0669444444444445e24 4.0669444444444445e24 … 4.0669444444444445e24 4.0669444444444445e24], 0.28, 0.28, [1.0000713826855053e17 4.887401299882576e17 … 9.769469918922297e18 2.0879495459931455e19; 2.3072199590312496e17 1.5586550105845824e17 … 6.973262275191022e18 8.023576542337122e18; … ; 6.221541465964311e17 1.1680542124861146e18 … 5.417351309025628e19 1.7470566997634925e19; 2.7742014858401277e17 1.0847783276129215e18 … 2.1177856440349252e20 1.2824900370343915e20;;; 2.874928181083441e19 2.4246887642920042e20 … 7.411139703483182e19 1.7100980036606878e19; 2.958475655940619e19 3.589021004339394e19 … 4.367440529494825e19 1.390384763742017e19; … ; 1.1488729389531903e20 1.6638493446882696e20 … 1.1997340918082454e20 4.2287596319460565e19; 8.03195006905522e19 1.2334143589099081e20 … 6.434339912110696e19 4.13404309897143e19;;; 3.664375746478333e20 3.07560812831034e21 … 5.4508471007474935e20 7.480196385395158e19; 1.1567113579074315e21 3.461266188711196e20 … 4.591028618714835e20 5.742513171477367e19; … ; 2.308718590351057e21 9.644794974369998e20 … 5.919595072485586e20 4.8772273440839015e20; 1.3479755037458741e21 6.939880623602302e20 … 1.1781890985100077e20 1.3603635913581722e20;;; 1.0e21 1.0e21 … 1.0e21 1.0e21; 1.0e21 1.0e21 … 1.0e21 1.0e21; … ; 1.0e21 1.0e21 … 1.0e21 1.0e21; 1.0e21 1.0e21 … 1.0e21 1.0e21], [88000.0 88000.0 … 88000.0 88000.0; 88000.0 88000.0 … 88000.0 88000.0; … ; 88000.0 88000.0 … 88000.0 88000.0; 88000.0 88000.0 … 88000.0 88000.0;;; 180000.0 180000.0 … 180000.0 180000.0; 180000.0 180000.0 … 180000.0 180000.0; … ; 180000.0 180000.0 … 180000.0 180000.0; 180000.0 180000.0 … 180000.0 180000.0;;; 280000.0 280000.0 … 280000.0 280000.0; 280000.0 280000.0 … 280000.0 280000.0; … ; 280000.0 280000.0 … 280000.0 280000.0; 280000.0 280000.0 … 280000.0 280000.0;;; 400000.0 400000.0 … 400000.0 400000.0; 400000.0 400000.0 … 400000.0 400000.0; … ; 400000.0 400000.0 … 400000.0 400000.0; 400000.0 400000.0 … 400000.0 400000.0])

To make this problem more exciting, we shift the center of the ice load to $ (-1000, -1000) \: \mathrm{km} $ where the viscosity field displays a less uniform structure. For the sake of simplicity, the data to fit is obtained from a FastIsostasy simulation with the ground-truth viscosity field.

R, H = 1000e3, 1e3
Hice = uniform_ice_cylinder(Omega, R, H, center = [-1000e3, -1000e3])
t_out = years2seconds.(1e3:1e3:2e3)
fip = FastIsoProblem(Omega, c, p, t_out, false, Hice)
solve!(fip)
true_viscosity = copy(p.effective_viscosity)

32×32 Matrix{Float64}:
 3.6704e20   7.00243e20  9.71794e20  …  1.06319e21  1.05202e21  1.04257e21
 5.32816e20  4.51145e20  7.40019e20     1.04048e21  1.03481e21  1.01124e21
 5.57151e20  4.16429e20  7.66751e20     9.32778e20  8.94699e20  9.31377e20
 9.1139e20   6.53154e20  7.00812e20     9.26984e20  9.6971e20   9.5188e20
 3.3541e20   1.25592e20  3.35984e20     9.06864e20  9.34526e20  9.1974e20
 8.28284e20  3.88763e20  3.39259e20  …  8.07857e20  6.3595e20   4.43712e20
 9.05803e20  9.81046e20  8.24306e20     8.01611e20  4.31973e20  2.52081e20
 7.8454e20   9.66944e20  9.4441e20      1.03061e21  9.56666e20  1.09844e21
 5.00079e20  9.17061e20  1.0228e21      1.00614e21  1.01645e21  1.07899e21
 8.30033e20  1.0151e21   1.043e21       9.20827e20  9.81348e20  1.04752e21
 ⋮                                   ⋱              ⋮           
 1.22891e20  9.065e20    1.081e21       1.08645e21  1.0922e21   1.09126e21
 3.97022e20  9.80223e20  1.00115e21     1.08934e21  1.09384e21  1.09371e21
 8.7232e20   1.00963e21  1.05744e21  …  1.08583e21  1.09006e21  1.08935e21
 7.9576e20   9.01289e20  1.03423e21     1.09219e21  1.09222e21  1.09094e21
 8.745e20    9.29279e20  1.0073e21      1.09271e21  1.09518e21  1.09452e21
 9.21016e20  9.70803e20  1.07117e21     1.09079e21  1.09503e21  1.09467e21
 8.61142e20  9.51535e20  1.07658e21     1.08915e21  1.09217e21  1.08869e21
 7.50371e20  8.67613e20  1.06016e21  …  1.08969e21  1.08286e21  1.06242e21
 5.74457e20  8.54557e20  1.05923e21     1.08163e21  1.07893e21  1.07486e21

Now that we have the displacement field, we can recover the viscosity field from which it results. We therefore pass an InversionConfig and an InversionData to an InversionProblem. Let's look at the initialized viscosity field:

config = InversionConfig(N_iter = 15)
data = InversionData(copy(fip.out.t[2:end]), copy(fip.out.u[2:end]), copy([Hice, Hice]), config)
paraminv = InversionProblem(deepcopy(fip), config, data)

function plot_viscfields(paraminv)
    estim_viscosity = copy(true_viscosity)
    estim_viscosity[paraminv.data.idx] .= 10 .^ get_ϕ_mean_final(paraminv.priors, paraminv.ukiobj)

    cmap = cgrad(:jet, rev = true)
    crange = (20, 21.2)
    fig = Figure()
    axs = [Axis(fig[1,i], aspect = DataAspect()) for i in 1:2]
    heatmap!(axs[1], log10.(true_viscosity), colormap = cmap, colorrange = crange)
    heatmap!(axs[2], log10.(estim_viscosity), colormap = cmap, colorrange = crange)
    Colorbar(fig[2, :], vertical = false, colormap = cmap, colorrange = crange,
        width = Relative(0.5))
    return fig
end
plot_viscfields(paraminv)

By performing a Kalman inversion, we can achieve a close match between ground truth and estimated viscosity field:

solve!(paraminv)
plot_viscfields(paraminv)

This remains an academic example, where we try to recover a known parameter field from data generated by the model itself. Nonetheless, the user should get a proof-of-concept and a scheme of how to implement such a procedure themself.