20 Years of SeisSol: Calculating Earthquake Waves and More

Earthquakes generally occur where tectonic plates are moving beneath the Earth’s surface. If these plates are unable to slide past one another, stress builds up along their boundaries. Layers of rock eventually give way, releasing the stored energy abruptly. In such cases, geophysicists refer to transform earthquakes. Subduction faults are even more destructive: following a collision, one plate is forced beneath another. Anyone investigating the causes of these phenomena examines not only seismic waves but also stress and fracture energy, as well as rock deformation, temperature, subsurface conditions and many other factors. And there are equations for all of these phenomena: “SeisSol uses a discontinuous Galerkin discretisation which, up to 2006, had been applied in computational fluid dynamics but not yet in geophysical applications,” explains geophysicist Alice-Agnes Gabriel, now a professor at the Scripps Institution of Oceanography at the University of San Diego.

She joined the SeisSol team at LMU in 2012 after completing her doctorate, developed a range of simulations, refined fundamental equations, and completed her habilitation in this field in 2022. The code itself, however, originally goes back to mathematician and geophysicist Martin Käser, who, together with numerical analyst Michael Dumbser at LMU, computationally analysed the propagation of seismic waves. Alongside differential equations used to solve questions of wave propagation and earthquake rupture dynamics, they employed the discontinuous Galerkin method (DG) and combined it with a new time-stepping approach, the Arbitrary High-Order DERivatives (ADER) method. The ADER-DG method is based on equations for representing seismic waves: first, these equations are broken down into many small spatial and temporal segments to allow for more precise calculations of wave behaviour; second, the numerous computational steps can be efficiently distributed across thousands of nodes in high-performance computing (HPC) systems. Non-specialists can imagine the difference between DG and ADER-DG results as similar to that between a photograph and a film: while DG calculations produce many individual snapshots, ADER-DG provides smooth transitions. Waves are represented more realistically, making it possible to see what accelerates ground motion.

“In SeisSol, the aim is to simulate the fundamental physics of earthquakes and tsunamis on the basis of sensor data and measurements, and to develop an ever deeper understanding of the processes involved,” Bader summarises. “When I took up my professorship in 2011, the opportunity arose to work on SeisSol. Right from the start, I was drawn to the combination of high-order numerics, high-performance computing and geophysics, which motivated me to engage more closely with SeisSol.”

In its early years, SeisSol made a name for itself as a numerical method. With the addition of ever more parameters, it became increasingly complex, even representing wave propagation in three dimensions: “In 2012, SeisSol was the most accurate, but also the most computationally expensive, computer-based method for solving the wave equation,” Gabriel reports. The code delivered precision, but required substantial computing time. This was set to change when the interdisciplinary team from LMU, TUM and LRZ—funded by the German Research Foundation—came together to revise it. At the time, the geophysicists at LMU were considering how to model fault rupture processes and integrate them into SeisSol. The computer science teams at TUM and LRZ also planned to release the code as an open-source programme so that more researchers could use SeisSol: “The previously closed, numerically focused code,” says Gabriel, “became an application-oriented one, with which we were able to model a real earthquake for the first time.”

To achieve this, the geophysics team expanded the set of equations and algorithms to include rupture processes. Meanwhile, computer scientists at TUM, working with specialists at LRZ, optimised the code: algorithms were improved using adaptive time-stepping methods, while performance-critical computations were replaced by new programming approaches. SeisSol’s input and output routines were modernised, and by 2023 the former Fortran code had finally evolved into a C++ package with elements of Python.

The tremors of the Landers earthquake were calculated as the first simulation at petascale level: SeisSol proved highly scalable and efficient, achieving several quadrillion floating-point operations per second (petaflops) on SuperMUC. With this work, the team reached the final of the prestigious Gordon Bell Prize of the Association for Computing Machinery in 2014. Three years later, it won the Best Paper Award at the Supercomputing Conference (SC) in the United States for the longest and most comprehensive simulation to date of the 2004 Sumatra–Andaman earthquake, also computed on SuperMUC.