Glacial cycle
The previous examples focused on benchmarking FastIsostasy against analytical, numerical 1D and numerical 3D solutions. However, these cases were largely idealised. To test the model on more realistic simulations, we now provide compute the GIA response in Antarctica over the last glacial cycle. This includes the use of:
- a heterogeneous lithosphere thickness (Pan et al., 2022),
- a heterogeneous upper-mantle viscosity (Pan et al., 2022),
- a stack of few viscous channels,
- a more elaborate load that evolves over time (Peltier et al., 2018),
- transient changes of the relative sea-level.
We start by generating a RegionalDomain with intermediate resolution for the sake of the example and load the ice history thanks to the convenience of load_dataset. To get an idea of the ICE6G forcing, the ice thickness is visualised at the last glacial maximum (LGM):
using CairoMakie, FastIsostasy
N = 140
domain = RegionalDomain(3500e3, 3500e3, N, N) # resolution = 50 km
Lon, Lat = domain.Lon, domain.Lat
(_, _, t), Hice, Hitp = load_dataset("ICE6G_D")
Hice_vec = [Hitp.(Lon, Lat, tk) for tk in t]
k_lgm = argmax([mean(Hice_vec[k]) for k in eachindex(Hice_vec)])
plot_load(domain, Hice_vec[k_lgm])This already looks like a much more exciting ice thickness field! Here again, the ice history is wrapped into an interpolator, which is passed to an instance of BoundaryConditions. We define the SeaLevel to include the gravitational response by making the surface a LaterallyVariableSeaSurface. Furthermore, we allow the changes in sea level to affect the deformational response of the solid Earth by setting InteractiveSealevelLoad. Finally, we compute the evolution of the barystatic sea level (BSL) according to a piece-wise constant approximation of the ocean surface as a function of the BSL:
it = TimeInterpolatedIceThickness(t .* 1e3, Hice_vec, domain)
bcs = BoundaryConditions(domain, ice_thickness = it)
sealevel = SeaLevel(
surface = LaterallyVariableSeaSurface(),
load = InteractiveSealevelLoad(),
bsl = PiecewiseConstantBSL(),
)SeaLevel{LaterallyVariableSeaSurface, InteractiveSealevelLoad, PiecewiseConstantBSL{Float32, ReferenceBSL{Float32, FastIsostasy.TimeInterpolation0D{Float32}}}, InternalUpdateBSL}(LaterallyVariableSeaSurface(), InteractiveSealevelLoad(), PiecewiseConstantBSL{Float32, ReferenceBSL{Float32, FastIsostasy.TimeInterpolation0D{Float32}}}(ReferenceBSL{Float32, FastIsostasy.TimeInterpolation0D{Float32}}(0.0f0, 3.625f14, Float32[-150.0, -149.9, -149.8, -149.7, -149.6, -149.5, -149.4, -149.3, -149.2, -149.1 … 69.1, 69.2, 69.3, 69.4, 69.5, 69.6, 69.7, 69.8, 69.9, 70.0], Float32[3.4401262f14, 3.4402312f14, 3.4403148f14, 3.440368f14, 3.4404507f14, 3.4404983f14, 3.4405882f14, 3.440671f14, 3.4407245f14, 3.4408074f14 … 3.852362f14, 3.8525696f14, 3.852777f14, 3.8529833f14, 3.8531923f14, 3.8534044f14, 3.853611f14, 3.853818f14, 3.854026f14, 3.8542305f14], FastIsostasy.TimeInterpolation0D{Float32}(Float32[-150.0, -149.9, -149.8, -149.7, -149.6, -149.5, -149.4, -149.3, -149.2, -149.1 … 69.1, 69.2, 69.3, 69.4, 69.5, 69.6, 69.7, 69.8, 69.9, 70.0], Float32[3.3774267f14, 3.3775297f14, 3.377612f14, 3.3776643f14, 3.377745f14, 3.377792f14, 3.3778804f14, 3.3779616f14, 3.3780143f14, 3.3780955f14 … 3.782149f14, 3.782353f14, 3.7825566f14, 3.7827592f14, 3.7829642f14, 3.7831726f14, 3.7833756f14, 3.7835786f14, 3.783783f14, 3.7839836f14], false)), 0.0f0, 3.625f14), InternalUpdateBSL())Finally, we load the interpolators of earth structure thanks to the convenience function load_dataset.
(_, _), Tpan, Titp = load_dataset("Lithothickness_Pan2022")
Tlitho = Titp.(Lon, Lat) .* 1e3 # convert from m to km
(_, _, _), _, logeta_itp = load_dataset("Viscosity_Pan2022")(([-179.5, -179.0, -178.5, -178.0, -177.5, -177.0, -176.5, -176.0, -175.5, -175.0 … 175.5, 176.0, 176.5, 177.0, 177.5, 178.0, 178.5, 179.0, 179.5, 180.0], [-90.0, -89.5, -89.0, -88.5, -88.0, -87.5, -87.0, -86.5, -86.0, -85.5 … 85.5, 86.0, 86.5, 87.0, 87.5, 88.0, 88.5, 89.0, 89.5, 90.0], [3.5298448275852134e6, 3.5796896551732332e6, 3.6300000000001593e6, 3.680600706716196e6, 3.731678445228479e6, 3.7827561837463e6, 3.833833922264341e6, 3.884911660773705e6, 3.9359893992920252e6, 3.9870671378111225e6 … 6.151000000000001e6, 6.171000000002125e6, 6.190999999973099e6, 6.2109999999932265e6, 6.231000000027917e6, 6.250999999971964e6, 6.274999999972597e6, 6.298499999991569e6, 6.322000000031837e6, 6.34659999997364e6]), Float32[21.898117 21.812355 … 22.85709 22.764662; 21.898117 21.812355 … 22.854446 22.764662; … ; 21.898117 21.812355 … 22.857298 22.764662; 21.898117 21.812355 … 22.857801 22.764662;;; 22.37607 22.264135 … 23.260921 23.19398; 22.37607 22.264135 … 23.254314 23.19398; … ; 22.37607 22.264135 … 23.269352 23.19398; 22.37607 22.264135 … 23.26836 23.19398;;; 22.761507 22.656908 … 23.116743 23.120834; 22.761507 22.656908 … 23.116238 23.120834; … ; 22.761507 22.656908 … 23.11527 23.120834; 22.761507 22.656908 … 23.115332 23.120834;;; … ;;; 22.42295 22.171955 … 20.865768 20.805326; 22.42295 22.166954 … 20.866772 20.805326; … ; 22.42295 22.181955 … 20.86479 20.805326; 22.42295 22.176954 … 20.864803 20.805326;;; 22.422972 22.17197 … 20.872458 20.811415; 22.422972 22.166971 … 20.873459 20.811415; … ; 22.422972 22.18197 … 20.871456 20.811415; 22.422972 22.176971 … 20.871456 20.811415;;; 22.42297 22.17197 … 20.87197 20.81097; 22.42297 22.16697 … 20.87297 20.81097; … ; 22.42297 22.18197 … 20.87097 20.81097; 22.42297 22.17697 … 20.87097 20.81097], 721×361×64 extrapolate(interpolate((::Vector{Float64},::Vector{Float64},::Vector{Float64}), ::Array{Float32, 3}, Gridded(Linear())), Flat()) with element type Float32:
[:, :, 1] =
21.8981 21.8124 21.7422 21.6872 … 22.893 22.9051 22.8571 22.7647
21.8981 21.8124 21.7422 21.6872 22.8855 22.9003 22.8544 22.7647
21.8981 21.8124 21.7422 21.6872 22.8833 22.8948 22.8508 22.7647
21.8981 21.8124 21.7422 21.6872 22.8818 22.8899 22.8477 22.7647
21.8981 21.8124 21.7422 21.6872 22.8769 22.8893 22.848 22.7647
21.8981 21.8124 21.7422 21.6872 … 22.8724 22.8899 22.849 22.7647
21.8981 21.8124 21.7422 21.6875 22.8696 22.8907 22.8495 22.7647
21.8981 21.8124 21.7422 21.6875 22.8668 22.8899 22.8493 22.7647
21.8981 21.8123 21.7422 21.6874 22.8572 22.8842 22.8475 22.7647
21.8981 21.8122 21.7422 21.6873 22.8488 22.8774 22.8435 22.7647
⋮ ⋱ ⋮
21.8981 21.8115 21.7422 21.6896 22.9433 22.9287 22.8659 22.7647
21.8981 21.8118 21.7422 21.6878 22.9359 22.9241 22.864 22.7647
21.8981 21.8118 21.7422 21.6861 22.9279 22.9207 22.8628 22.7647
21.8981 21.8118 21.7422 21.6861 … 22.9258 22.9212 22.8632 22.7647
21.8981 21.8118 21.7422 21.6872 22.923 22.9207 22.8627 22.7647
21.8981 21.8121 21.7422 21.6872 22.9148 22.9145 22.8598 22.7647
21.8981 21.8124 21.7422 21.6872 22.9064 22.9083 22.857 22.7647
21.8981 21.8124 21.7422 21.6872 22.903 22.9081 22.8573 22.7647
21.8981 21.8124 21.7422 21.6872 … 22.9007 22.9083 22.8578 22.7647
[:, :, 2] =
22.3761 22.2641 22.1733 22.094 … 23.2252 23.2759 23.2609 23.194
22.3761 22.2641 22.1733 22.094 23.2104 23.2617 23.2543 23.194
22.3761 22.2641 22.1733 22.094 23.2072 23.259 23.2526 23.194
22.3761 22.2641 22.1733 22.094 23.2069 23.2582 23.2527 23.194
22.3761 22.2641 22.1733 22.094 23.2005 23.258 23.2535 23.194
22.3761 22.2641 22.1733 22.094 … 23.1946 23.2582 23.2545 23.194
22.3761 22.2641 22.1733 22.0949 23.1957 23.2582 23.2535 23.194
22.3761 22.2641 22.1733 22.0949 23.1969 23.2582 23.2525 23.194
22.3761 22.2643 22.1733 22.0951 23.192 23.2606 23.2546 23.194
22.3761 22.2644 22.1733 22.0953 23.1849 23.2592 23.2549 23.194
⋮ ⋱ ⋮
22.3761 22.2613 22.1733 22.104 23.2769 23.3149 23.2789 23.194
22.3761 22.2623 22.1733 22.0986 23.2634 23.3019 23.2724 23.194
22.3761 22.2623 22.1733 22.0933 23.2498 23.2878 23.2658 23.194
22.3761 22.2623 22.1733 22.0932 … 23.2479 23.2868 23.265 23.194
22.3761 22.2623 22.1733 22.094 23.2488 23.2878 23.2661 23.194
22.3761 22.2632 22.1733 22.094 23.2428 23.2883 23.2673 23.194
22.3761 22.2641 22.1733 22.094 23.2372 23.2888 23.2684 23.194
22.3761 22.2641 22.1733 22.094 23.2393 23.2907 23.2694 23.194
22.3761 22.2641 22.1733 22.094 … 23.2395 23.2888 23.2684 23.194
[:, :, 3] =
22.7615 22.6569 22.5638 22.4852 … 22.9725 23.0703 23.1167 23.1208
22.7615 22.6569 22.5638 22.4852 22.9725 23.0693 23.1162 23.1208
22.7615 22.6569 22.5638 22.4852 22.9559 23.0545 23.1089 23.1208
22.7615 22.6569 22.5638 22.4852 22.94 23.0396 23.1015 23.1208
22.7615 22.6569 22.5638 22.4852 22.9396 23.0377 23.1005 23.1208
22.7615 22.6569 22.5638 22.4852 … 22.943 23.0396 23.1014 23.1208
22.7615 22.6569 22.5638 22.4852 22.943 23.0396 23.1016 23.1208
22.7615 22.6569 22.5638 22.4852 22.9429 23.0396 23.1018 23.1208
22.7615 22.6578 22.5638 22.4861 22.9441 23.0391 23.1014 23.1208
22.7615 22.6588 22.5638 22.487 22.9435 23.0393 23.1013 23.1208
⋮ ⋱ ⋮
22.7615 22.6571 22.5638 22.4845 22.9805 23.0658 23.1125 23.1208
22.7615 22.657 22.5638 22.4847 22.9833 23.0667 23.1129 23.1208
22.7615 22.657 22.5638 22.485 22.9842 23.0676 23.1133 23.1208
22.7615 22.657 22.5638 22.4851 … 22.9771 23.0676 23.1143 23.1208
22.7615 22.657 22.5638 22.4852 22.9686 23.0676 23.1152 23.1208
22.7615 22.657 22.5638 22.4852 22.9679 23.0676 23.1152 23.1208
22.7615 22.6569 22.5638 22.4852 22.9689 23.0676 23.1153 23.1208
22.7615 22.6569 22.5638 22.4852 22.9686 23.0675 23.1153 23.1208
22.7615 22.6569 22.5638 22.4852 … 22.9687 23.0676 23.1153 23.1208
;;; …
[:, :, 62] =
22.423 22.172 21.913 22.3061 … 21.1364 20.9679 20.8658 20.8053
22.423 22.167 21.904 22.2871 21.1413 20.9699 20.8668 20.8053
22.423 22.162 21.893 22.2681 21.1463 20.9719 20.8678 20.8053
22.423 22.156 21.884 22.2491 21.1503 20.9739 20.8688 20.8053
22.423 22.151 21.874 22.2301 21.1552 20.9768 20.8698 20.8053
22.423 22.146 21.864 22.2111 … 21.1592 20.9798 20.8698 20.8053
22.423 22.141 22.6312 22.1931 21.1632 20.9808 20.8707 20.8053
22.423 22.136 22.6182 22.1751 21.1672 20.9828 20.8717 20.8053
22.423 22.131 22.6052 22.1571 21.1721 20.9848 20.8728 20.8053
22.423 22.126 22.5912 22.1391 21.1761 20.9868 20.8737 20.8053
⋮ ⋱ ⋮
22.423 22.22 22.004 22.4841 21.0927 20.947 20.8598 20.8053
22.423 22.214 21.994 22.4631 21.0986 20.9499 20.8608 20.8053
22.423 22.209 21.984 22.4431 21.1036 20.9519 20.8618 20.8053
22.423 22.203 21.974 22.4241 … 21.1086 20.9549 20.8618 20.8053
22.423 22.198 21.963 22.4031 21.1125 20.9569 20.8628 20.8053
22.423 22.193 21.953 22.3831 21.1185 20.9589 20.8638 20.8053
22.423 22.187 21.944 22.3641 21.1235 20.9609 20.8648 20.8053
22.423 22.182 21.933 22.3451 21.1284 20.9629 20.8648 20.8053
22.423 22.177 21.923 22.3251 … 21.1334 20.9659 20.8648 20.8053
[:, :, 63] =
22.423 22.172 21.913 22.306 … 21.1456 20.9755 20.8725 20.8114
22.423 22.167 21.904 22.287 21.1507 20.9775 20.8735 20.8114
22.423 22.162 21.893 22.268 21.1557 20.9795 20.8745 20.8114
22.423 22.156 21.884 22.249 21.1597 20.9815 20.8755 20.8114
22.423 22.151 21.874 22.23 21.1647 20.9845 20.8765 20.8114
22.423 22.146 21.864 22.211 … 21.1687 20.9875 20.8765 20.8114
22.423 22.141 22.631 22.193 21.1727 20.9885 20.8775 20.8114
22.423 22.136 22.618 22.175 21.1767 20.9905 20.8785 20.8114
22.423 22.131 22.605 22.157 21.1817 20.9925 20.8795 20.8114
22.423 22.126 22.591 22.139 21.1857 20.9945 20.8805 20.8114
⋮ ⋱ ⋮
22.423 22.22 22.004 22.484 21.1016 20.9545 20.8665 20.8114
22.423 22.214 21.994 22.463 21.1076 20.9575 20.8675 20.8114
22.423 22.209 21.984 22.443 21.1126 20.9595 20.8685 20.8114
22.423 22.203 21.974 22.424 … 21.1176 20.9625 20.8685 20.8114
22.423 22.198 21.963 22.403 21.1216 20.9645 20.8695 20.8114
22.423 22.193 21.953 22.383 21.1276 20.9665 20.8705 20.8114
22.423 22.187 21.944 22.364 21.1326 20.9685 20.8715 20.8114
22.423 22.182 21.933 22.345 21.1376 20.9705 20.8715 20.8114
22.423 22.177 21.923 22.325 … 21.1426 20.9735 20.8715 20.8114
[:, :, 64] =
22.423 22.172 21.913 22.306 22.513 … 21.145 20.975 20.872 20.811
22.423 22.167 21.904 22.287 22.482 21.15 20.977 20.873 20.811
22.423 22.162 21.893 22.268 22.453 21.155 20.979 20.874 20.811
22.423 22.156 21.884 22.249 22.424 21.159 20.981 20.875 20.811
22.423 22.151 21.874 22.23 22.395 21.164 20.984 20.876 20.811
22.423 22.146 21.864 22.211 22.366 … 21.168 20.987 20.876 20.811
22.423 22.141 22.631 22.193 22.337 21.172 20.988 20.877 20.811
22.423 22.136 22.618 22.175 22.31 21.176 20.99 20.878 20.811
22.423 22.131 22.605 22.157 22.282 21.181 20.992 20.879 20.811
22.423 22.126 22.591 22.139 22.254 21.185 20.994 20.88 20.811
⋮ ⋱ ⋮
22.423 22.22 22.004 22.484 22.141 21.101 20.954 20.866 20.811
22.423 22.214 21.994 22.463 22.116 21.107 20.957 20.867 20.811
22.423 22.209 21.984 22.443 22.093 21.112 20.959 20.868 20.811
22.423 22.203 21.974 22.424 22.068 … 21.117 20.962 20.868 20.811
22.423 22.198 21.963 22.403 22.046 21.121 20.964 20.869 20.811
22.423 22.193 21.953 22.383 22.021 21.127 20.966 20.87 20.811
22.423 22.187 21.944 22.364 21.999 21.132 20.968 20.871 20.811
22.423 22.182 21.933 22.345 22.573 21.137 20.97 20.871 20.811
22.423 22.177 21.923 22.325 22.543 … 21.142 20.973 20.871 20.811)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
layerboundary_vec = range(mindepth, stop = 300e3, length = 3)
lb = cat(Tlitho, [fill(lbval, domain.nx, domain.ny)
for lbval in layerboundary_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)
eta_lowerbound = 1e16
lv_3D[lv_3D .< eta_lowerbound] .= eta_lowerbound
maskactive = gaussian_smooth(Hice_vec[k_lgm], domain, 0.05, 0) .> 10
solidearth = SolidEarth(
domain,
layer_boundaries = lb,
layer_viscosities = lv_3D,
maskactive = maskactive, # required when using `InteractiveSealevelLoad`
)
fig = plot_earth(domain, solidearth)This already looks like a much more exciting solid Earth structure!
Finally, we define the output struct, the Simulation and run! it.
nout = NativeOutput(vars = [:u, :ue, :dz_ss, :z_ss, :H_ice], t = [-26f3, -12f3, 0])
sim = Simulation(domain, bcs, sealevel, solidearth; nout = nout)
run!(sim)
println("Computation time: ", sim.nout.computation_time)Ok, that was fast! Let's visualise three snapshots of displacements that roughly correspond to LGM, the end of the meltwater pulse 1A and the present-day:
copts = (colormap = :PuOr, colorrange = (-500, 500))
fig = plot_out_over_time(sim, :u_tot, [-26f3, -12f3, 0], copts)This looks very much like what is obtained by Seakon (Swierczek-Jereczek et al., 2024, Fig.9.g), a 3D GIA model.