SEISMIC WAVE PROPAGATION WITH THE ADER-DG METHOD 

The ADER-DG (Arbitrary high-order DErivative Riemann Discontinuous Galerkin) method is essentially a finite element method which makes use of numerical fluxes to solve discontinuities of the variables at element interfaces. It has been applied to seismological problems (Käser & Dumbser 2006, Dumbser & Käser 2006), where it is combined with triangular or tetrahedral meshes, thus making it easier to describe complex geological structures.

Some outline of the capabilities of this method is described here: 

Triangular or tetrahedral cell decomposition: The space is discretized in a conforming mesh of triangles (or tetrahedra in 3D). This makes the method well suited for handling complex geometrical features (topography, material discontinuities...) while it is very fast and easy to set a mesh up.
Arbitrary high-order basis functions set: A complete basis of functions is used to describe the unknowns (stresses and velocities) inside each element. An arbitrary polynomial degree can be used to increase the resolution of the method.
Information is shared through „fluxes“: The functions are not forced to be continuous between the cells. Information is shared between neighbouring cells using "numerical fluxes", a generalisation of Gauss' Theorem which is widely used in Finite Volume methods.
Single-step high-order time integration: The accuracy of the scheme can be higher than in conventional methods due to a time-integration scheme which has the same order as the space-integration. This scheme does not require the storage of intermediate substeps as, for example, in Runge-Kutta time integration schemes.
Local time stepping /p-adaptivity:The polynomial order and the time step of each cell can be set individually for optimal resolution and efficiency. CFL stability conditions become local.

Mesh model of a reservoir inspired by the SEG/EAGE overthrust model. The color scale represents different seismic wave speeds. Tetrahedral meshes can honor all discontinuities and allow for mesh refinement there were the wavelengths are expected to be small or at regions of interest.

Gather of synthetic seismograms recorded at a receiver line along the y-axis of the model above. The line crosses the source point.	Gather of synthetic seismograms recorded at a receiver line along the y-axis of the model above, running parallel to the gather shown on the left at a distance of 1000m from it.

RHEOLOGICAL MODELS

The assumption that the Earth is completely elastic and isotropic has lead to many important findings in seismology, as are for example the existence of a liquid inner core or the location of most major discontinuities. However, at certain scales, the Earth shows a slightly different behavior. In order to account for more complex Earth models, three major rheologies have been included in the ADER-DG method: viscoelasticity, anisotropy and poroelasticity, as well as all possible mixed combinations of the three.

Viscoelasticity is handled by the use of a Generalized Maxwell Body approach, in which a combination of springs and dashpots model the physics of viscoelastic wave propagation with its two main effects: attenuation of the amplitudes and wave dispersion. The increase in size of the number of equations and variables to be solved had to be handled by separating the Jacobian and reaction matrices into matrix blocks which present similar structures, thus allowing for a remarkable saving in the amount of operations to be done by the scheme at every time step. The resulting scheme is very efficient and has shown excellent performance in comparisons with other state-of-the-art approaches.

Snapshot of an isotropic-anisotropic material interface

Anisotropy is a very complex problem when using unstructured grids. Possible symmetries of the medium most likely will not coincide with the structure of the mesh and therefore we cannot make use of them. In addition an accurate representation of the flux through a surface with an arbitrary orientation poses a problem not been handled before. We have successfully implemented the general most tricilinic anisotropy with both exact and approximate fluxes and have shown that the order of the scheme is unaffected by the increase of complexity of the physics of the problem. The coupling of anisotropic and viscoelastic effects has been carefully derived and implemented in order to get a consistent and very general solution for cases in which both rheologies are important.

The poroelasticity theory established by Biot in the early 40s describes the mechanical properties of fluid-filled porous materials. These materials appear often in nature, especially in reservoir, and the most characteristic effect they show in wave propagation is the apparition of a pressure wave of the second-kind, also known as slow P-wave, in addition to the classic elastic waves. The nature of this wave comes from out-of-phase solid and fluid particle motion. In general P-waves of the second kind travel even slower than S-waves. Another further new aspect is the fact that Biot's equations, in the low-frequency range (often below the KHz, in practice), the slow wave becomes a diffusive mode, only noticeable close to the source or material heterogeneities. This effect is very difficult to model with explicit time-integration numerical methods due to unstability problems and therefore has been rarely modelled.

A classical approach to describe the dynamic response of porous materials is to use elastic- and viscoelastic-equivalent models, so that elastic theory is used but adjusting the material parameters in order to produce waves which propagate at the same velocities as those observed in the porous rocks.

Diffusion-type wave at source location in poroelastic materials at seismic frequencies

This however is only valid to model homogeneous rocks and, additionally, fails to capture the energy balance present at porous/porous and porous/elastic interfaces. We have combined the DG method with a space-time discontinuous time integration to build up a scheme (ADER-DG(ST)) able to accurately model such effects without the strong computational load required by globally implicit methods. In ADER-DG(ST) schemes the time integration is also handled by numerical fluxes, producing a method which captures the timescales of the wave and diffusion phenomena simultaneously. High-order convergence properties are achieved and, furthermore, the solution is correct in the assymptotic limit. The possibility of solving Biot's equations in the whole frequency range combined with the use of tetrahedral meshes might allow to improve modelling of reservoir fields and to further understand the seismic signature of poroelastic materials in observed data

This study of these rheologies and their effects on the sesimic wavefield have led to my Ph.D. thesis, available here.

FAULT DYNAMICS

A further interesting problem in seismology is understanding the mechanical and thermodynamical processes involved in the slip of faults during earthquakes. Unfortunately, the conditions present at seismogenic faults are difficult to reproduce in laboratories. Therefore, numerical modelling is one of the most useful tools to explore the validity of assumptions regarding the fault dynamics at seismological time scales.

The ADER-DG method offers an alternative methodology for the modelling of the rupture process and the wave field generated by it. The discontinuous representation of the variables at element interfaces makes it simple to represent dislocation phenomena and the refinement / coarsening possibilities of tetrahedral meshes allows to capture small scale phenomena at the fault itself. A further advantage is that arbitrary fault geometries can be used as long as they are honored by the computational mesh.

Example of mesh honoring a curved fault. A strong mesh refinement has been applied to better resolve the fault geometry.	Wave field generated by the rupture of the fault. The effect of the fault tips, which radiate energy as the rupture stops suddenly, can be clearly observed.