Parameter inversion

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 (Huang et al., 2021; Huang and Huang, Feb 2021), 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(Wiens et al., 2022) can be loaded by using load_dataset with appropriate depths of the layer boundaries:

using CairoMakie
using FastIsostasy

Omega = ComputationDomain(3000e3, 5)
c = PhysicalConstants()
lb = [100e3, 200e3, 300e3]
_, eta, eta_itp = load_dataset("Wiens2022")
loglv = cat([eta_itp.(Omega.X, Omega.Y, z) for z in lb]..., dims = 3)
lv = 10 .^ loglv
p = LayeredEarth(Omega, layer_boundaries = lb, layer_viscosities = lv)
heatmap(p.effective_viscosity)

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
Hcylinder = uniform_ice_cylinder(Omega, R, H, center = [-1000e3, -1000e3])
Hice = [zeros(Omega.Nx, Omega.Ny), Hcylinder, Hcylinder]
t_out = years2seconds.(collect(1e3:1e3:2e3))
pushfirst!(t_out, t_out[1]-1e-8)
t_Hice = copy(t_out)

true_viscosity = copy(p.effective_viscosity)
fip = FastIsoProblem(Omega, c, p, t_out, t_Hice, Hice)
solve!(fip)

┌ Warning: Duplicated knots were deduplicated. Use Interpolations.deduplicate_knots!(knots) explicitly to avoid this warning.
│   k1 =
│    3-element Vector{Float64}:
│     3.15576e10
│     3.15576e10
│     6.31152e10
└ @ Interpolations ~/.julia/packages/Interpolations/91PhN/src/gridded/gridded.jl:77

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[2:end]),
    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 = (19.5, 21.5)
    fig = Figure(size = (1800, 1000), fontsize = 40)
    axs = [Axis(fig[1,i], aspect = DataAspect()) for i in 1:2]
    [hidedecorations!(ax) for ax in axs]
    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
fig1 = plot_viscfields(paraminv)

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

solve!(paraminv)
fig2 = 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.