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.