Abstract

Numerical simulations of seismic cycles with rate-, state-, and temperature-dependent friction explain the various modes of seismic and aseismic ruptures in the brittle section of the lithosphere. However, the effects of viscoelastic flow in the ductile layers remain challenging to incorporate due to the wide range of length scales involved, from extremely localized within fault zones to widely distributed in the lower crust and asthenosphere. Here, we describe simulations of seismic cycles in a viscoelastic half-space using the integral method that combines discrete surface and volume elements to capture the coupling between brittle and ductile deformation. Viscoelastic flow is captured by cuboidal and tetrahedral volume elements within rectilinear and curvilinear meshes, respectively. The model resolves all phases of the seismic cycle under the radiation-damping approximation, including the nucleation and propagation of earthquake ruptures, but also the viscoelastic relaxation that follows in the ductile layers. We illustrate the approach in three dimensions with numerical simulations of seismic cycles on finite strike-slip and thrust faults overlying a viscoelastic lower crust with linear and nonlinear rheology. In two-dimensional models of subduction zones with the in-plane strain approximation, the ductile regions are meshed with triangle volume elements. The use of Green's functions only requires the discretization of the actively deforming region, resulting in a relatively small mesh. We provide open-source software implementing the method with parallel computing in a distributed architecture. The approach allows increasingly realistic representations of the lithosphere-asthenosphere system with nonlinear constitutive laws in structurally complex tectonic settings.