Page Navigation

Overview

This program is built on the fast nonlinear gravity inversion method proposed by Uieda & Barbosa (2017) in spherical coordinates, but with the forward modelling engine replaced by the modern tesseroid implementation from Fatiando a Terra. This preserves the original clear inversion framework while significantly improving computation speed and code maintainability.

At its core, the inversion solves a regularized nonlinear least-squares problem: given observed gravity $g_{\mathrm{obs}}$, find Moho undulations $x$ such that the forward-modelled gravity $g(x)$ best fits the observations while maintaining spatial smoothness.

Mathematical Background

1. Objective Function

The inversion treats Moho undulations $x \in \mathbb{R}^{N}$ ($N = n_{\mathrm{lat}} \times n_{\mathrm{lon}}$) as unknowns. The objective function uses classic Tikhonov regularization:

where $g(x)$ is the tesseroid forward operator, $x_0$ is the initial model, $L$ is the first-order finite-difference operator, and $\mu \ge 0$ is the regularization parameter. When $\mu = 0$, this reduces to pure data fitting.

2. First-Order Difference Operator $L$

$L$ encodes the difference between each pair of adjacent grid nodes in one matrix row, covering both longitudinal and latitudinal directions. For node $(i, j)$, define flattened index $k(i, j) = i \cdot n_{\mathrm{lon}} + j$, then

The final $L$ is a sparse matrix of shape $(n_{\mathrm{edges}},\, N)$ formed by stacking both blocks. The regularization term $\lVert L(x - x_0) \rVert^2$ acts as a discrete "gradient energy," penalizing large jumps between neighbouring cells.

3. Bott Diagonal Jacobian

Computing the full Jacobian $\partial g_i / \partial x_j$ requires partial derivatives for every voxel and is prohibitively expensive. Bott (1960) approximated this by noting that each observation point's gravity is dominated by the voxel directly below it, yielding a diagonal Jacobian:

Taking the absolute value of the density contrast is critical — it ensures $J^\top J$ is always a positive diagonal matrix, guaranteeing the linear system is positive definite. The factor $10^{5}$ converts m/s² to mGal.

4. Gauss-Newton Iteration

A second-order approximation of the objective gives the iteration formula:

where the residual $r = g(x) - g_{\mathrm{obs}}$. Each update uses Armijo backtracking line search to determine step size $\alpha$:

Defaults: $c_1 = 10^{-4}$, initial $\alpha = 1$, halved on each failure, up to 30 attempts. After each iteration, a safety clipping $x \le h - r_{\mathrm{ref}}$ ensures the Moho stays below the observation surface.

5. Convergence Criteria

Iteration stops when any of the following is satisfied:

Geometry & Density Conventions

The geometric conventions are central to inversion correctness — a sign error will drive convergence in the wrong direction. The figure below shows the spatial relationships among the observation surface, reference Moho, and actual Moho:

Density Convention

The density_contrast parameter is positive, meaning the material below the Moho is denser than above (mantle > crust). The sign of the density array fed to the forward engine follows these rules:

Inversion Flowchart

The diagram below shows the complete workflow from input data to final results. Blue boxes are user interfaces; grey boxes are internal steps:

Workflow

• Prepare regular grid lon, lat, and corresponding 2D field data.
• Set observation height height, reference Moho reference, initial undulation surface, and density_contrast.
• Create a MohoLayer object; run one forward computation to check magnitude and sign.
• Pass observed gravity go to MohoInversion, set mu and max_iter, then run the inversion.
• Evaluate quality via residuals, convergence curves, and test-point results.

API Quick Reference

MohoLayer

MohoInversion

Evaluation & Seismic Constraints

Independent Test-Point Evaluation

If you held out some observation points from the inversion, use evaluate_testset to check extrapolation quality:

go_pred, residual_test, scores = inv.evaluate_testset(
    lon_pts=test_lon,
    lat_pts=test_lat,
    go_test=test_go,
)
# scores = {'rmse': ..., 'mae': ..., 'bias': ..., 'r2': ...}

Seismic Crustal Thickness Evaluation

The evaluate_seismic method interpolates the inverted Moho undulation to seismic station locations and compares with seismically determined crustal thickness. The geometric relation is:

where topo is the topographic height (relative to mean sphere radius), reference is typically negative (Moho below mean sphere), and $s$ is the inverted undulation. The returned diff = thickness_inv − seismic_thickness; positive values indicate the inversion overestimates thickness.

thickness_inv, diff, scores = inv.evaluate_seismic(
    lon_pts=seismic_lon,
    lat_pts=seismic_lat,
    seismic_thickness=seismic_h,   # in metres
    topo_pts=topo_at_seismic,      # in metres
)
# Prints MAE / RMSE / Bias / R² (in km)

Tips & Troubleshooting

• Start with a closed-loop synthetic test (forward → add noise → invert) to verify signs, units, and density contrast direction.
• If results oscillate or show excessive high-frequency undulations, increase $\mu$ first before modifying the initial model.
• For large grids, do a low-resolution trial run (e.g., halve the grid size) before scaling up.
• When applying to different bodies, check mean radius, reference surface definition, and density contrast meaning.
• Don't rely on RMS alone — combine spatial residual maps, independent test points, and seismic constraints.
• Persistent line search failures usually mean the initial model is too far from the true solution, or $\mu$ is poorly chosen — check whether the initial forward magnitude matches observations.
• If the residual histogram deviates significantly from Gaussian, there may be systematic bias or edge effects; consider detrending before inversion.