Earth modeling is extraordinarily complex – and when the earth quakes, it’s no different. The mechanisms of many earthquakes consistently flummox seismologists, making it difficult to understand why those earthquakes occurred, let alone perform the most important and elusive function in seismology: earthquake prediction. In a talk during ISC21, Alice-Agnes Gabriel, a professor of seismology at the Ludwig Maximilian University of Munich, explained how HPC-enabled earthquake codes have worked to close these gaps and take advantage of quickly accelerating supercomputing capabilities.

“Computational seismology has been a pioneering field for high-performance computing,” Gabriel said, “and [it] has itself been pioneered by HPC in order to image Earth’s interior to understand the dynamics of the metal, for example, and also to track down energy resources.” Computational seismology, she explained, just made sense for HPC: seismology is data-rich, and many of its dynamics could often be treated as linear systems.

On the other hand: “If you’re thinking about computational earthquake seismology, we’re facing a very different picture. Computational earthquake seismology is trying to understand earthquakes in a physics-based manner – we’re not imaging earthquakes, but we’re trying to model earthquakes based on physical first-order principles. We’re solving for spontaneous dynamic earthquake rupture as non-linear interaction of frictional failure across prescribed zones of weakness that are geological faults, that can take complicated shapes and interact with each other … and these frictional failure processes are linked to seismic wave propagation in a non-linear manner.”

So, in layperson’s terms: hunting for oil deposits, for instance, is easier than figuring out how many earthquakes occurred. And, indeed, there are many baffling earthquakes: earthquakes that do more damage than seismologists expect, earthquakes that occur along seemingly low-tension faults and even earthquakes that jump.

To model earthquakes, Gabriel explained, seismologists combine knowledge from a wide range of fields: seismology, geodesy, geology, tectonophysics, hydrology, numerical computing, data science, machine learning, applied mathematics, continuum mechanics, tribology, rock mechanics, materials science, and engineering. They amalgamate these processes in bootstrapped models, developing solvers that are able to read data from fault stresses and output rupture dynamics, ground deformation and synthetic observables (like simulated seismograms).

These models are robust, but three major challenges remain: first, earthquake source processes are very ill-constrained and highly non-linear; second, researchers are still working to identify which physical processes are dominant and relevant in real earthquakes and how to computationally budget for those processes’ inclusion in models; and third, teams struggle to assimilate all available knowledge in a manner both hardware- and software-suitable for modern HPC systems.

Progress on these challenges has, thankfully, been rapidly accelerating. Gabriel cited some of her team’s prior research from 2014 – itself a Gordon Bell Prize finalist – that modeled a large 1992 earthquake in California that confused researchers by seemingly jumping from fault to fault without a continuous connection. Ever since, she said, researchers had been modeling earthquakes at a petaflop scale, linking spontaneous dynamic rupture simulations with fault geometry and a range of other physical fields.

In the intervening seven years, of course, the scale is increasing: researchers are working to port those models to ever-larger spatiotemporal ranges in order to capture megathrust earthquakes (which represent all recorded earthquakes over magnitude 9.0) and tsunamis. “The extreme scale of megathrust earthquakes and tsunamis,” Gabriel said, “also requires extreme-scale runs.”

The investments required for such runs are rapidly decreasing. A “hero run” simulation of the 2004 Sumatra earthquake and tsunami, Gabriel said, took 14 hours on LRZ’s SuperMUC Phase 2 (2.8 Linpack petaflops) in 2017; now, typical runs take under six hours on just 16 nodes of LRZ’s SuperMUC-NG (19.5 Linpack petaflops).

Gabriel’s team uses a seismic solver called – appropriately – SeisSol. Using SeisSol, the researchers can solve seismic wave equations on unstructured tetrahedral meshes, combining complex geometries with heterogeneous inputs and high resolutions. Initially, SeisSol was simply a Fortran-based solver with MPI parallelization. Over time, it evolved: in 2014, for simulations of the 1992 California quake, the team implemented hybrid MPI/OpenMP parallelization, parallel I/O and a multiphysics offload function for many-core architectures, allowing the simulations to reach 96 billion degrees of freedom at 200,000 time steps.

For the Sumatra quake, SeisSol was further optimized through a bevy of features, including asynchronous I/O, overlapping computation and communication, and a feature called cluster-based local time-stepping that allowed each element to have its own time step. These optimizations allowed SeisSol to expand to 111 billion degrees of freedom at 3.3 million time steps, offering a 14-fold speedup and allowing the team to complete Sumatra simulations in hours rather than days. Essentially, Gabriel said, these speedups were the difference between a feasible simulation and an impossible simulation.

“We’ve been using this approach to shed light on puzzling earthquake observations,” Gabriel said. Lately, for instance, the researchers have been studying the 2016 Kaikōura earthquake, a magnitude 7.8 event that struck New Zealand to devastating – and befuddling – effect, perhaps most notably when the earthquake leapt from one fault to another fault 15km away, three times the previously thought maximum distance for a “jumping” earthquake. Other questions stumped researchers, too: who did the earthquake move so slowly? Why did the Hope Fault, which was thought to pose the greatest seismic hazard in the region, barely react?

SeisSol, it turns out, was able to answer many of these questions that would have dead-ended researchers a decade ago. Using the solver, the team was able to show, for instance, that another fault – the Point Kean fault – acted as a crucial link between the two endpoints of the jump, and that unfavorable dynamic stresses on the Hope Fault led to its lack of activity.

The researchers are excited to develop even more functionality, exploring the crevices of earthquake modeling. Gabriel elaborated on how her team is modeling manmade earthquakes, such as the 2017 magnitude 5.5 earthquake in South Korea, and modeling the 2018 Sulawesi earthquake that devastated Indonesia – which was complicated by the unavailability of close seismometers recording the event, leading Gabriel and her colleagues to turn to sparse data from satellites and social media.

Improvements in the code, of course, are complemented by improvements in supercomputing. Gabriel’s team has made extensive use of top-tier supercomputers, including modern leaders like SuperMUC-NG, the 5.5-Linpack petaflops Shaheen II at KAUST, the 5.4-Linpack petaflops Mahti at CSC and the 23.5-Linpack petaflops Frontera at TACC.

All of this will be amplified by the forthcoming GPU port of SeisSol – and, of course, the imminent shift to the exascale era will open up a wide range of further possibilities. Gabriel mused on what they could learn from running the models a thousand times faster: higher resolutions, uncertainty quantification, rapid earthquake response simulations, more multi-physics components – a paradigm shift toward more holistic earthquake system modeling. “The interplay of advances in high-performance computing and dense observations,” Gabriel said, “will allow us to go beyond scenario-based analysis.”