Tutorial

In this section, we will present some examples with idealised loads and solid-Earth parameters. This should give the user a basic understanding of FastIsostasy's basic functions.

3D ➡ 2D Earth

FastIsostasy relies on a (polar) stereographic projection. Let's first create Omega::ComputationDomain and visualise how this relates to a domain on a spherical Earth:

using CairoMakie, FastIsostasy

W = 3000e3      # (m) half-width of the domain Wx = Wy
n = 7           # implies an Nx x Ny grid with Nx = Ny = 2^n = 128.
Omega = ComputationDomain(W, n)
fig = Figure(size = (1600, 800), fontsize = 24)
ax1 = Axis3(fig[1, 1], title = "Original grid")
ax2 = Axis3(fig[1, 2], title = "Projected grid")
wireframe!(ax1, Omega.X .* Omega.K, Omega.Y .* Omega.K,
    Omega.R .* cos.(deg2rad.(Omega.Lat)), color = :gray10, linewidth = 0.1)
wireframe!(ax2, Omega.X .* Omega.K, Omega.Y .* Omega.K,
    Omega.null, color = :gray10, linewidth = 0.1)
for ax in [ax1, ax2]
    zlims!(ax, (0, 5e6))
    hidedecorations!(ax)
    hidespines!(ax)
end
fig

The projection 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 LayeredEarth with parameters corresponding to this situation automatically computes the conversion from a 3D to a 2D problem. Since we will compare our solution to an analytical one of a flat Earth, we exceptionally switch off the distortion correction, which accounts for the distortion factor Omega.K in the computation. This can be simply executed by running:

Omega = ComputationDomain(W, n, correct_distortion = false)
c = PhysicalConstants()
lv = [1e19, 1e21]       # viscosity layers (Pa s)
lb = [88e3, 400e3]      # depth of layer boundaries (m)
p = LayeredEarth(Omega, layer_viscosities = lv, layer_boundaries = lb, rho_litho = 0.0)
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 $\rho_{litho} = 0$ to prevent the lithosphere from contributing to the hydrostatic pressure term, such that the numerical solution obtained here is comparable with the analytical one provided in Bueler et al. (2007). The default value however yields $\rho_{litho} = 3200 \, \mathrm{kg/m^3}$ and contributes to the hydrostatic pressure term.

The next section shows how to use the now obtained p::LayeredEarth 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 LayeredEarth. 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 lv_elva!. The integration is performed according to FastIsoProblem.diffeq::NamedTuple, which contains the algorithm and optionally tolerances, maximum iteration number... etc.

t_out = [0.0, 2e2, 6e2, 2e3, 5e3, 1e4, 5e4]     # vector of output time steps (yr)
εt = 1e-8
pushfirst!(t_out, -εt)                          # append step to have Heaviside at t=0

R = 1000e3                                      # ice disc radius (m)
H = 1e3                                         # ice disc thickness (m)
Hcylinder = uniform_ice_cylinder(Omega, R, H)   # field representing ice disk

t_Hice = [-εt, 0.0, t_out[end]]                 # ice history = Heaviside at t=0
Hice = [zeros(Omega.Nx, Omega.Ny), Hcylinder, Hcylinder]

fip = FastIsoProblem(Omega, c, p, t_out, t_Hice, Hice, output = "sparse")
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
fig = plot3D(fip, [lastindex(t_out) ÷ 2, lastindex(t_out)])

The figure above shows the total displacement at $t = 0.6 \, \mathrm{kyr}$ and $t = 50 \, \mathrm{kyr}$. Since we defined output = "sparse", we can now access the elastic and viscous displacement at time t_out[k] by 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]

for k in eachindex(t_out)[2:end]
    analytic_solution_r(r) = analytic_solution(r, years2seconds(t_out[k]),
        c, p, H, R)
    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 = %$(t_out[k]) \, \mathrm{yr})$")
    lines!(ax, x, u_numeric, color = cmap[k], linewidth = 5, linestyle = :dash,
        label = L"$u_{num}(t = %$(t_out[k]) \, \mathrm{yr})$")
end
axislegend(ax, position = :rb, patchsize = (50.0f0, 20.0f0))
fig

GPU support

For about $n \geq 7$, the present 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:

Omega = ComputationDomain(W, n, use_cuda = true)

We then pass Omega to a LayeredEarth and a FastIsoProblem, which we solve as done above: that's it!

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, n)
p = LayeredEarth(Omega)
fip = FastIsoProblem(Omega, c, p, t_out, t_Hice, Hice, output = "sparse")

update_diagnostics!(fip.now.dudt, fip.now.u, fip, 0.0)
write_out!(fip.out, fip.now, 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.out, fip.now, k)
end
fig = plot3D(fip, [lastindex(t_out) ÷ 2, lastindex(t_out)])

Using different deformation models

ELRA is a GIA model that is commonly used in ice-sheet modelling. For the vast majority of applications, it is less accurate than LV-ELVA without providing any significant speed up (Swierczek-Jereczek et al., 2024). However, it can be used by specifying adequate options:

p = LayeredEarth(Omega, tau = 3e3)
opts = SolverOptions(deformation_model = :elra)
fip = FastIsoProblem(Omega, c, p, t_out, t_Hice, Hice, opts = opts, output = "sparse")
solve!(fip)
fig = plot3D(fip, [lastindex(t_out) ÷ 2, lastindex(t_out)])

The reader may have noticed that the equilibrium displacement (at $t = 50 \, \mathrm{kyr}$) given by ELRA is the same as the one given by LV-ELVA, although their transient behaviour differ (e.g. at $t = 0.6 \, \mathrm{kyr}$). This is expected when the lithosphere yields a constant thickness but breaks down when the lithosphere thickness varies.