Dr Sia Ghelichkhan

I work on numerical methods for geophysical continua and on the inverse problems that let us learn about the Earth's interior from observations made at its surface. My research covers mantle convection and its reconstruction through deep time, the viscoelastic response of the solid Earth to ice sheet loading and its consequences for sea level, and variably saturated groundwater flow at continental scale. The unifying method is adjoint-based optimisation applied to large, non-linear, time-dependent forward models.

Research overview

Most of what we want to know about the Earth's interior is not directly observable. We can image the mantle seismically as it is today, but not as it was fifty million years ago. We can measure how the ground surface moves, but not the viscosity structure that governs that motion. We can weigh a continent's water from orbit, but not see where within the subsurface that water sits. In each case the observation is an indirect, incomplete and noisy functional of the thing we care about, and recovering one from the other is an inverse problem.

Posing these problems as PDE-constrained optimisation makes them tractable, provided the gradient of the misfit with respect to the model's unknowns can be computed. That is what the adjoint gives, and it is the thread running through everything below.

G-ADOPT and adjoint methods

I am a co-lead of G-ADOPT, the Geoscientific ADjoint Optimisation PlaTform, a computational platform for simulating geoscientific flows and solving the inverse problems that accompany them. It is funded by the ARDC, AuScope and the ARC, and developed with Imperial College London and the University of Sydney.

The technical premise is that adjoints should be generated rather than written. G-ADOPT is built on Firedrake, an automated finite element system, so a model is expressed in Python in a form close to its variational statement, and the compiled forward code and the discrete adjoint both follow from that description through code generation and algorithmic differentiation. The consequence is practical: the physics can be changed, a new rheology introduced or a new observable added, without re-deriving and re-implementing an adjoint by hand. The foundations are set out in Davies, Kramer, Ghelichkhan and Gibson (Geoscientific Model Development, 2022), and the inversion machinery in Ghelichkhan et al. (Geoscientific Model Development, 2024).

The platform spans incompressible Boussinesq convection, compressible anelastic formulations and viscoelastic Maxwell rheology, in two and three dimensions including full spherical shells, with implicit and explicit time stepping through Irksome, multi-material level-set methods, and PDE-constrained optimisation through the Rapid Optimisation Library. It runs at scale on NCI's Gadi. Around it sits the supporting infrastructure: G-DRIFT for curated geodynamic and seismological datasets, gtrack for plate reconstruction boundary conditions, and SRTS and LLNL-ToFi for applying the resolution operators of published tomography models to model predictions, so that comparisons are made on equal terms.

Mantle convection

Seismic tomography constrains the present-day state of the mantle. It says nothing directly about how that state arose, and forward models started from an assumed initial condition can only test one hypothesis at a time. The inverse formulation asks instead which initial condition, tens of millions of years in the past, evolves into the mantle we observe now, subject to the geological and plate-kinematic record along the way.

I derived the adjoint equations for compressible and thermochemical mantle convection and verified them through twin experiments (Ghelichkhan and Bunge, Proceedings of the Royal Society A, 2018), and applied the method to global retrodictions of early Cenozoic mantle flow, tracking the evolution of dynamic topography, deep mantle structure and sublithospheric stress (Ghelichkhan, Bunge and Oeser, Geophysical Journal International, 2021). Current work relaxes the assumptions that made those early inversions tractable: non-linear rheology, compositional heterogeneity through level-set methods, and additional observables including intra-plate lava geochemistry.

A three-dimensional adjoint reconstruction of mantle flow. A recovered initial condition is evolved forward through the model's time window and compared with present-day mantle structure from seismic tomography.

Glacial isostatic adjustment and sea level

The redistribution of ice and water deforms the solid Earth, perturbs the gravity field and the rotation vector, and produces a sea-level response that is strongly non-uniform. Inverting an observed sea-level or geodetic signal for its causes is under-determined in a specific way: the signal depends jointly on the ice loading history and on the viscosity structure of the mantle, and conventional practice explores that space one candidate Earth model at a time.

Building the adjoint of the viscoelastic problem removes that restriction and allows joint inference of both. This capability now exists in G-ADOPT (Scott et al., Geoscientific Model Development, 2026) and underpins the CoastRI work with AuScope and Geoscience Australia. My ARC DECRA project applies it to the mid-Pliocene warm period, where palaeo-shoreline elevations carry a dynamic topography signal accumulated over three million years that must be separated from the glacial isostatic and eustatic contributions before the record can be used as an analogue for future warming.

Viscoelastic deformation of the solid Earth under a changing ice load, simulated in G-ADOPT. This is the work of Will Scott, and is described in Scott et al. (Geoscientific Model Development, 2026).

Continental groundwater

Groundwater is the largest accessible freshwater store on the Australian continent and the most sparsely observed. Variably saturated flow through porous media, described by the Richards equation, is strongly non-linear and degenerate; the hydraulic properties that control it are largely unknown at continental scale; and the available observations, satellite gravity, altimetry, InSAR and a sparse bore network, constrain integrals of the state rather than the state itself. GRACE weighs a region without resolving the vertical distribution of the mass change. InSAR sees the surface subside without identifying which unit compacted.

This is the same class of problem as the two above, and I am developing adjoint capability for the Richards equation within G-ADOPT so that continental groundwater models can assimilate these observations directly, and so that the resolution of the inversion can be assessed rather than assumed. The work is supported by the SmartSat CRC and the ANU Institute for Water Futures, and runs in collaboration with the Environmental Geodesy group.

Moisture content in a three-dimensional soil column during pumping. The cone of depression develops as water is withdrawn, and the unsaturated zone above the water table responds non-linearly.

The same physics in cross-section. Colours show moisture content and the lines are flow paths as the water table draws down.

Group

Postdoctoral researchers: Dr Liam Morrow; Dr Hamish Brown, with Rhodri Davies; Dr Will Scott, with Rhodri Davies and Mark Hoggard.

Students I have supervised or co-supervised: Ruby Turner (ANU); Tom New (University of Sydney, with Maria Seton, Rhodri Davies and Ben Mather); Hojatollah Shirmardgouravan (University of Sydney, with Dietmar Muller); Ishan Tandon (ANU).