Gamra is a parallel, adaptive 2D and 3D finite difference C++ code for modeling processes in the Earth such as earthquakes, mountain building, and mantle convection.


This the calculated displacement in cm on the surface and through a cross section of a 3D model of the 2011 Mw 9.1 Tohoku-Oki earthquake. Click on the image for a larger version.

Tohoku displacement

This is a 3D model of the 1992 Landers earthquake progressively zoomed in. There are 426 individual fault segments that make up the entire model.

Landers faults 1 Landers faults 2 Landers faults 3 Landers faults 4

For more details, see the paper.

Crust and Mantle Dynamics

Gamra is also designed to model the dynamics of the crust and mantle of the Earth as a whole. The crust and mantle can be approximated as viscous fluids that flow very, very slowly. The mantle flows relatively easily, while the crust it much more rigid. A simple test case is to model a dense rigid, block of crust falling through the mantle.


This setup gives solutions for the pressure in 2D and 3D.

Block 2D Block 3D

We can see sharp spikes in the pressure near the interface between the crust and mantle. These spikes are not real. They are numerical artifacts, and they do not improve with resolution.

To remedy this, I have been experimenting with the Immersed Interface Method. This method eliminates the numerical artifacts, but it is rather complicated to implement. So I have created a simple, serial code to work out the details.

The results are encouraging. One popular benchmark is to model a block with an embedded sphere being sqeezed.

Inclusion setup

Unlike the falling block test, it has an analytic solution so we can know exactly how large our errors are. A side-by-side plot of the numerical and analytic solutions for the pressure shows no artifacts.


While most methods struggle to model viscosity ratios of $10^3$, I can model jumps of $10^10$ with little difficulty.