Instantaneous mantle convection is governed by the Stokes equations for a highly viscous, incompressible fluid. Geodynamic studies solve that boundary-value problem repeatedly, across viscosity structures, plate velocities and slab configurations, and every one of those solves is a saddle-point finite-element problem. A neural operator learns the solution map directly, which turns a repeated solve into a repeated forward pass.
Background
Stokes flow drops the inertial terms from Navier-Stokes, leaving a balance between viscous stress, pressure gradient and buoyancy, with incompressibility as a constraint. The constraint is what makes it expensive: velocity and pressure must be solved together, and the resulting saddle-point system needs careful preconditioning to converge at all.
For a single model that cost is acceptable. For a parameter sweep it dominates, and sweeps are the normal mode of geodynamic work, because the interesting questions are about how flow responds to structure rather than about any single configuration.
Method
The surrogate learns a mapping between function spaces rather than between vectors: forcing and boundary conditions in, velocity field out. That distinction matters practically, because a network that learns the operator can be evaluated on a different discretisation from the one it trained on, where a network that learns a fixed grid cannot [1].
Incompressibility can be handled two ways. Enforcing it softly, as a penalty in the loss, is simple and leaves small divergence errors. Enforcing it structurally, by predicting a stream function and taking its curl, makes the constraint exact by construction at the cost of a less direct parameterisation. The second is the more defensible choice when downstream use depends on mass conservation.
Results
Two forcing modes dominate the sweep: plate velocity, which drives corner flow beneath the surface, and slab pull, which drives downwelling at the hinge. The surrogate reproduces both, and a forward pass costs milliseconds against a full saddle-point solve.
The more useful result is where it fails. Error is not distributed evenly across the domain: it concentrates in the boundary layers and at the slab hinge, precisely where velocity gradients are steepest. That is the expected behaviour of a smooth function approximator on a field with thin high-gradient structures, and it means an aggregate error figure understates the risk. The regions where the surrogate is least reliable are the regions a geodynamicist most wants to interrogate.
Surrogate forward pass
Dominant forcing modes swept
Incompressibility residual
Where error concentrates
Discussion
Three limits belong with any use of this surrogate.
It returns an instantaneous snapshot. Time evolution requires chaining snapshots, and chained surrogate calls accumulate error in a way a single call does not, so a trajectory is a much stronger claim than a snapshot.
Its accuracy is spatially structured rather than uniform, which argues for reporting error fields rather than scalars, and for escalating to the full solve when a question concerns the hinge or a boundary layer specifically.
And it inherits its training distribution. A viscosity structure unlike anything in the training set is out of distribution, and the surrogate has no way to know that unless it is asked to report uncertainty alongside its prediction.
Limitations
The demonstration sweeps two forcing modes over a single geometry family. Generalisation across geometry families is untested here and should not be assumed. The error field shown is an illustrative proxy with the correct qualitative structure, not a measured residual against a reference solve.
Key takeaways
- Stokes flow is expensive because incompressibility couples velocity and pressure into a saddle-point system, and geodynamic work sweeps parameters rather than solving once.
- A neural operator learns a mapping between function spaces, so it can be evaluated on a discretisation different from the one it trained on.
- Enforcing incompressibility through a stream function makes the constraint exact by construction, which matters when downstream use depends on mass conservation.
- Surrogate error concentrates in boundary layers and at the slab hinge, so an aggregate error figure understates the risk in exactly the regions of interest.
- The surrogate returns an instantaneous snapshot; chaining for time evolution accumulates error and is a substantially stronger claim.
References
[1] Li, Z. et al. Fourier Neural Operator for Parametric Partial Differential Equations. arXiv:2010.08895.
