Real-time Bayesian Inference With Full Physics Models: A Digital Twin for Tsunami Early Warning

Although rare, tsunamis can be extremely deadly and cause catastrophic socioeconomic losses that amount to hundreds of billions of dollars in damage. State-of-the-art tsunami early warning systems typically infer earthquake moment magnitude and hypocenter location, then use precomputed scenarios or simplified fault models to trigger tsunami forecasts within minutes. However, the seismic information on which these systems rely is based on simplified assumptions that often fail to capture the complexities of earthquake rupture dynamics [6], which can lead to delayed, missed, or false warnings.

To address some of the shortcomings of existing systems, we present a full-physics Bayesian inversion and prediction framework—a so-called digital twin [8]—that enables accurate, data-driven tsunami forecasting with dynamic adaptivity to complex source behavior. Our approach aims to improve warning capabilities by directly integrating real-time ocean bottom pressure data to constrain tsunami sources without relying on assumed fault geometry or magnitude. We utilize the acoustic-gravity wave equations—a partial differential equation (PDE) model that couples the propagation of ocean acoustic waves with surface gravity waves [7]—to model the tsunami dynamics of our digital twin. Based on this physics model, the inversion framework employs ocean bottom pressure data to statistically infer earthquake-induced spatiotemporal seafloor motion in real time. The solution of this Bayesian inverse problem is then used to simulate tsunami propagation toward coastlines and forecast wave heights at specified locations with quantified uncertainties.

This work falls into an important class of problems called inverse problems, which are fundamental to many fields of science and engineering. In a forward problem, the input parameters (e.g., initial or boundary conditions, coefficients, or source terms) are specified; the governing equations’ solution then yields the state variables that describe the system’s behavior. However, forward models typically contain uncertain parameters that are often infinite-dimensional fields that can vary in space (and possibly time). Inverse problems hence seek to infer uncertain parameter fields from measurements or observations of system states.

Infinite-dimensional inverse problems are usually ill posed, in that many parameter fields may be consistent with the observational data within the noise [2]. The Bayesian approach treats the uncertain parameter as a random field, and the inverse problem’s solution is the posterior probability distribution, which characterizes the probability that any parameter field is consistent with the data and prior knowledge. As such, the Bayesian framework is a very powerful technique for the characterization of uncertainty in the inverse problem’s solution.

But this power comes at a price: computing the full Bayesian solution is generally intractable for infinite-dimensional inverse problems that are governed by PDEs. For example, computing just the mean of the posterior formally necessitates numerical integration in the discretized parameter dimension; each evaluation of the integrand requires solution of the forward problem—which can take hours—and we need numerous forward solutions. The specific Bayesian inverse problem in question is the inference and prediction of tsunamis for early warning, which adds a real-time requirement to the aforementioned challenges. Our forecasting goal—to provide early warning in less than a minute—may seem futile, as state-of-the-art inversion methods would call for many years of supercomputing time to solve this problem. However, we achieve real-time performance through a combination of novel parallel inversion algorithms, advanced PDE discretization libraries, and the exploitation of leading-edge graphics processing unit (GPU) supercomputers.

Our digital twin framework [5] combines the following key ideas: (i) offline-online decomposition of the inverse solution that precomputes various mappings to enable the real-time application of the inverse operator and circumvent the need for PDE solutions during inversion; (ii) representation of the inverse operator in the data space instead of the high-dimensional parameter space; (iii) extension of this real-time inference methodology to the goal-oriented setting to predict quantities of interest (QoI) under uncertainty; (iv) exploitation of the shift invariance for the discretized differential operator to reduce the number of PDE solutions by several orders of magnitude; and (v) exploitation of this same property to perform efficient fast-Fourier-transform-based and multi-GPU-accelerated Hessian matrix-vector products (matvecs) [9].

The solution of inverse problems requires numerous Hessian matvecs, each of which comes at the cost of a pair of forward/adjoint PDE solutions. PDE solutions are computationally expensive for large-scale problems, which means that we cannot perform them in real time. One key innovation that enables real-time inference is the recognition and exploitation of the fact that the parameter-to-observable (p2o) map is governed by an autonomous dynamical system, which takes the seafloor velocity (parameter field) as input, solves the acoustic-gravity wave equations, and extracts the pressures (observables) at the sensor locations and observation times. Specifically, the acoustic-gravity model represents a linear, time-invariant dynamical system for which we can write the discrete p2o map as a block lower triangular Toeplitz matrix that corresponds to a time-shift invariance of its blocks. This invariance implies that (i) we can precompute the p2o map from only as many adjoint PDE solutions as the number of sensors; (ii) we can embed the block Toeplitz matrix within a block circulant matrix that is block-diagonalized by the discrete Fourier transform; and (iii) the p2o matvec becomes a block-diagonal matvec operation in Fourier space.

Scalability is critical for the most expensive part of the offline computations: the adjoint PDE solutions. We discretized these adjoint wave propagation problems with the MFEM finite element library [1]. Next, we ran MFEM-based GPU simulations on large computer systems, including Lawrence Livermore National Laboratory’s El Capitan — the highest-ranked supercomputer on the Top500 list. Here, the MFEM solver demonstrated excellent scalability with up to 43,520 GPUs. The biggest problem involved more than 55.5 trillion degrees of freedom, which appears to be the largest reported unstructured mesh finite element computation.

We applied our digital twin framework to the Cascadia Subduction Zone (CSZ), which stretches 1,000 kilometers from northern California to British Columbia (see Figure 1). Because its last catastrophic seismic event was 300 years ago, seismologists consider the CSZ to be overdue for a magnitude 8.0-9.0 megathrust earthquake; paleoseismic and geodetic evidence indicates a recurrence interval of roughly 250 to 500 years. To assess our framework’s feasibility for near-field, real-time tsunami forecasting, we generated realistic seafloor displacements that capture the nonlinear interaction of seismic wave propagation and frictional rock failure during a magnitude 8.7 CSZ dynamic rupture scenario that we computed with the earthquake simulation code SeisSol [3]. We then used these “true” seafloor displacements as a source for our tsunami model to generate synthetic observational data with added noise at 600 different hypothesized seafloor pressure sensors. The sparse, noisy observations of acoustic pressure allowed us to solve the PDE-based inverse problem for the spatiotemporal parameter field with more than one billion parameters—which represent seafloor velocities during the earthquake—and make predictions of the QoI (sea surface wave height) at 21 forecast locations.

The offline phase, which we performed just once, required 621 adjoint wave propagations and took 538 hours on 512 GPUs. Given real-time data, the online phase requires just 0.2 seconds to exactly solve the billion-parameter inverse problem. Figures 2 and 3 illustrate the results of the parameter inference and QoI prediction under uncertainty.

Given a sufficient number of ocean bottom sensors, this work ultimately enables real-time probabilistic tsunami forecasts for near-field events where conventional systems may fail due to delayed or inaccurate source characterization, such as a future magnitude 8.0–9.0 CSZ megathrust earthquake. Our digital twin framework accounts for complex rupture dynamics and tsunami wavefields in real time, under uncertainty, and without sacrificing the predictive power of the full three-dimensional wave propagation equations. This structure enables early predictions of inundation and evacuation zones within seconds, which is critical for life-saving responses where warning time is otherwise measured in minutes.

Autonomous dynamical systems arise in many different settings beyond geophysical inversion. Our Bayesian inversion-based digital twin framework is more broadly applicable to acoustic, electromagnetic, and elastic inverse scattering; source inversion for the transport of atmospheric or subsurface hazardous agents; satellite inference of emissions; and treaty verification, among other settings.

This article is based on Omar Ghattas’ Ivo & Renata Babuška Prize lecture at the 2025 SIAM Conference on Computational Science and Engineering, which took place last year in Fort Worth, Texas. The corresponding paper [4] was awarded the Association for Computing Machinery’s Gordon Bell Prize at the 2025 International Conference for High Performance Computing, Networking, Storage, and Analysis (SC25), which was held in St. Louis, Mo., last November.

Acknowledgments: This research was supported by U.S. Department of Energy (DOE) Advanced Scientific Computing Research grant DE-SC0023171; U.S. Department of Defense Multidisciplinary University Research Initiatives grants FA9550-21-1-0084 and FA9550-24-1-0327; U.S. National Science Foundation grants DGE-2137420, EAR-2225286, EAR-2121568, OAC-2139536, and RISE-2425922; and European Union 2020 Horizon grants 101093038, 101058518, and 101058129. Work was performed under the auspices of the DOE by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344 (LLNL-PROC-2004753). This research utilized resources from the Swiss National Supercomputing Centre, the Texas Advanced Computing Center, the National Energy Research Scientific Computing Center (ALCC-ERCAP0030671), and Livermore Computing.

References 
[1] Anderson, R., Andrej, J., Barker, A., Bramwell, J., Camier, J.-S., Cerveny, J., ...  Zampini, S. (2021). MFEM: A modular finite element methods library. Comput. Math. Appl., 81, 42-74.
[2] Ghattas, O., & Willcox, K. (2021). Learning physics-based models from data: Perspectives from inverse problems and model reduction. Acta Numer., 30, 445-554.
[3] Glehman, J., Gabriel, A., Ulrich, T., Ramos, M., Huang, Y., & Lindsey, E. (2025). Partial ruptures governed by the complex interplay between geodetic slip deficit, rigidity, and pore fluid pressure in 3D Cascadia dynamic rupture simulations. Seismica, 2(4).
[4] Henneking, S., Venkat, S., Dobrev, V., Camier, J., Kolev, T., Fernando, M., Gabriel, A.-A., &  Ghattas, O. (2025). Real-time Bayesian inference at extreme scale: A digital twin for tsunami early warning applied to the Cascadia subduction zone. In SC ’25: Proceedings of the international conference for high performance computing, networking, storage and analysis (pp. 60-71). St. Louis, MO: Association for Computing Machinery.
[5] Henneking, S., Venkat, S., & Ghattas, O. (2025). Goal-oriented real-time Bayesian inference for linear autonomous dynamical systems with application to digital twins for tsunami early warning. Preprint, arXiv:2501.14911.
[6] LeVeque, R.J., Bodin, P., Cram, G., Crowell, B.W., González, F.I., Harrington, M., … Wilcock, W.S.D. (2018). Developing a warning system for inbound tsunamis from the Cascadia subduction zone. In OCEANS 2018 MTS/IEEE Charleston (pp. 935-944). Charleston, SC: Institute of Electrical and Electronics Engineers.
[7] Lotto, G.C., & Dunham, E.M. (2015). High-order finite difference modeling of tsunami generation in a compressible ocean from offshore earthquakes. Comput. Geosci., 19(2), 327-340.
[8] National Academies of Sciences, Engineering, and Medicine. (2024). Foundational research gaps and future directions for digital twins. Washington, D.C: National Academies Press. Retrieved from https://nap.nationalacademies.org/catalog/26894/foundational-research-gaps-and-future-directions-for-digital-twins.
[9] Venkat, S., Fernando, M., Henneking, S., & Ghattas, O. (2025). Fast and scalable FFT-based GPU-accelerated algorithms for block-triangular Toeplitz matrices with application to linear inverse problems governed by autonomous dynamical systems. SIAM J. Sci. Comput., 47(5), B1201-B1226.