How to Simulate an Earthquake

It goes without saying that the physics of an earthquake is complicated. This, in combination with the complex structure of the earth’s crust, makes the effects of earthquakes difficult to predict. Luckily, high performance computing and numerical models are the perfect tools to simulate complex earthquake scenarios!

In my first blog post I briefly mentioned that I would work on this project – an HPC wave propagation solver called SEM3D. This software is developed at my host institution and efficiently solves the 3D wave-propagation problem using the Spectral Element Method (a type of Finite Element Method) together with a time-stepping algorithm.

So far so good, but how does this actually work? In this blog post, I will try to summarize my first few weeks of learning, and hopefully give a good high-level explanation of how earthquakes can be simulated in an HPC environment!

Note: The Finite Element Method (FEM) is a popular method for solving differential equations numerically by splitting a large domain of interest (such as the earth’s crust) into smaller, simpler elements (for example cubes). The solution to the differential equation is then approximated by simpler equations over the finite elements. In this blog post, I will assume that the reader has some basic knowledge of FEM and numerical methods for solving (partial) differential equations.

Semi-Discretization of the Wave Equation

The wave equation is a time-dependent second-order partial differential equation, meaning that in addition to the spatial variables, one must also take time into account. Anyone that has ever written FEM code is fully aware of the increased implementational efforts required when going from 1D to 2D, or worse, 2D to 3D.

Since we’re already working with a 3D domain (the earth’s crust), we will not treat time as the fourth dimension in our FEM code, but instead, opt for semi-discretization and a time-stepping scheme. This essentially means that we formulate the spatial dimensions of the equation as a FEM-problem, but express time propagation as an ordinary differential equation which can be solved with a simpler numerical time-stepping scheme.

Spectral Element Method

The Spectral Element Method (SEM) is very popular for simulating seismic wave propagation. Unlike the simpler finite difference method, it can handle the complexity of 3D earth models, while (for this problem) offering better spatial discretization, lower calculation costs, and a more straightforward parallel implementation than normal FEM.

SEM and FEM are very similar in many aspects – both, for instance, uses conforming meshes to map the physical grid to a reference element using nodes on the elements and basis functions. These are some of the key features that differ SEM from FEM:

• SEM typically uses higher-degree polynomials for the vector fields, but lower-degree polynomials for the geometry (whereas FEM often uses the same lower-degree polynomials for both)
• The interpolation functions are Legendre polynomials defined in the Gauss-Lobatto-Legendre (GLL) points (see figure).
• SEM typically uses hexahedral elements in combination with the tensor product of 1D basis functions
• For numerical integration, SEM uses the GLL rule instead of Gaussian quadrature. The GLL rule is not exact, but accurate enough.

Although it may sound strange to limit the implementation to hexahedra, GLL points, and inexact integration, the combination of the three leads to an enormous benefit: a fully diagonal mass matrix. I really want to emphasize that this simplifies the algorithm and complexity drastically.

Time Stepping

After obtaining the mass and stiffness matrices M and K from the semi-discretization and SEM-formulation, we’re left with a second-order ordinary differential equation Mu”+Ku=F which allows us to propagate over time.

There exist several different time-stepping schemes, most of which are relatively straightforward. SEM3D uses a second-order accurate velocity Newmark scheme in which the velocity (u’) at the next time step is a function of the known velocity and displacement at the previous time step.

Stability

As in all numerical time-stepping schemes, stability (which depends on the time step, dt) is crucial. Since a large time step means less number of iterations and fewer calculations, we want to use the largest possible time step without losing stability. This is where I come in!

The classical CFL condition can be estimated for homogenous materials and structured meshes. When adapted to heterogeneous media, however, the method often becomes unstable. To avoid instabilities, SEM3D uses a safety factor which results in an extremely small time step, making propagation much slower than necessary.

The better alternative, and my job to implement, is an advanced stability condition based on the largest eigenvalue of the matrix M^-1K. This advanced stability condition will allow us to calculate the maximum allowed time step without the safety factor, which means that the overall performance of the wave propagation solver will increase significantly!

---

That’s it for now – hopefully, you have a rough understanding of how seismic wave propagation can be simulated! In my next blog post, I will focus more on my own work and dive deeper into the implementation aspects of the advanced stability condition in Fortran 90 and an HPC environment.