Forward and Inverse Mantle Convection with Neural Operators

Chenxi Kong; Michael Gurnis; Zachary E. Ross

REVIEW 4 major objections 4 minor 3 references

Forward and Inverse Mantle Convection with Neural Operators

T0 review · 4 major / 4 minor · reviewed 2026-08-03 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read Forecasting and reversing mantle convection with learned neural surrogates, the paper claims a joint inversion of the present-day thermal field plus surface velocity history recovers past mantle states further back in time than any other te

desk verdict Useful 2D benchmark for neural-operator thermal-state inversion, but without a numerical-solver adjoint control, the 'replace the solver' claim is not yet demonstrated. read the letter →

arxiv 2601.23178 v2 pith:WS52M5XX submitted 2026-01-30 physics.geo-ph

Chenxi Kong , Michael Gurnis , Zachary E. Ross This is my paper · ORCID

classification physics.geo-ph

keywords mantleconvectionthermalstatereconstructionneuraloperatorsFourieroperatorinverseproblemadjoint-freeinversionsurfacekinematicsauto-differentiation

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper tries to establish that Fourier neural operators can replace numerical solvers for mantle convection, both for stepping thermal states forward in time and for the ill-posed task of stepping them backward. It further claims that the most reliable way to reconstruct a past thermal state is not to reverse the dynamics directly but to run the learned forward operator in an iterative inversion that fits two sets of observations: the terminal (present-day) temperature field and the time series of surface horizontal velocity. On synthetic convection sequences at Rayleigh number 10^7, this joint inversion recovers nearly all thermal structures that existed within a one-transit-time window, and remains useful back to about 0.86 transit times when the observations are polluted with 5% pink noise. If correct, the result makes computationally feasible a class of geophysical inversions — recovering earlier mantle structure from seismic tomography and plate-kinematic reconstructions — that was previously prohibitive because each adjoint solve costs orders of magnitude more than a forward model.

What carries the argument

The Fourier neural operator (FNO), a neural operator that parameterizes its integral kernel in Fourier space, is the load-bearing object. Three instances are trained: a Stokes operator S_phi mapping temperature to velocity and pressure via purely physics-informed PDE losses; a forward convection operator F^{+nΔt}_phi mapping a thermal state to one nΔt later (and to its instantaneous velocity), trained data-driven with time steps about 100-300 times the CFL limit; and a reverse operator F^{-nΔt}_phi trained by swapping the same data pairs. Auto-differentiation through the surrogate's computational graph supplies gradients for the inversion, replacing adjoint-state solves; the joint inversion

What would settle it

Take a convection sequence computed with the finite-element solver, degrade the terminal observation to a realistic tomographic image (blurred, restricted to certain depth ranges, with 5% pink noise), and run the joint inversion; if the reconstructed initial state's correlation with ground truth drops below the clean-case threshold before -0.86 transit times, the promise for real geophysical data fails. Alternatively, test the surrogate's gradients by perturbing the initial state in a direction that leaves the objective unchanged but changes the true past state; if the optimization cannot move

Watch

Extended reading notes

Core claim

The central discovery is that a data-driven Fourier neural operator can approximate the mapping between two convecting thermal states separated by a time interval hundreds of times larger than the Courant-Friedrichs-Lewy step, and that the same architecture, trained on reversed input-output pairs, approximates the ill-posed reverse mapping without the blow-up that plagues direct numerical reversal of the anti-diffusion equation. The paper uses these surrogates to compare four reconstruction strategies and shows that the reverse convection operator, while accurate on noiseless inputs, is destroyed by 5% pink noise, whereas an inversion that minimizes misfit to the terminal thermal field alone

Load-bearing premise

The inversion experiments assume perfect observability of the full terminal thermal field and the complete surface velocity time series; if real tomographic images and plate reconstructions are too partial or indirect, the demonstrated robustness to 5% pink noise may not transfer, and the surrogate gradients could lead to surrogate-specific local minima.

Editorial extensions

If this is right

Thermal state reconstruction becomes computationally tractable: the cost of training the surrogate and running a joint inversion is comparable to performing one time-dependent inversion with conventional adjoint methods, and the advantage grows with grid resolution.
Present-day seismic tomography plus plate-kinematic histories can, in principle, be fed into this joint inversion to recover mantle structure over roughly the last transit time (~150-200 Myr) with the dominant long-wavelength features preserved.
The learned forward operator accelerates forward convection modeling by factors of 5,000-20,000 at 257x257 resolution, and the acceleration scales with problem size.
Direct reverse-convection operators provide a fast approximate backward map when inputs are clean, but require a noise-filtering preprocessing step before they can be used on real observations; the paper explicitly proposes such a denoising mapping.
Gradients computed via auto-differentiation reach machine precision relative to the surrogate and avoid the numerical errors of solving adjoint equations.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

If the joint-inversion scheme carries over to non-linear viscosity and spherical geometry, it would allow the plate-tectonic record to be assimilated as time-dependent boundary data on a whole-mantle model, effectively replaying mantle history — a step the paper leaves implicit but its cost scaling directly invites.
The Green's-function sensitivity analysis implies surface horizontal velocity constrains upper-mantle structure far more strongly than deep structure; a testable extension is to add dynamic-topography or anisotropy data to the joint objective and ask whether the -0.86 transit-time noise limit extends deeper.
The paper's cost comparison assumes training from scratch; since the same trained forward operator can serve many inversions, the economic case becomes much stronger in a production setting where repeated reconstructions are made.
A natural falsifiable prediction: on a sequence whose initial state contains thermal structure below the depth sensitivity of surface velocity, the joint inversion's reconstruction of that deep structure should degrade; if it does not, the surrogate is learning more from the terminal field than the stated mechanism suggests.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

Short version: this is a well-put-together 2D proof-of-concept for using Fourier neural operators in mantle-convection inversion, and the four-method comparison is genuinely informative. The joint inversion with surface velocity is the clear winner under noise. But the paper does not include the one control that would actually establish its strongest claim: an inversion using the numerical solver's own adjoint gradients on the same objective. Until that is done, we can't tell whether the surrogate gradients are producing a faithful solution to the inverse problem or an implicitly regularized version of it.

What is actually new: training a reverse convection operator by swapping input-output pairs is a standard trick, and FNO surrogates for convection already exist. The new bit is the systematic comparison of reverse buoyancy, reverse operator, terminal-state inversion, and joint inversion, and the finding that adding surface-velocity time series makes the inversion stable to about -0.86 transit times under 5% pink noise. The Stokes operator trained purely from physics with no data is a nice demonstration of the super-resolution property. The forward operators achieve one-step errors in the 0.1-0.9% range and stay on-attractor for several transit times. That is useful evidence that the approach has legs.

The soft spots are real but not fatal. The missing adjoint baseline is the main one. The paper even reports that terminal-state-only inversion falls into non-physical minima before t=-0.52; without a solver-based benchmark we cannot attribute that failure to the inverse problem itself or to surrogate artifacts. The cost model in the supplement has a factor-of-25 inconsistency between Eq. S30 and the 25,000 iterations in Fig. S6 — it actually makes the neural-operator workflow look even cheaper, but it should be fixed. The Green's function validation in S4 compares a no-slip neural operator to a free-slip analytical kernel, which is an apples-to-oranges comparison. No code or data are provided, and there are no repeated-seed uncertainty intervals; both should be expected for an ML paper.

Bottom line: this deserves a serious referee and will likely be useful to the geodynamics-ML community once the baseline experiment and reproducibility items are added. I'd engage with it, but cite with care until it's revised.

Referee Report

4 major / 4 minor

Summary. The paper trains tensorized Fourier neural operators (TFNO/FNO) for three elements of 2-D bottom-heated Rayleigh–Bénard convection: a purely physics-informed Stokes operator S_φ, a data-driven forward convection operator F^{+nΔt}_φ, and a data-driven reverse convection operator F^{-nΔt}_φ. All surrogates are trained and evaluated against the Underworld finite-element solver at Ra = 10^5–10^7. The authors then compare four thermal-state reconstruction methods—reverse buoyancy, reverse convection neural operator, terminal-state-only inversion, and joint inversion with terminal thermal state plus surface velocity time series—on synthetic convection sequences. Reported results include one-step forward relative L2 errors of 0.1–0.9% (Table 2), speedups of about 5000×–20000× for forward integration, and a joint inversion that remains informative to about t = -0.86 transit times under 5% pink noise (Section 3.3). The paper concludes that neural-operator workflows can make thermal-state reconstruction computationally feasible and that the total workflow cost is comparable to one traditional adjoint inversion.

Significance. If the central claims hold, the paper makes a useful contribution to surrogate-modeled mantle convection: it demonstrates that data-driven FNOs can step over CFL-limited time steps by orders of magnitude, that a physics-informed Stokes operator can be trained without precomputed training data, and that joint inversion with time-dependent surface kinematics materially stabilizes an ill-posed reversal problem. Strengths include the systematic error documentation across Rayleigh numbers and step sizes (Table 2), the independent Green's-function check of the learned surface-velocity sensitivity (Fig. S5), and the clear visualization of why reverse-operator methods fail under noise. However, the 'replacement of numerical solvers' claim is not yet fully demonstrated because no numerical-adjoint inversion baseline is provided, and the current submission lacks code, data, and any uncertainty quantification from repeated training runs.

major comments (4)

[Section 3.3, Eq. (12)–(14)] The four-method comparison lacks a control inversion in which the same objective and observations are optimized using gradients from the numerical solver (e.g., an adjoint-based inversion within Underworld). All inverse variants differentiate through the FNO forward surrogate, whose one-step relative L2 errors are 0.1–0.9% (Table 2) and which is spectrally truncated. In an ill-posed problem, these model errors can bias or implicitly regularize gradients, so the observed robustness of the joint inversion—and the failure of the terminal-state-only inversion before t = -0.52—may be partly surrogate-specific rather than intrinsic to the inverse problem. Please add a numerical-adjoint baseline for at least the joint inversion, or explicitly restrict the claim that FNOs 'replace' numerical solvers for inversion.
[Supplementary S7, Eq. (S30) and Fig. S6] The cost comparison relies on Eq. (S30), AI·AF ~ AT (~10^3), to conclude that the neural-operator workflow cost is comparable to one traditional inversion. However, Fig. S6 shows that the converged demonstration used the 25,000th optimization iteration, and the inversion window is about one transit time (AF ~ 1), so AI·AF ~ 2.5×10^4, an order of magnitude larger than AT. The equation is therefore inconsistent with the paper's own optimization history. Re-derive the cost comparison with the actual iteration count; the qualitative conclusion may survive or be strengthened, but the arithmetic as written is not self-consistent.
[Section 2.2, Eq. (12)] The synthetic observables are idealized: the terminal thermal field is known over the whole domain, and full surface horizontal velocity is known at every boundary point at N+1 discrete times. The abstract's geophysical prospect—applying the joint technique to seismic tomography and plate reconstructions—implicitly assumes that these fields provide complete coverage with only 5% noise. Real tomographic models are partial, indirect, and spatially heterogeneously resolved, and plate reconstructions sample limited, uncertain surface motions. Robustness to 5% pink noise on complete fields does not guarantee robustness to missing or indirect observations. Please discuss how the method would be reformulated for partial observations, or add experiments with partial/noisy boundary data.
[Data Availability / Reproducibility] The submission provides no code or training data, and the Data Availability section only promises them 'for the final version.' The paper also reports no repeated-seed or retraining variability for any of the key quantitative comparisons (Table 2, Fig. 5, Fig. S6). Since all central results depend on trained neural operators, this prevents verification of the results and assessment of training stochasticity. Please release the code and data, or at a minimum provide full training details, seeds, and evaluation statistics (mean/standard deviation over repeated runs) in the supplement.

minor comments (4)

[Global] Typos and small errors: 'propogation' (Introduction), 'Courant Fredrich Lewy' (Abstract and Section 2.1, should be 'Courant–Friedrichs–Lewy'), 'Rayleigh–Bernard' (Introduction, should be 'Bénard'), and 'Inverstion' in the Fig. 4 axis label.
[Figures 4 and 6] The red/blue predicted surface-velocity curves may be difficult to distinguish for color-blind readers. Please add line styles (dashed/solid) or distinct markers in addition to color.
[Section 3.2 / Table S6] The speedup discussion quotes 5000×, 10000×, and 20000× for F^{+1}, F^{+2}, F^{+4}. These are instantaneous per-transit-time comparisons on a CPU core vs. a GPU; it would be helpful to state once in the main text that the comparison is single-core CPU, and to include the training wall-clock cost for the forward operator in the main workflow cost discussion.
[Section 2.2 / Eq. (12)] The objective function is described as 'Modified from Li et al. 2017'. The modifications (addition of the Laplacian smoothing term to P, the temporal surface-velocity term, and the normalization by domain measures) should be stated explicitly relative to the prior formulation, since the regularization choices directly affect inversion performance.

Circularity Check

0 steps flagged · score 1.0 of 10

No significant circularity: the FNO surrogates are trained on and benchmarked against Underworld in the standard surrogate-validation sense, and the inversion claims do not reduce to fitted inputs or self-citations.

full rationale

The paper's load-bearing claims are (i) that FNO surrogates can approximate the Underworld Stokes/advection-diffusion dynamics, and (ii) that inversions through those surrogates can recover synthetic past thermal states from terminal-state and surface-velocity observations. Neither claim is circular by the paper's own equations. The Stokes operator S_phi is trained with a purely physics-informed loss (Eq. 21) on random thermal fields, not on Underworld outputs, and is then benchmarked against Underworld; this is independent validation, not a fitted quantity relabeled as a prediction. The forward and reverse convection operators F^{±nΔt}_phi are trained on data pairs generated by Underworld, and evaluating them on additional Underworld trajectories is the standard way to validate a learned surrogate; it does not make the benchmarked error a tautology. The inversion experiments use synthetic observations of the terminal thermal field and surface horizontal velocities as data, and the reconstructed initial state is not among the objective's inputs; minimizing Eq. 12 through the surrogate is a genuine inverse benchmark, even though the observations are synthetic and generated from the same solver family. The absence of a traditional-adjoint inversion baseline is a legitimate correctness/benchmarking concern, but it is not circularity: the joint inversion's reconstruction is not equal by construction to any fitted parameter. Self-citations such as Li et al. (2017) for the objective-function form and Conrad & Gurnis (2003) for the reverse-buoyancy baseline provide context, baselines, and cost estimates, but none of the central claims is forced by an unverified self-citation. Overall, the derivation chain is self-contained in the sense required here: the predicted quantities are not equivalent to the training targets or objective terms by definition.

Assumptions & free parameters 6 free parameters · 6 assumptions · 0 invented entities

No new physical entities are introduced. The central experiment relies on hand-tuned loss and objective weights, chosen architecture hyperparameters, and domain assumptions about training-data coverage and observability. Code and data are not yet released, so these premises cannot currently be audited independently.

free parameters (6)

Objective weights beta1-beta4 = 1.0, 2.5e-1, 1.0e-5, 2.0e-8 (Table S4)
Hand-tuned weights in Eq. 12-13 control terminal-state, surface-velocity, Laplacian-smoothing, and mean-penalty terms; all inversion results depend on them.
PDE-loss weights for Stokes S_phi = Table S2 multi-stage values (betaC, betaC1, betaC2, betaM, betaB, betaN, betaU, betaP)
Hand-scheduled weights in Eq. 21 tune the physics-informed Stokes training; the reported ~5% velocity error depends on them.
Convection-operator step sizes n*Dt = 2.5e-5 to 1e-2 depending on Ra (Table 2)
Design choice of integration interval per recursion; balances accuracy, GPU memory, and chronological constraint density, and is not derived from the physics.
Gaussian covariance length scale l for random initial fields = Not fixed numerically; 'slightly larger than boundary layer thickness' (Section 2.3.2)
Training distribution is generated from these fields; operator generalization depends on this ad hoc choice.
Low-pass filter for correlation metric = Cutoff ratio = 0.2, transition width = 0.4, smoothing power = 4 (Algorithm S1)
Selected post hoc to emphasize long wavelengths; inflates reported reconstruction correlations.
TFNO/FNO architecture hyperparameters = Table S1: modes 65/129, hidden channels 128, blocks 5/6, rank 0.1/0.25
Chosen by hand; determine model capacity and therefore the reported accuracies.

assumptions (6)

domain assumption Underworld solves the governing equations accurately enough to serve as ground truth
All training labels and validation targets come from this FEM solver; if the solver is not a valid reference, the surrogate and inversion benchmarks inherit the error (used throughout Section 3).
domain assumption Random Gaussian initial fields span the relevant thermal-state function space
Training pairs are drawn only from convection sequences started from these fields (Section 2.3.2); generalization to other states is assumed.
domain assumption Boussinesq, constant-viscosity, 2D Cartesian approximation is an adequate mantle-convection model
Eq. 1-6 and Section 2.1; the paper acknowledges 3D geometry and nonlinear rheology as future work (Section 4).
domain assumption Boundary conditions: no-slip top/bottom, periodic sides, fixed temperature top/bottom
Section 2.1; the supplementary sensitivity kernel (Eq. S13) switches to free-slip top/bottom without explicit discussion.
domain assumption Surface velocity observations are available on the entire top boundary over the reconstruction interval
The inversion formulation in Eq. 12 assumes known v_obs on L_s for all N+1 times; real plate kinematics are sparser and noisier.
domain assumption Gradients through the learned surrogate are accurate enough for optimization
Auto-differentiation through the FNO replaces adjoint solutions (Section 2.2); this relies on surrogate approximation error being small enough not to bias the descent direction.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Forward and Inverse Mantle Convection with Neural Operators." pith.science (2026). https://pith.science/paper/WS52M5XX

@misc{pith2026260123178,
  author       = {Pith},
  title        = {Pith review of: Forward and Inverse Mantle Convection with Neural Operators},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/WS52M5XX}},
  note         = {Machine review of arXiv:2601.23178}
}

read the original abstract

Thermal state reconstruction--reversing convection to recover the thermal structure of the mantle at an earlier geologic time--is an important tool to understand the evolution of mantle convection and its relation to seismic tomographic images and observations at the surface. Thermal state reconstructions are computationally expensive. Here we transformed the basic computational element, numerical solvers, into neural operators, a class of machine learning models for learning mappings between function spaces. Focusing on a specific architecture, Fourier Neural Operators, we demonstrate that they can represent not only a surrogate model like the Stokes system of equations using a purely physics informed approach, but also discover operators without explicit mathematical formulations or even ill-posedness from data, including the direct mapping between two convecting thermal states separated by a long time interval much larger than the Courant-Friedrichs-Lewy condition and its reversal. These neural operators significantly accelerate forward and inverse convection modelling by transforming forward physical processes into surrogate models with lower complexity while utilizing auto-differentiation to calculate gradients. With this framework, we demonstrate the strengths and weaknesses of four methods for thermal state reconstructions: reverse buoyancy, reverse convection operator, an inversion with only the terminal thermal state, and a joint inversion with the terminal thermal state and surface velocity evolution. The reverse convection operator is shown to perform poorly in the presence of observational noise, but the joint inversion overcomes this limitation. The joint technique could probably become a solution to large-scale thermal state inversion problems using seismic tomography and plate tectonic reconstructions.

Figures

Figures reproduced from arXiv: 2601.23178 by the authors.

**Figure 1.** An example of forward convection computed with Underworld, which is used as an evaluation data sequence in this study. The computation is initiated from a Gaussian random initial thermal field that satisfies the prescribed boundary conditions. The total integration time is around 34 transit times, during which the convection pattern evolves from a chaotic one to a relatively steady state, as shown by the tracked Nu… view at source ↗

**Figure 2.** The architecture of three neural operators described in this study. Detailed MLP and FB parameters are listed in Table S1. and (3) application of neural operators to thermal state reconstruction, including a fast and direct method using F −n∆t ϕ , and a robust inverse method using F +n∆t ϕ . 2.3.1 Stokes neural operator The Stokes neural operator, Sϕ, solves the Stokes equation for velocity and pressure from buoyanc… view at source ↗

**Figure 3.** Comparisons between forward computations using F +n ϕ7 and Underworld, Ra = 107 . Row 1 to 4 shows the forward evolution snapshots of the systems. Column 1: temperature snapshots and velocity streamlines computed by F +1 ϕ7 , system integrated by F +1 ϕ7 ; Column 2: temperature and velocity computed by Underworld, system integrated by Underworld; Column 3: velocity computed by Sϕ based on thermal fields integrated b… view at source ↗

Figures from the paper (3 more)

**Figure 4.** Figure 4: Reconstruction performances of different methods. Reversal time steps shown by row (from top to bottom backwards in time). Column 1: A ground-truth thermal convection sequence within time window 3 computed by Underworld forward in time from bottom to top; Column 2-5: R…

**Figure 5.** Figure 5: Correlation coefficient of reconstructed thermal fields with ground-truth fields versus backwards time. Colored lines denote different reconstruction methods. (a) Reconstruction with synthesized observations (no noise); (b) Reconstruction with synthesized observations…

**Figure 6.** Figure 6: Reconstruction performances of different methods. Reversal time steps shown by row (from top to bottom backwards in time). Column 1: A ground-truth thermal convection patterns within time window 3 computed by Underworld forward in time from bottom to by two inversion m…

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

3 extracted references · 2 linked inside Pith

[2017]

(though with a shorter integration time) while the latter is derived from this study (Fig. S6), we find that: AI AF ∼A T (S30) Forward and Inverse Mantle Convection with Neural Operators55 showing that the total cost of an adjoint time-dependent inversion is about the same or larger than than the dominant term in a neural operator–based cost—the training ...
[2021]

Li, Z., Kovachki, N., Choy, C., Li, B., Kossaifi, J., Otta, S., Nabian, M

Physics-informed neural operator for learning partial differential equations,arXiv preprint arXiv:2111.03794. Li, Z., Kovachki, N., Choy, C., Li, B., Kossaifi, J., Otta, S., Nabian, M. A., Stadler, M., Hundt, C., Azizzadenesheli, K., et al., 2023. Geometry-informed neural operator for large-scale 3d pdes, Advances in Neural Information Processing Systems,...

arXiv 2023
[2023]

Neural operator: Learning maps between function spaces with applications to pdes,Journal 34C. Kong, M. Gurnis and Z. Ross of Machine Learning Research,24(89), 1–97. LeCun, Y., Bottou, L., Bengio, Y., & Haffner, P., 2002. Gradient-based learning applied to document recognition,Proceedings of the IEEE,86(11), 2278–2324. Li, D., Gurnis, M., & Stadler, G., 20...

arXiv 2002

Pith tools

Reviewed August 3, 2026 · model on record in the stance chip above.

[1] [2017]

(though with a shorter integration time) while the latter is derived from this study (Fig. S6), we find that: AI AF ∼A T (S30) Forward and Inverse Mantle Convection with Neural Operators55 showing that the total cost of an adjoint time-dependent inversion is about the same or larger than than the dominant term in a neural operator–based cost—the training ...

[2] [2021]

Li, Z., Kovachki, N., Choy, C., Li, B., Kossaifi, J., Otta, S., Nabian, M

Physics-informed neural operator for learning partial differential equations,arXiv preprint arXiv:2111.03794. Li, Z., Kovachki, N., Choy, C., Li, B., Kossaifi, J., Otta, S., Nabian, M. A., Stadler, M., Hundt, C., Azizzadenesheli, K., et al., 2023. Geometry-informed neural operator for large-scale 3d pdes, Advances in Neural Information Processing Systems,...

arXiv 2023

[3] [2023]

Neural operator: Learning maps between function spaces with applications to pdes,Journal 34C. Kong, M. Gurnis and Z. Ross of Machine Learning Research,24(89), 1–97. LeCun, Y., Bottou, L., Bengio, Y., & Haffner, P., 2002. Gradient-based learning applied to document recognition,Proceedings of the IEEE,86(11), 2278–2324. Li, D., Gurnis, M., & Stadler, G., 20...

arXiv 2002