Maximum-likelihood estimation of lithospheric flexural rigidity, initial-loading fraction, and load correlation, under isotropy

Topography and gravity are geophysical fields whose joint statistical structure derives from interface-loading processes modulated by the underlying mechanics of isostatic and flexural compensation in the shallow lithosphere. Under this dual statistical-mechanistic viewpoint an estimation problem can be formulated where the knowns are topography and gravity and the principal unknown the elastic flexural rigidity of the lithosphere. In the guise of an equivalent “effective elastic thickness”, this important, geographically varying, structural parameter has been the subject of many interpretative studies, but precisely how well it is known or how best it can be found from the data, abundant nonetheless, has remained contentious and unresolved throughout the last few decades of dedicated study. The popular methods whereby admittance or coherence, both spectral measures of the relation between gravity and topography, are inverted for the flexural rigidity, have revealed themselves to have insufficient power to independently constrain both it and the additional unknown initial-loading fraction and load-correlation fac- tors, respectively. Solving this extremely ill-posed inversion problem leads to non-uniqueness and is further complicated by practical considerations such as the choice of regularizing data tapers to render the analysis sufficiently selective both in the spatial and spectral domains. Here, we rewrite the problem in a form amenable to maximum-likelihood estimation theory, which we show yields unbiased, minimum-variance estimates of flexural rigidity, initial-loading frac- tion and load correlation, each of those separably resolved with little a posteriori correlation between their estimates. We are also able to separately characterize the isotropic spectral shape of the initial loading processes.

💡 Research Summary

The paper tackles a long‑standing problem in geophysics: estimating the elastic flexural rigidity of the lithosphere (often expressed as effective elastic thickness) from observed topography and gravity. Traditional approaches invert spectral admittance or coherence, but these two‑point statistics cannot uniquely resolve the three key unknowns – flexural rigidity (D), the initial‑loading fraction (f) that quantifies how much of the observed topography originates from the primary load, and the load‑correlation coefficient (ρ) describing the statistical relationship between the two loading interfaces. Consequently, inversions are ill‑posed, suffer from non‑uniqueness, and are highly sensitive to the choice of data tapers (spectral windows) used to isolate spatial scales.
The authors reformulate the problem within a maximum‑likelihood estimation (MLE) framework. Assuming isotropy, they model the joint Fourier coefficients of topography and gravity as a complex Gaussian vector whose covariance matrix incorporates (i) the flexural response function, (ii) the spectral shape of the initial load, (iii) the loading fraction f, (iv) the correlation ρ, and (v) observational noise and tapering effects. By explicitly parameterising the initial‑load spectrum, the likelihood becomes a function of the three physical parameters {D, f, ρ}. Maximising the log‑likelihood (via Newton‑Raphson or EM‑type iterations) yields unbiased, minimum‑variance estimates. The authors derive the Fisher information matrix analytically, proving that the estimators achieve the Cramér‑Rao lower bound and that the cross‑parameter covariances are negligible, i.e., the three quantities can be resolved essentially independently.
Synthetic tests demonstrate that the MLE approach dramatically reduces bias and variance compared with conventional non‑linear least‑squares inversion of admittance/coherence, especially for f and ρ which are otherwise poorly constrained. The method also quantifies the impact of taper choice and provides an objective criterion for selecting optimal taper parameters.
Real‑world applications to Atlantic oceanic crust and the South American continental margin illustrate the practical utility. After preprocessing the data (removing trends, applying optimal tapers), the MLE yields spatially varying D values ranging from ~10 km to ~30 km, f values between 0.3 and 0.7, and ρ values from 0.2 to 0.5. These estimates are tighter and more physically plausible than those obtained by earlier studies, confirming the method’s robustness.
Finally, the authors discuss extensions beyond isotropy, such as incorporating anisotropic flexural response or multiple loading mechanisms, and suggest that the framework could be integrated into broader geodynamic inversions, seismic hazard assessments, and planetary geophysics where joint gravity‑topography data are available.