Spontaneous Rupture Propagation Research Articles

Spectral boundary integral equation method for simulation of 2D and 3D slip ruptures at bi‐material interfaces

AbstractThis paper presents a 3D spectral numerical scheme that can accurately model planarcracks of any shape residing on an interfacial plane between two elastic solids. The cracks propagate dynamically under the action of interfacial loads and interfacial constitutive laws. The scheme is a spectral form of the boundary integral equations (BIE) that relate the displacement discontinuity fields along the interfacial plane to the stress fields at the plane. In the BIE, the displacement discontinuities are expressed as a spatio‐temporal convolution of the traction components at the interface. This is in contrast with most previous studies which use a spatio‐temporal convolution of the displacement discontinuities or of the displacements at the interface. Due to the continuity of tractions across the interface, the present formulation is simpler, resulting in convolution kernels that can be written in closed‐form. Previous studies have relied on numerically obtained convolution kernels. The accuracy of the proposed scheme is validated with the known analytical solution of the 3D Lamb problem. Furthermore, new algorithms are developed for coupling the BIE with a slip‐weakening friction law, both for homogeneous materials as well as for bi‐material interfaces. This results in a widely applicable formulation for studying spontaneous rupture propagation. The algorithms are validated by comparing rupture simulation exercises to benchmark problems from Southern California Earthquake Center (SCEC). Finally, the methodology is used to simulate 3D frictional rupture propagation along a bi‐material interface.

Dynamic Rupture Simulations of the 2008 7.9 Wenchuan Earthquake: Implication for Heterogeneous Initial Stress and Complex Multifault Geometry

AbstractThe 2008 Mw7.9 Wenchuan earthquake ruptured northwest‐dipping imbricate oblique reverse faults along the Longmenshan thrust belt at the eastern margin of the Tibetan Plateau, and developed one of the most complex coseismic rupture patterns for reverse faulting events ever reported in intraplate settings. To evaluate to what extent complex fault geometry and heterogeneities in initial stress can explain the characteristics of coseismic observations, we simulate spontaneous dynamic rupture propagation governed by slip‐weakening friction law on the geometrically complex multifault system loaded by heterogeneous stress. Our model reproduces many observed features including the multifault ruptures, locations of the maximum slip, Interferometric Synthetic Aperture Radar observations, source time functions and peak ground velocities. Our results suggest that the maximum horizontal principal stress in the northern part of Beichuan Fault (BFC) is rotated approximately counterclockwise compared to that in the southern part of BFC, yet this rotation is not sufficient to produce a supershear rupture. We infer that the fault core zone in the southwestern Beichuan Fault (BCF) may have been more severely damaged than the northeastern counterpart, which provides an explanation that the rupture in the northern BCF propagates faster than in the southern portion. Besides, we find that Coulomb failure stress changes on Wenchuan‐Maoxian Fault (WMF) may be counteracted by the aftershock slips of the Lixian Fault, which suggests Wenchuan earthquake could not significantly push the WMF closer to failure.

Free‐Surface‐Induced Supershear Transition in 3‐D Simulations of Spontaneous Dynamic Rupture on Oblique Faults

AbstractSeveral observations and simulations suggested that the supershear rupture induced by free surface invariably occurs in strike‐slip faults rather than in dip‐slip faults. To explore the mechanism of supershear transition induced by free surface on oblique faults, we simulate spontaneous dynamic rupture propagation governed by slip‐weakening friction on nonvertical planar faults. We find that oblique thrust faults are more likely to produce free‐surface‐induced supershear transition than oblique normal faults and significant asymmetry arises between the two opposite directions along fault strike. The supershear rupture may only occur within shallow depth or over the entire seismogenic depth, depending on geometric relations between the slip and rupture propagation vectors. Our results help explain the lack of universal observation of supershear earthquakes. Besides, the asymmetric ground shaking on the oblique faults allows us to predict the possible spatial distribution of the high seismic risk.

The Community Code Verification Exercise for Simulating Sequences of Earthquakes and Aseismic Slip (SEAS)

Abstract Numerical simulations of sequences of earthquakes and aseismic slip (SEAS) have made great progress over past decades to address important questions in earthquake physics. However, significant challenges in SEAS modeling remain in resolving multiscale interactions between earthquake nucleation, dynamic rupture, and aseismic slip, and understanding physical factors controlling observables such as seismicity and ground deformation. The increasing complexity of SEAS modeling calls for extensive efforts to verify codes and advance these simulations with rigor, reproducibility, and broadened impact. In 2018, we initiated a community code-verification exercise for SEAS simulations, supported by the Southern California Earthquake Center. Here, we report the findings from our first two benchmark problems (BP1 and BP2), designed to verify different computational methods in solving a mathematically well-defined, basic faulting problem. We consider a 2D antiplane problem, with a 1D planar vertical strike-slip fault obeying rate-and-state friction, embedded in a 2D homogeneous, linear elastic half-space. Sequences of quasi-dynamic earthquakes with periodic occurrences (BP1) or bimodal sizes (BP2) and their interactions with aseismic slip are simulated. The comparison of results from 11 groups using different numerical methods show excellent agreements in long-term and coseismic fault behavior. In BP1, we found that truncated domain boundaries influence interseismic stressing, earthquake recurrence, and coseismic rupture, and that model agreement is only achieved with sufficiently large domain sizes. In BP2, we found that complexity of fault behavior depends on how well physical length scales related to spontaneous nucleation and rupture propagation are resolved. Poor numerical resolution can result in artificial complexity, impacting simulation results that are of potential interest for characterizing seismic hazard such as earthquake size distributions, moment release, and recurrence times. These results inform the development of more advanced SEAS models, contributing to our further understanding of earthquake system dynamics.

Recent Milestones in Unraveling the Full-Field Structure of Dynamic Shear Cracks and Fault Ruptures in Real-Time: From Photoelasticity to Ultrahigh-Speed Digital Image Correlation

AbstractThe last few decades have seen great achievements in dynamic fracture mechanics. Yet, it was not possible to experimentally quantify the full-field behavior of dynamic fractures, until very recently. Here, we review our recent work on the full-field quantification of the temporal evolution of dynamic shear ruptures. Our newly developed approach based on digital image correlation combined with ultrahigh-speed photography has revolutionized the capabilities of measuring highly transient phenomena and enabled addressing key questions of rupture dynamics. Recent milestones include the visualization of the complete displacement, particle velocity, strain, stress and strain rate fields near growing ruptures, capturing the evolution of dynamic friction during individual rupture growth, and the detailed study of rupture speed limits. For example, dynamic friction has been the biggest unknown controlling how frictional ruptures develop but it has been impossible, until now, to measure dynamic friction during spontaneous rupture propagation and to understand its dependence on other quantities. Our recent measurements allow, by simultaneously tracking tractions and sliding speeds on the rupturing interface, to disentangle its complex dependence on the slip, slip velocity, and on their history. In another application, we have uncovered new phenomena that could not be detected with previous methods, such as the formation of pressure shock fronts associated with “supersonic” propagation of shear ruptures in viscoelastic materials where the wave speeds are shown to depend strongly on the strain rate.

Physical mechanism for severe seismic hazard in the 2010 Yushu, China, earthquake (Mw= 6.9): insights from FEM simulations

The moderate-size 2010 Yushu earthquake (Mw 6.9) caused severe seismic damage around the source region which is situated on the sparsely populated hinterland of the Tibetan Plateau. But until now the mechanisms have not been well understood. To this end, we constructed the model to simulate the fault spontaneous rupture propagation, in which the realistic fault with a curved bend is imbedded. Our modeling results show that the special fault geometry controlled the rupture behavior, which encouraged rupture propagation speed transition from subshear to supershear with the speed of 5.17 km/s larger than the local shear wave velocity. Moreover, calculation results demonstrated that strong ground motion acceleration was greatly intensified by supershear ruptures, leading to widespread destruction. This may be the main reason why serious earthquake damage happened in the Yushu earthquake. In particular, we can see from numerous numerical experiments that rupture styles will be different if geometries of fault are varied. It is confirmed that it is the special geometry of the seismogenic fault of the Yushu mainshock produced grave seismic hazard. Thus, deeply investigating fault geometry will be helpful to better understand seismic source process and seismic hazard assessment. Key Points Curved bend in the seismogenic fault for the 2010 Yushu earthquake encouraged the supershear rupture transition. Strong ground motion acceleration was greatly intensified by supershear ruptures. Special geometry of the fault of the Yushu mainshock lead to seriouos seismic damage.

Modeling Megathrust Earthquakes Across Scales: One‐way Coupling From Geodynamics and Seismic Cycles to Dynamic Rupture

AbstractTaking the full complexity of subduction zones into account is important for realistic modeling and hazard assessment of subduction zone seismicity and associated tsunamis. Studying seismicity requires numerical methods that span a large range of spatial and temporal scales. We present the first coupled framework that resolves subduction dynamics over millions of years and earthquake dynamics down to fractions of a second. Using a two‐dimensional geodynamic seismic cycle (SC) model, we model 4 million years of subduction followed by cycles of spontaneous megathrust events. At the initiation of one such SC event, we export the self‐consistent fault and surface geometry, fault stress and strength, and heterogeneous material properties to a dynamic rupture (DR) model. Coupling leads to spontaneous dynamic rupture nucleation, propagation, and arrest with the same spatial characteristics as in the SC model. It also results in a similar material‐dependent stress drop, although dynamic slip is significantly larger. The DR event shows a high degree of complexity, featuring various rupture styles and speeds, precursory phases, and fault reactivation. Compared to a coupled model with homogeneous material properties, accounting for realistic lithological contrasts doubles the amount of maximum slip, introduces local pulse‐like rupture episodes, and relocates the peak slip from near the downdip limit of the seismogenic zone to the updip limit. When an SC splay fault is included in the DR model, the rupture prefers the splay over the shallow megathrust, although wave reflections do activate the megathrust afterward.

Dynamic weakening and amorphization in serpentinite during laboratory earthquakes

Abstract The mechanical properties of serpentinites are key factors in our understanding of the dynamics of earthquake ruptures in subduction zones, especially intermediate-depth earthquakes. Here, we performed shear rupture experiments on natural antigorite serpentinite, which showed that friction reaches near-zero values during spontaneous dynamic rupture propagation. Rapid coseismic slip (>1 m/s), although it occurs over short distances (<1 mm), induces significant overheating of microscale asperities along the sliding surface, sufficient to produce surface amorphization and likely some melting. Antigorite dehydration occurs in the fault walls, which leaves a partially amorphized material. The water generated potentially contributes to the production of a low-viscosity pressurized melt, explaining the near-zero dynamic friction levels observed in some events. The rapid and dramatic dynamic weakening in serpentinite might be a key process facilitating the propagation of earthquakes at intermediate depths in subduction zones.

Spectral formulation of the elastodynamic boundary integral equations for bi-material interfaces

A spectral formulation of the plane-strain boundary integral equations for an interface between dissimilar elastic solids is presented. The boundary integral equations can be written in two equivalent forms: (a) the tractions can be written as a space–time convolution of the displacement continuities at the interface (Budiansky and Rice, 1979) (b) the displacement discontinuities can be written as a space–time convolution of the tractions at the interface (Kostrov, 1966). Prior work on spectral formulation of the boundary integral equations has adopted the former as the starting point. The present work has for its basis the latter form based on a space–time convolution of the tractions. Tractions and displacement components are given a spectral representation in the spatial coordinate along the interface. The radiation damping term is then explicitly extracted to avoid singularities in the convolution kernels. With the spectral forms introduced, the space–time convolutions reduce to convolutions in time for each Fourier mode. Due to continuity of tractions at the interface, this leads to a simpler formulation and form of the convolution kernels in comparison to the formulation involving convolutions over the slip and opening history at the bi-material interface. The convolution kernels are validated by studying some model problems to which analytical solutions are known. When coupled with a cohesive law or a friction law at the interface, the formulation proposed here is of wide applicability for studying spontaneous rupture propagation.

On the initiation of sustained slip-weakening ruptures by localized stresses

Numerical simulations of dynamic earthquake rupture require an artificial initiation procedure, if they are not integrated in long-term earthquake cycle simulations. A widely applied procedure involves an ‘overstressed asperity’, a localized region stressed beyond the static frictional strength. The physical properties of the asperity (size, shape and overstress) may significantly impact rupture propagation. In particular, to induce a sustained rupture the asperity size needs to exceed a critical value. Although criteria for estimating the critical nucleation size under linear slip-weakening friction have been proposed for 2-D and 3-D problems based on simplifying assumptions, they do not provide general rules for designing 3-D numerical simulations. We conduct a parametric study to estimate parameters of the asperity that minimize numerical artefacts (e.g. changes of rupture shape and speed, artificial supershear transition, higher slip-rate amplitudes). We examine the critical size of square, circular and elliptical asperities as a function of asperity overstress and background (off-asperity) stress. For a given overstress, we find that asperity area controls rupture initiation while asperity shape is of lesser importance. The critical area obtained from our numerical results contrasts with published theoretical estimates when background stress is low. Therefore, we derive two new theoretical estimates of the critical size under low background stress while also accounting for overstress. Our numerical results suggest that setting the asperity overstress and area close to their critical values eliminates strong numerical artefacts even when the overstress is large. We also find that properly chosen asperity size or overstress may significantly shorten the duration of the initiation. Overall, our results provide guidelines for determining the size of the asperity and overstress to minimize the effects of the forced initiation on the subsequent spontaneous rupture propagation.

A domain decomposition approach to implementing fault slip in finite‐element models of quasi‐static and dynamic crustal deformation

We employ a domain decomposition approach with Lagrange multipliers to implement fault slip in a finite‐element code, PyLith, for use in both quasi‐static and dynamic crustal deformation applications. This integrated approach to solving both quasi‐static and dynamic simulations leverages common finite‐element data structures and implementations of various boundary conditions, discretization schemes, and bulk and fault rheologies. We have developed a custom preconditioner for the Lagrange multiplier portion of the system of equations that provides excellent scalability with problem size compared to conventional additive Schwarz methods. We demonstrate application of this approach using benchmarks for both quasi‐static viscoelastic deformation and dynamic spontaneous rupture propagation that verify the numerical implementation in PyLith.

Dynamic rupture of the 2011 Mw 9.0 Tohoku‐Oki earthquake: Roles of a possible subducting seamount

Using a hybrid MPI/OpenMP parallel finite element method for spontaneous rupture and seismic wave propagation simulations, we investigate features in rupture propagation, slip distribution, seismic radiation, and seafloor deformation of the 2011 Mw 9.0 Tohoku‐Oki earthquake to gain physical insights into the event. With simplified shallow dipping (10°) planar fault geometry, 1D velocity structure, and a slip‐weakening friction law, we primarily investigate initial stress and strength conditions that can produce rupture and seismic radiation characteristics of the event revealed by kinematic inversions, and seafloor displacements observed near the epicenter. By a large suite of numerical experiments aided by parallel computing on modern supercomputers, we find that a seamount of a dimension of ∼70 km by 23 km just updip of the hypocenter on the subducting plane, parameterized by higher static friction, lower pore fluid pressure, and higher initial stress than surrounding regions, may play a dominant role in the 2011 event. Its high strength stalls updip rupture for tens of seconds, and its high stress drop generates large slip. Its failure drives the rupture to propagate into the shallow portion that is likely velocity‐strengthening, resulting in significant slip near the trench within a limited area. However, the preferred model suggests that the largest slip in the event occurs near the hypocenter. High‐strength patches along the downdip portion of the subducting plane are most effective among several possible factors in generating high‐frequency seismic radiations, suggesting the initial strength distribution there is very heterogeneous.

Spontaneous Dynamic Rupture Propagation beyond Fault Discontinuities: Effect of Thermal Pressurization

Abstract A range of several fault segments often sequentially ruptures during an earthquake. We investigated the effects of thermal pressurization (TP) on dynamic rupture propagation beyond fault discontinuities by simulating spontaneous rupture propagation on two vertical strike‐slip fault segments. We revealed that a rupture can jump wider stepovers owing to TP, and that TP on a primary (nucleating) fault enables a rupture to jump at deep portions. In previous numerical studies on dry fault systems, it was found that a rupture sometimes fails to propagate to an unconnected fault, which is observed in the case of real earthquakes, and a rupture that successfully propagates is usually triggered near the surface of the Earth, unlike rupture evolution images obtained by seismic waveform modeling. TP can explain the inconsistencies between the previous numerical simulations and the observations, without depending on the heterogeneity of the initial stress and/or friction. Under depth‐dependent stress, we showed that TP enables a rupture to jump much wider stepovers at deep portions. If TP is in effect on faults, hydraulic diffusivity along with fault geometry can strongly control the characteristics of rupture propagation at fault discontinuities.

Three‐dimensional dynamic rupture simulation with a high‐order discontinuous Galerkin method on unstructured tetrahedral meshes

Accurate and efficient numerical methods to simulate dynamic earthquake rupture and wave propagation in complex media and complex fault geometries are needed to address fundamental questions in earthquake dynamics, to integrate seismic and geodetic data into emerging approaches for dynamic source inversion, and to generate realistic physics‐based earthquake scenarios for hazard assessment. Modeling of spontaneous earthquake rupture and seismic wave propagation by a high‐order discontinuous Galerkin (DG) method combined with an arbitrarily high‐order derivatives (ADER) time integration method was introduced in two dimensions by de la Puente et al. (2009). The ADER‐DG method enables high accuracy in space and time and discretization by unstructured meshes. Here we extend this method to three‐dimensional dynamic rupture problems. The high geometrical flexibility provided by the usage of tetrahedral elements and the lack of spurious mesh reflections in the ADER‐DG method allows the refinement of the mesh close to the fault to model the rupture dynamics adequately while concentrating computational resources only where needed. Moreover, ADER‐DG does not generate spurious high‐frequency perturbations on the fault and hence does not require artificial Kelvin‐Voigt damping. We verify our three‐dimensional implementation by comparing results of the SCEC TPV3 test problem with two well‐established numerical methods, finite differences, and spectral boundary integral. Furthermore, a convergence study is presented to demonstrate the systematic consistency of the method. To illustrate the capabilities of the high‐order accurate ADER‐DG scheme on unstructured meshes, we simulate an earthquake scenario, inspired by the 1992 Landers earthquake, that includes curved faults, fault branches, and surface topography.

MPI-OpenMP hybrid simulations using boundary integral equation and finite difference methods for earthquake dynamics and wave propagation: Application to the 2007 Niigata Chuetsu-Oki earthquake (Mw6.6)

Abstract For the past two decades, the combination between the boundary integral equation method (BIEM) and the finite difference method (FDM) has provided a powerful approach for simulating the dynamic rupture process along a complex fault system and the resultant seismic wave radiation and propagation for the purpose of understanding the physics of seismology and the generation of strong ground motion by large earthquakes. In the many-core architecture currently in use, it can be useful to take advantage of hybrid MPI-OpenMP programming. In this study, we are investigating the Mw6.6, 2007 Niigata Chuetsu-Oki, Japan, earthquake. In the BIEM simulation of the dynamic rupture process, MPI-OpenMP hybrid programming is 24 to 38% faster than MPI programming on the same resource using 8 to 16 nodes (with 8 cores). In the FDM simulation on the wave propagation, hybrid programming is 25 to 37% faster. This performance allows us to investigate this earthquake using various parameter settings. First we obtain the frictional parameters based on the proposed kinematic finite source model using BIEM, then model both the spontaneous dynamic rupture propagation and wave radiation by combining BIEM and FDM. It is found that on the given asperities, the peak stress can reach 20 MPa and reduced stress, -20 MPa. Dynamic rupture simulations require a barrier around the southernmost asperity so that this asperity is ruptured in a reverse manner. This effect can be observed very locally in the wave radiation. Furthermore, we observe that the dynamic rupture propagation between two cross-cutting segments generates a more complex wavefield. The change in fault geometry may be dominant over the asperities on a single plane for this earthquake.

3D dynamic simulations of spontaneous rupture propagation governed by different constitutive laws with rake rotation allowed

In this work we present a 3D Finite Difference numerical method to model the dynamic spontaneous propagation of an earthquake rupture on planar faults in an elastic half-space. We implement the Traction-at-Split-Nodes fault boundary condition for a system of faults, either vertical or oblique, using different constitutive laws. We can adopt both a slip-weakening law to prescribe the traction evolution within the breakdown zone or rate- and state-dependent friction laws, which involve the choice of an evolution relation for the state variable. Our numerical procedure allows the use of oblique and heterogeneous distribution of initial stress and allows the rake rotation. This implies that the two components of slip velocity and total dynamic traction are coupled together to satisfy, in norm, the adopted constitutive law. The simulations presented in this study show that the rupture acceleration to super-shear crack speeds occurs along the direction of the imposed initial stress; the rupture front velocity along the perpendicular direction is slower than that along the pre-stress direction. Depending on the position on the fault plane the orientation of instantaneous total dynamic traction can change with time with respect to the imposed initial stress direction. These temporal rake rotations depend on the amplitude of initial stress and on its distribution on the fault plane. They also depend on the curvature and direction of the rupture front with respect to the imposed initial stress direction: this explains why rake rotations are mostly located near the rupture front and within the cohesive zone.

Estimation of the Breakdown Slip from Strong-Motion Seismograms: Insights from Numerical Experiments

Fukuyama and Mikumo (2007) presented a method to extract informa- tion about the friction controlling earthquake rupture directly from near-field seismo- grams recorded out to several kilometers from a fault. This method computes a parameter (D 00) approximating the breakdown slip (Dc) from the rake-parallel dis- placement at the time of the peak particle velocity. The validation of the method was based on simple 2D steady-state Green's functions, an approach lacking sufficient physics of the rupture process to support their conclusions. Here, we use 3D simula- tions of subshear strike-slip spontaneous rupture propagation to demonstrate that D 00 is almost always controlled by rupture and wave propagation effects, rather than the fault friction process. Only if rupture reaches the Earth's surface and within a short distance from the fault (Rc) is it possible to extract information about Dc from strong-motion data, due to the fast decay with distance from the fault of the seismic energy related to the stress breakdown process. Beyond Rc, D 00 is controlled by the dynamic stress drop from the earthquake rupture without information about Dc. Rc is comparable to the length of the fault cohesive zone where the breakdown process takes place during rupture and approximately equal to 80% of the wavelength associated with the break- down frequency, defined as the reciprocal of the time span required by the stress to drop to the dynamic level. Our findings suggest that the Dc estimate by Fukuyama and Mikumo (2007) for the 2000 Mw 6.6 Tottori earthquake is unrelated to the breakdown slip. Moreover, the very narrow widths for Rc obtained from our simulations are in agreement with previous studies challenging the resolution of bandlimited strong- motion data to estimate Dc from kinematic-rupture models (Guatteri and Spudich, 2000; Spudich and Guatteri, 2004).

ShakeOut‐D: Ground motion estimates using an ensemble of large earthquakes on the southern San Andreas fault with spontaneous rupture propagation

We simulate ground motion in southern California from an ensemble of 7 spontaneous rupture models of large (Mw7.8) northwest‐propagating earthquakes on the southern San Andreas fault (ShakeOut‐D). Compared to long‐period spectral accelerations from the Next Generation Attenuation (NGA) empirical relations, ShakeOut‐D predicts similar average rock‐site values (i.e., within roughly their epistemic uncertainty), but significantly larger values in Los Angeles and Ventura basins due to wave‐guide focusing effects. The ShakeOut‐D ground motion predictions differ from those of a kinematically parameterized, geometrically similar, scenario rupture: (1) the kinematic rock‐site predictions depart significantly from the common distance‐attenuation trend of the NGA and ShakeOut‐D results and (2) ShakeOut‐D predictions of long‐period spectral acceleration within the basins of the greater Los Angeles area are lower by factors of 2–3 than the corresponding kinematic predictions. We attribute these differences to a less coherent wavefield excited by the complex rupture paths of the ShakeOut‐D sources.

Heterogeneous rupture on homogenous faults: Three‐dimensional spontaneous rupture simulations with thermal pressurization

To understand role of fluid on earthquake rupture processes, we investigated effects of thermal pressurization on spatial variation of dynamic rupture by computing spontaneous rupture propagation on a rectangular fault. We found thermal pressurization can cause heterogeneity of rupture even on a fault of uniform properties. On drained faults, tractions drop linearly with increasing slip in the same way everywhere. However, by changing the drained condition to an undrained one, the slip‐weakening curves become non‐linear and depend on locations on faults with small shear zone thickness w, and the dynamic frictional stresses vary spatially and temporally. Consequently, the super‐shear transition fault length decreases for small w, and the final slip distribution can have some peaks regardless of w, especially on undrained faults. These effects should be taken into account of determining dynamic rupture parameters and modeling earthquake cycles when the presence of fluid is suggested in the source regions.

A 3-D hybrid finite-difference-finite-element viscoelastic modelling of seismic wave motion

SUMMARY We have developed a new hybrid numerical method for 3-D viscoelastic modelling of seismic wave propagation and earthquake motion in heterogeneous media. The method is based on a combination of the fourth-order velocity–stress staggered-grid finite-difference (FD) scheme, that covers a major part of a computational domain, with the second-order finite-element (FE) method which can be applied to one or several relatively small subdomains. The FD and FE parts causally communicate at each time level in the FD–FE transition zone consisting of the FE Dirichlet boundary, FD–FE averaging zone and FD Dirichlet zone. The implemented FE formulation makes use of the concept of the global restoring-force vector which significantly reduces memory requirements compared to the standard formulation based on the global stiffness matrix. The realistic attenuation in the whole medium is incorporated using the rheology of the generalized Maxwell body in a definition equivalent to the generalized Zener body. The FE subdomains can comprise extended kinematic or dynamic models of the earthquake source or the free-surface topography. The kinematic source can be simulated using the body-force term in the equation of motion. The traction-at-split-node method is implemented in the FE method for simulation of the spontaneous rupture propagation. The hybrid method can be applied to a variety of problems related to the numerical modelling of earthquake ground motion in structurally complex media and source dynamics.