ENVIRONMENT AND ENERGY

The earthquake simulation software “SeisSol” developed at LMU Munich enables researchers and civil engineers to explore the fundamentals of earthquake physics as well as the assessment of seismic hazards for fault systems with complex geometries, normally inhibited by the lack of sufficiently dense data networks.

Modeling the frictional failure of rock (on meter scale) coupled to seismic wave propagation (on scales of hundreds of kilometers) requires a multi-scale and multi-physics approach for which high-performance computing is indispensable. A pioneering petascale simulation of the 1992 Landers, California, earthquake, accompanied by detailed analysis of rock deformation effects in segmented fault systems, yields complex rupture dynamics and ground motion frequencies up to 10 Hz relevant for civil engineering purposes.

A petascale scenario of the 1992 Landers earthquake

Understanding the physics of earthquake rupture, occurring on multiple scales and at depths that cannot be probed directly, is a ‘Grand Challenge’ of Earth sciences. Geophysicists at LMU use the in-house developed SeisSol earthquake simulation software (www.seissol.org) to improve fundamental comprehension of earthquake dynamics by numerical simulation of complicated wave and rupture phenomena. The research is motivated by the need to understand and, if possible, anticipate consequences of earthquakes, which not only strike without warning, but also have disastrous after-effects.

In the absence of any realistic hope to be able to deterministically predict earthquakes, the forecasting of strong ground motions for likely earthquake scenarios will become the method of choice in the coming decades to mitigate earthquake-related damage to our infrastructure, society and economy. Also, exploration seismology relies on seismic simulations for model validation and survey design. Furthermore, modern forward simulation tools provide controlled testing environments of various plausible modeling parameters, as for example poorly constrained details of stress and friction conditions governing natural fault zones

However, the calculations involved in this kind of simulation are so complex that they push even supercomputers to their limits. The complete simulation of earthquake dynamics involves the simultaneous computation of physics-based fault slip, the subsequent seismic wave propagation and the resulting ground motions.

Modeling dynamic earthquake rupture (on meter scale) coupled to seismic wave propagation (on scales of hundreds of kilometers) requires a multiscale and non-linear multi-physics approach, and the incorporation of very complex model geometries. Realistic model setups should acknowledge topography, 3D geological structures, rheology, and fault geometries with appropriate stress and frictional parameters (Fig. 1), all of which contribute to complex ground motion patterns. Thus, LMU scientists incorporate information about the structure of the Earth, the geometry of existing fracture planes, laboratory experiments of frictional sliding between rocks, as well as the present stress conditions inside the Earth into their numerical models. The numerical models constrained in this way quickly exhibit millions of elements, making the use of high-performance computing infrastructure indispensable.

A fundamental challenge is to model the high frequency content of the three-dimensional wave field, affected by small-scale features inside the Earth. Given the lack of in-situ observations from the seismogenic zone, tens of lower resolution test scenarios (up to 1 Hz) were performed with varying fault strengths and tectonic initial conditions. The model outputs were compared to earthquake observation yielding in a best model, which was subsequently ran as a multi-petaflop simulation. This enabled the observation of high-detail rupture evolution and ground motion frequencies up to 10 Hz, relevant for civil engineering purposes (Fig. 2).

The software package SeisSol (e.g., Pelties et al., 2014) is based on an Arbitrary High­Order Derivative Discontinuous Galerkin (ADER­-DG) method using a purely local scheme which involves many matrix-matrix multiplications with sparsity patterns ranging from very sparse to nearly dense, all known a-priori.

The use of tetrahedral elements facilitates the mesh generation, specifically for complex fault geometries and topography.

In cooperation with Prof. Michael Bader of the Technical University of Munich, and co-workers, SeisSol has been optimized for the latest CPU architectures. It delivers high performance on a single compute node as well as on thousands of compute nodes in parallel. For example, it was the first simulation software to exceed 1 PFlop/s performance on the SuperMUC supercomputer (Breuer et al. 2014), and it scaled up to 1.6 million compute cores of the Tianhe-2 supercomputer (Heinecke et al. 2014).

For the Landers earthquake scenario, SeisSol was optimized to account for heterogeneous supercomputing architectures comprised of a vast array of co-processing and processing hardware via a heterogeneous solver structure. Very high scalability is achieved by avoiding complicated or even global communication and focusing on direct neighbor data exchanges. For the geometrically highly complicated 1992 Landers fault system, this pioneering dynamic rupture simulation could be performed on an adaptive mesh with 191,098,540 tetrahedrons, executing 234,456 time steps to reach 42 seconds in simulation time. The total runtime of the simulation was 7 hours and 15 minutes on all 9,216 nodes of SuperMUC Phase 2. The production run on SuperMUC sustained a performance of 1.25 Pflops (96.7% of the strong scaling performance) including all overhead, such as the MPI- initialization, mesh-reading, setup of element-local matrix structures and output. The respective 39.3% peak efficiency demonstrates that supercomputers based on commodity CPUs fully profit from the Xeon-Phi-motivated code optimizations (Breuer et al., 2014, Heinecke et al., 2014).

Coupling faulting to rock deformation

The high-resolution scenario simulation of the 1992 Landers earthquake reproduces many of the main characteristics, such as the overall rupture duration, the activation of all fault segments and the seismic moment of the real earthquake. However, the material surrounding the fault was assumed to behave purely elastically, although in nature, rocks are exposed to permanent deformation, e.g. the formation of cracks, whenever stresses reach a critical value. During a seismic event, high stresses around the fault are typically observed, potentially leading to inelastic failure of the material (e.g., Gabriel et al., 2013)

Extending simulation software to include realistic physical assumption on rheology is crucial to improve earthquake scenarios for seismic hazard analysis. However, including off-fault plasticity additionally increases the computational demands. LMU PhD student Stephanie Wollherr recently incorporated off-fault plasticity in the form of visco-plastic material response, using a non-associated Drucker-Prager plastic yield criterion, into the software package SeisSol successfully and efficiently (Wollherr & Gabriel, 2016), enabling now the analysis of parts of seismic energy being fed into the plastic deformation process, altering the source mechanics and the resulting ground motions.

The Landers model is discretized here using 22 million tetrahedral elements, thus resolving a lower frequency range in the domain than above but high-resolution fault dynamics. The elastic simulation with hybrid MPI/OpenMP parallelization ran for 3 hours and 35 minutes on 240 MPI nodes with each 28 OpenMP threads on supercomputer SuperMUC Phase 2 for a simulated time of 40 s. An additional 16 nodes were used for output only (Rettenberger et al. 2016). The computational cost for exactly the same setup including plasticity increased the computational time by only 7.2% in comparison to the purely elastic one. In comparison, comparable state-of-the-art codes (Roten et al. (2016)) report increased computational cost by plastic yielding by 67%.

A direct comparison of elastic and plastic earthquake scenarios exhibit distinct mechanisms of rupture transferring between fault segments: In comparison to fully elastic simulations, plastic rock deformation does not only influences the step-over timing and location, but also favours triggering of rupture by seismic waves instead of direct fault branching. Field observations report “flower-like” damage zones surrounding the thin slipping interfaces of faults. The plastic strain in SeisSol’s simulation consistently increases close to the free surface (Fig. 3, cross section). Furthermore, the results indicate that the extension of plastic strain around the fault system increases in areas of higher structural complexity and at locations where earthquakes “jump” to adjacent segments (Fig. 3). These results may directly shed light on earthquake mechanics from observations.