GitHub

Glacial cycle

We now want to provide an example that presents:

  • a heterogeneous lithosphere thickness,
  • a heterogeneous upper-mantle viscosity,
  • a stack of few viscous channels,
  • a more elaborate load that evolves over time,
  • transient changes of the relative sea-level.

For this we run a glacial cycle of Antarctica with lithospheric thickness and upper-mantle viscosity from Wiens et al. (2022) and the ice thickness history from Peltier et al. (2018). We start by generating a ComputationDomain with intermediate resolution for the sake of the example and load the heterogeneous lithospheric from Pan et al. (2022) thanks to the convenience of load_dataset:

using CairoMakie, FastIsostasy

N = 140     # corresponds to 50 km resolution
Omega = ComputationDomain(3500e3, 3500e3, N, N)
(; Lon, Lat) = Omega
c = PhysicalConstants()

(_, _), Tpan, Titp = load_dataset("Lithothickness_Pan2022")
Tlitho = Titp.(Lon, Lat) .* 1e3                     # convert from m to km

function nicer_heatmap(X)
    fig = Figure(size = (800, 700))
    ax = Axis(fig[1, 1], aspect = DataAspect())
    hidedecorations!(ax)
    hm = heatmap!(ax, X)
    Colorbar(fig[1, 2], hm, height = Relative(0.6))
    return fig
end
nicer_heatmap(Tlitho)
Example block output

In a similar way, we can load the log-viscosity field from Pan et al. (2022) and plot it at about 300 km depth

(_, _, _), _, logeta_itp = load_dataset("Viscosity_Pan2022")
logeta300 = logeta_itp.(Lon, Lat, c.r_equator - 300e3)
nicer_heatmap(logeta300)
Example block output

The number of layers and the depth of viscous half-space are arbitrary parameters that have to be defined by the user. We here use a relatively shallow model (half-space begins at 300 km depth) with 1 equalisation layer and 3 intermediate layers:

mindepth = maximum(Tlitho) + 1e3
lb_vec = range(mindepth, stop = 300e3, length = 3)
lb = cat(Tlitho, [fill(lbval, Omega.Nx, Omega.Ny) for lbval in lb_vec]..., dims=3)

rlb = c.r_equator .- lb
nlb = size(rlb, 3)
lv_3D = 10 .^ cat([logeta_itp.(Lon, Lat, rlb[:, :, k]) for k in 1:nlb]..., dims=3)
size(lv_3D)
(140, 140, 4)

To prevent extreme values of the viscosity, we require it to be larger than a minimal value, fixed to be $10^{16} \, \mathrm{Pa \, s} $. We subsequently generate a LayeredEarth that embeds all the information that has been loaded so far:

eta_lowerbound = 1e16
lv_3D[lv_3D .< eta_lowerbound] .= eta_lowerbound
p = LayeredEarth(Omega, layer_boundaries = lb, layer_viscosities = lv_3D)
nicer_heatmap(log10.(p.effective_viscosity))
Example block output

We now load the ice thickness history from ICE6G, again helped by the convenience of load_dataset. We then create a vector of anomaly snapshots, between which FastIsostasy automatically interpolates linearly. To get an idea of ICE6G, the ice thickness anomaly is then visualised at LGM:

(lon, lat, t), Hice, Hitp = load_dataset("ICE6G_D")
Hice_vec = [Hitp.(Lon, Lat, tk) for tk in t]
nicer_heatmap(Hitp.(Lon, Lat, -26) - Hitp.(Lon, Lat, 0))
Example block output

Finally, we define and solve the resulting FastIsoProblem. We hereby choose the verbose=true option to track the progress of the computation.

opts = SolverOptions(verbose = true)
tyr = t .* 1e3
fip = FastIsoProblem(Omega, c, p, tyr, tyr, Hice_vec, output = "sparse", opts = opts)
solve!(fip)
Computing until t = -122400 years...
Computing until t = -122300 years...
Computing until t = -122200 years...
Computing until t = -122100 years...
Computing until t = -122000 years...
Computing until t = -120000 years...
Computing until t = -118000 years...
Computing until t = -116000 years...
Computing until t = -114000 years...
Computing until t = -112000 years...
Computing until t = -110000 years...
Computing until t = -108000 years...
Computing until t = -106000 years...
Computing until t = -104000 years...
Computing until t = -102000 years...
Computing until t = -100000 years...
Computing until t = -98000 years...
Computing until t = -96000 years...
Computing until t = -94000 years...
Computing until t = -92000 years...
Computing until t = -90000 years...
Computing until t = -88000 years...
Computing until t = -86000 years...
Computing until t = -84000 years...
Computing until t = -82000 years...
Computing until t = -80000 years...
Computing until t = -78000 years...
Computing until t = -76000 years...
Computing until t = -74000 years...
Computing until t = -72000 years...
Computing until t = -70000 years...
Computing until t = -68000 years...
Computing until t = -66000 years...
Computing until t = -64000 years...
Computing until t = -62000 years...
Computing until t = -60000 years...
Computing until t = -58000 years...
Computing until t = -56000 years...
Computing until t = -54000 years...
Computing until t = -52000 years...
Computing until t = -50000 years...
Computing until t = -48000 years...
Computing until t = -46000 years...
Computing until t = -44000 years...
Computing until t = -42000 years...
Computing until t = -40000 years...
Computing until t = -38000 years...
Computing until t = -36000 years...
Computing until t = -34000 years...
Computing until t = -32000 years...
Computing until t = -31000 years...
Computing until t = -30000 years...
Computing until t = -29000 years...
Computing until t = -28000 years...
Computing until t = -27000 years...
Computing until t = -26000 years...
Computing until t = -25000 years...
Computing until t = -24000 years...
Computing until t = -23000 years...
Computing until t = -22000 years...
Computing until t = -21000 years...
Computing until t = -20500 years...
Computing until t = -20000 years...
Computing until t = -19500 years...
Computing until t = -19000 years...
Computing until t = -18500 years...
Computing until t = -18000 years...
Computing until t = -17500 years...
Computing until t = -17000 years...
Computing until t = -16500 years...
Computing until t = -16000 years...
Computing until t = -15500 years...
Computing until t = -15000 years...
Computing until t = -14750 years...
Computing until t = -14500 years...
Computing until t = -14250 years...
Computing until t = -14000 years...
Computing until t = -13750 years...
Computing until t = -13500 years...
Computing until t = -13250 years...
Computing until t = -13000 years...
Computing until t = -12750 years...
Computing until t = -12500 years...
Computing until t = -12250 years...
Computing until t = -12000 years...
Computing until t = -11750 years...
Computing until t = -11500 years...
Computing until t = -11250 years...
Computing until t = -11000 years...
Computing until t = -10750 years...
Computing until t = -10500 years...
Computing until t = -10250 years...
Computing until t = -10000 years...
Computing until t = -9750 years...
Computing until t = -9500 years...
Computing until t = -9250 years...
Computing until t = -9000 years...
Computing until t = -8750 years...
Computing until t = -8500 years...
Computing until t = -8250 years...
Computing until t = -8000 years...
Computing until t = -7500 years...
Computing until t = -7000 years...
Computing until t = -6500 years...
Computing until t = -6000 years...
Computing until t = -5500 years...
Computing until t = -5000 years...
Computing until t = -4500 years...
Computing until t = -4250 years...
Computing until t = -4000 years...
Computing until t = -3750 years...
Computing until t = -3500 years...
Computing until t = -3250 years...
Computing until t = -3000 years...
Computing until t = -2750 years...
Computing until t = -2500 years...
Computing until t = -2250 years...
Computing until t = -2000 years...
Computing until t = -1750 years...
Computing until t = -1500 years...
Computing until t = -1250 years...
Computing until t = -1000 years...
Computing until t = -750 years...
Computing until t = -500 years...
Computing until t = -250 years...
Computing until t = 0 years...
Computing until t = 250 years...

For a resolution of 50 km, the computation time of this last step is less than a minute on a modern i7 (Intel i7-10750H CPU @ 2.60GHz)! We visualise three snapshots of displacements that roughly correspond to LGM, the end of the meltwater pulse 1A and the present-day:

tplot = [-26, -12, 0]
fig = Figure(size = (1200, 400))
opts = ( colormap = :PuOr, colorrange = (-400, 400) )
for k in eachindex(tplot)
    kfi = argmin( abs.(tplot[k] * 1e3 .- tyr) )
    ax = Axis(fig[1, k], aspect = DataAspect(), title = "t = $(t[kfi]) kyr")
    hidedecorations!(ax)
    heatmap!(ax, fip.out.u[kfi] + fip.out.ue[kfi]; opts...)
    println(kfi)
end
Colorbar(fig[1, 4], height = Relative(0.6); opts...)
fig
Example block output

The displayed fields are displacement anomalies w.r.t. to the last interglacial, defined as the reference for the ice thickness anomalies. In Swierczek-Jereczek et al. (2024) these computations are performed on a finer grid, with an interactive sea level, and show great agreement with a 3D GIA model that runs between 10,000-100,000 slower (however at with advantage of obtaining a global and richer output).