A multigrid solver for modeling complex interseismic stress fields

Abstract

We develop a multigrid, multiple time stepping scheme to reduce computational efforts for calculating complex stress interactions in a strike-slip 2D planar fault for the simulation of seismicity. The key elements of the multilevel solver are separation of length scale, grid-coarsening, and hierarchy. In this study the complex stress interactions are split into two parts: the first with a small contribution is computed on a coarse level, and the rest for strong interactions is on a fine level. This partition leads to a significant reduction of the number of computations. The reduction of complexity is even enhanced by combining the multigrid with multiple time stepping. Computational efficiency is enhanced by a factor of 10 while retaining a reasonable accuracy, compared to the original full matrix-vortex multiplication. The accuracy of solution and computational efficiency depend on a given cut-off radius that splits multiplications into the two parts. The multigrid scheme is constructed in such a way that it conserves stress in the entire half-space.

Introduction

Multiplications of a vector by a dense matrix demand high computational expense for half-space elastodynamic solutions in a fault model for the simulation of seismicity (Ben-Zion and Rice, 1993, Ben-Zion, 1996, Zöller et al., 2004, Zöller et al., 2005) Such fault models calculate the evolution of slip, stress, and other related quantities as a response on long-term accumulation of stress in the Earth's crust resulting from the motion of the tectonic plates. The outcome of those models is an earthquake catalog including time, hypocenter and magnitude of each event. These data are useful for purposes of earthquake statistics (e.g. frequency-size distributions, recurrence times) and seismic hazard studies. In this way, the poor statistics of observational earthquake data can be overcome to some extent. The interseismic build-up of stress period is related to plate motion with constant velocity. The release of stress comes from power-law creep (interseismic) accounting for aseismic processes and from earthquakes (coseismic). On average there is a balance between build-up and release of stress (backslip model). Recurrence times of large earthquakes are tens to hundreds of years, while the earthquake itself occurs on a time-scale of a few seconds. In simple fault models, different regimes of a fault are loaded independently during the periods between earthquakes. More realistic models include complex spatio-temporal interactions at each time. This leads to expensive multiplications as in many-body simulations in other physical problems. There have been a number of efforts to reduce the computational cost for many-body calculations such as the Barnes–Hut method, the parallel tree methods, and the fast multipole expansion method in tree algorithms for long-range potentials. Meanwhile mesh-based fast algorithms include the particle–particle particle–mesh method, the particle–mesh Ewald summation, and multigrid methods and adaptive refinement (Griebel et al., 2007 and references therein). These methods attempted to reduce the complexity of N×N to $\mathcal{O}(N)$ or $\mathcal{O}(N\log N)$ for interactions of a set of N particles or bodies. Moreover, these approaches can be combined with a multiple time stepping to further enhance the computational speed. To develop efficient solvers for elastodynamics some studies use, for example, fast multipole boundary element methods (e.g., Chaillat et al., 2008) or parallel computations on grids with different spacings (Aoi and Fujiwara, 1999).

Iterative multigrid methods are used as fast numerical methods for the solution of a linear equation arising for partial differential equations (Trottenberg et al., 2007). The idea of using multiple grids was adopted for an efficient multiplication by a dense matrix, in a non-iterative way (e.g., Brandt and Lubrecht, 1990). Amongst fast multiplication methods, multigrid methods have advantages in that they are relatively simple to implement and applicable to general potentials. In this work we test multigrid methods combined with a multiple time stepping, based on the multigrid approach in Skeel et al. (2002), hereafter STH02. Then we compare the results with those from the original full matrix-vector multiplications. In preliminary tests we found that the method of STH02 can lead to a lower error than using the multigrid method suggested by Brandt and Lubrecht (1990), hereafter BL90, for the problem in the fault model. The difference between the algorithms in BL90 and STH02 will be discussed in the section of algorithm description.

In the following section, we describe the fault in earthquake modeling briefly and the kernel used in the multiplications. In Section 3, the idea of the multigrid algorithm in STH02 is explained and the advantage of this approach over that in BL90. In Sections 4 and 5, we present the results from the multigrid multiplications in the fault model. Finally we summarize this study and deliver outlook for implementations in operational models.