Dynamic earthquake rupture modelled with an unstructured 3-D spectral element method applied to the 2011 M9 Tohoku earthquake

An important goal of computational seismology is to simulate dynamic earthquake rupture and strong ground motion in realistic models that include crustal heterogeneities and complex fault geometries. To accomplish this, we incorporate dynamic rupture modelling capabilities in a spectral element solver on unstructured meshes, the 3-D open source code SPECFEM3D, and employ state-of-the-art software for the generation of unstructured meshes of hexahedral elements. These tools provide high flexibility in representing fault systems with complex geometries, including faults with branches and non-planar faults. The domain size is extended with progressive mesh coarsening to maintain an accurate resolution of the static field. Our implementation of dynamic rupture does not affect the parallel scalability of the code. We verify our implementation by comparing our results to those of two finite element codes on benchmark problems including branched faults. Finally, we present a preliminary dynamic rupture model of the 2011 M_w 9.0 Tohoku earthquake including a non-planar plate interface with heterogeneous frictional properties and initial stresses. Our simulation reproduces qualitatively the depth-dependent frequency content of the source and the large slip close to the trench observed for this earthquake.


I N T RO D U C T I O N
3-D numerical methods for earthquake rupture dynamics and ground motion simulation capable of incorporating complex nonplanar fault systems, rough surface topography, non-linear rheologies and the heterogeneous structure of the Earth interior (e.g. Ma et al. 2007;Barall 2009;Ely et al. 2009Ely et al. , 2010Tago et al. 2012;Pelties et al. 2012) are gaining increasing importance in the study of the physics of earthquakes. Because rupture dynamics involves small-scale processes that need to be accurately resolved (e.g. Day et al. 2005), such numerical simulations pose high demands in terms of memory and running time and need parallel computation on thousands of processors to achieve an accurate numerical solution of the dynamic rupture process. Today's supercomputer resources are allowing earthquake scientists to use such numerical models to investigate the physics of earthquakes at high resolution and large scales that were previously beyond hardware capabilities (Olsen et al. 2009). This new era of large-scale highresolution 3-D numerical calculations allows to unveil new features of the rupture propagation, contributing to a better understanding of the mechanics and physics of earthquakes, which in turn provides * Now at: SwissNuclear, Aarauerstrasse 55, 4601 Olten, Switzerland. useful insights for improving our capability to predict ground motion for assessment of seismic hazard.
Dynamic earthquake models usually idealize the rupture process as a dynamically running shear crack on a frictional interface embedded in an elastic continuum. The spatio-temporal evolution of stress and slip during fault rupture is determined by solving the elastodynamic equation coupled to frictional sliding, leading to a highly non-linear mixed boundary value problem (e.g. Andrews 1976; Das & Aki 1977;Day 1982). These dynamic models have been implemented in several volumetric 3-D numerical algorithms based on finite difference methods (FDMs), the different classes of finite element methods (FEMs) and finite volume methods (e.g. Dalguer & Day 2006Kaneko et al. 2008;Dalguer 2012;Tago et al. 2012;Pelties et al. 2012;Kozdon & Dunham 2013, and references therein). Standard FDM, though widely used for wave propagation, are limited to planar faults and face serious difficulties to be extended to complex fault geometries. A notable exception in 3-D is the FDM of Ely et al. (2009Ely et al. ( , 2010 that uses different operators with irregular geometries. In contrast, the boundary element method (BEM) for dynamic rupture problems can handle non-planar faults and is highly accurate. Aochi & Fukuyama (2002) applied this methodology to the 1992 Landers earthquake taking into account a complex fault system in homogenous media. However, BEM is not suitable for heterogeneous velocity structures or inelastic media, and few implementations account for the free surface. FEM overcomes these difficulties naturally owing to its capability to mesh general geometries in heterogeneous and non-linear media. However, traditional, low-order FEM with mass lumping produces dispersion that penalizes the accuracy of wave propagation and dynamic stress transfer in rupture simulations. High-order FEMs, such as spectral element methods (SEMs; e.g. Festa & Vilotte 2005;Kaneko et al. 2008), and discontinuous Galerkin methods (DGMs; e.g. Tago et al. 2012;Pelties et al. 2012), with diagonal mass matrices by construction, are very accurate and maintain the geometrical flexibility.
Here we build upon the unstructured 3-D open source spectral element code SPECFEM3D (http://www.geodynamics. org/cig/software/specfem3d), a highly scalable Fortran 90 code parallelized using the Message Passing Interface (MPI) for large-scale simulations. This tool was introduced in seismology as a solver for the elastic wave equation by Komatitsch & Vilotte (1998) and Komatitsch & Tromp (1999, 2002. The main characteristic of SEM, compared to the standard FEM, is that it uses high-order basis functions that make the method accurate enough to solve the wave equations only with four to five nodes per wavelength in most practical situations (Komatitsch & Tromp 1999). The mass matrix is naturally diagonalized by using the Gauss-Lobatto-Legendre (GLL) nodes inside the elements for both quadrature and interpolation, while preserving the accuracy (De Basabe & Sen 2010). In the current version of SPECFEM3D (Peter et al. 2011), complex geometries are handled with unstructured meshes of hexahedral elements, which can be constructed using a mesh generation tool such as CUBIT (http://cubit.sandia.gov/). The code can also take advantage of GPU-enhanced architectures (Rietmann et al. 2012). SEM can handle non-linear bulk rheologies, such as continuum damage (e.g. Ampuero et al. 2008;Lyakhovsky et al. 2009) and plasticity (Kaneko & Fialko 2011;Xu et al. 2012a,b;Gabriel et al. 2013). The accuracy of SEM for dynamic rupture simulations with off-fault plastic deformation on planar faults and structured grids has been successfully tested in the Southern California Earthquake Center (SCEC) benchmarks TPV13 and TPV27 (Harris et al. 2011;Gabriel et al. 2013).
Here we present the implementation of the boundary conditions for spontaneous dynamic rupture into SPECFEM3D. Our implementation follows the principles introduced by Kaneko et al. (2008) and involves encapsulated modules plugged into the SPECFEM3D code. It provides the capability to model dynamic rupture for multiple, non-planar faults governed by slip-weakening friction and rate-and-state friction. We verify the efficiency and accuracy of our implementation. We show that the parallel computation is scalable to thousands of processors, enabling high-performance execution for large-scale dynamic rupture calculations. The accuracy of the code is successfully verified through benchmark problems developed by the SCEC/USGS Dynamic Rupture Code Validation Project (Harris et al. 2009), including 3-D problems with branched faults. We finally apply our new tool to develop a preliminary dynamic model of the rupture process of the 2011 M w 9.0 Tohoku earthquake incorporating a 3-D non-planar geometry of the megathrust interface.

Statement of the problem
We consider the problem of spontaneous earthquake rupture propagation on a pre-existing fault surface embedded in an elastic medium enclosed by a surface ∂ . The evolution of slip is controlled by a friction law and the initial stresses on the fault, see Fig. 1. The problem is governed by the linear elastodynamic equations: where ρ is the density of the medium, σ stress tensor and u the incremental displacement field. We assume linear elasticity and small displacements (Hooke's law): where c is the elastic tensor and ε the strain defined as (∇u + ∇u T )/2. We also assume zero initial conditions on displacements and velocities and free stress boundary conditions at the surface of the Earth (σ · n = 0, where n is the vector normal to the free surface). In practice, the model domain is truncated to a finite size and approximate absorbing boundary conditions are applied on the artificial exterior boundaries (Komatitsch & Tromp 1999). However, to simplify this presentation we will ignore the absorbing boundary conditions. We treat the fault as a surface of displacement discontinuity. We represent the fault as a 2-D interface composed of two matching surfaces in contact, = + − (Fig. 1). The slip is defined as the displacement discontinuity across the fault, where u + and u − denote the displacements on + and − , respectively. The traction on the fault surface − is denoted by: where n is the normal vector of − , pointing towards + (see Fig. 1). To simplify this presentation, we ignore the possibility of fault opening: s · n = 0. We denote by T T and T N the tangential and normal tractions on − . The normal traction is negative in compression. The friction boundary conditions on the fault interface are: whereṡ is the slip velocity vector and μ is the friction coefficient, which can depend on slip, slip rate and other fault state variables.
Here we adopt the linear slip weakening friction law (e.g. Ida 1973;Palmer & Rice 1973;Andrews 1976): where μ s and μ d are the static and dynamic friction coefficients, respectively, D c the critical slip distance and δ andδ are the magnitude of the slip and slip rate, respectively. Despite its simplicity, this friction law represents key features of fault strength: a finite friction coefficient μ, progressive weakening (μ s − μ d ) and finite fracture energy Just like FDM and FEM, also SEM suffers from non-physical oscillations at the high frequencies that are not well resolved by the Figure 1. Representation of a fault as two split surfaces. The fault interface, , is embedded in an elastic medium, , and is composed by two matching surfaces, ± , that can deform independently. For clarity, the two surfaces are plotted as separated in the zoomed in view, although they are most typically considered to be in frictional contact. The vector normal to the surface − , pointing towards + , is denoted byn. On each side of the fault tractions are denoted by τ ± and displacements by u ± . mesh. To reduce these spurious oscillations, we introduce Kelvin-Voigt damping (Day & Ely 2002) in a narrow layer of elements surrounding the fault. This amounts to replace eq. (2) by where η is a viscous relaxation time.

Variational formulation
Like in the FEM, the SEM discretizes the weak (variational) form of the governing eq. (1) by doting it with an arbitrary test vector w and integrating over a finite volume . After integrating by parts and applying the free surface boundary condition on ∂ , we get : The fault is viewed as an infinitely thin closed hole, a slit, whose surface is naturally portioned into two surfaces in contact, = + − . The solution and test functions are described by smooth fields inside the domain with a slit . This naturally allows for a displacement discontinuity across the fault . The left-hand side of eq. (12) can be decomposed over the two faces in contact, where the traction in each surface satisfies T − = −T + . Taking as a reference the − fault side, we define T = T − and obtain: where w = w + − w − is the difference of the test function across the fault. Finally, replacing (14) in (12) we get the elastodynamic equation with the fault term included: The problem is to find u that satisfies (15) for all w together with Hooke's law (2) or the Kelvin-Voigt constitutive eq. (11), the friction conditions (5)-(7), the friction law (8) and (9) and the given initial conditions.

Discrete formulation
The discretization of the weak form of the elastodynamic eq. (15) by the SEM leads to the matrix equation: where M and K are the mass and stiffness matrix, respectively, given by Komatitsch & Tromp (2002), C is a viscosity matrix and B is the fault boundary matrix given by Kaneko et al. (2008). For Kelvin-Voigt viscosity we set C = K H, where H is a diagonal matrix whose only non-zero components are the diagonal terms associated with the nodes of the near-fault viscous elements and are equal to the viscosity η. The relative fault traction is τ = T − T o , where T o is the initial stress on the fault. The time discretization is done with a central explicit Newmark algorithm: where u n ,u n andü n represent the particle displacement, velocity and acceleration at the nth time step, respectively, and t is the time step size. To preserve the explicit nature of the time stepping scheme, the viscous term in eq. (18) involves a semi-updated velocity defined aṡ To update the values of fault traction, we manipulate eqs (18) and (19) on the fault split nodes to obtain the following expression: where u =u + −u − is the slip rate, is the fault impedance matrix, is the 'stick traction' that would prevail if the fault node suddenly arrested anḋ Subscripts ± denote values on the nodes lying on one of the two sides of the fault, ± . Here we do not allow for fault opening, hence the slip velocity normal to the fault vanishes, u N n+1 = 0, and the normal fault traction remains as: The explicit Newmark algorithm (17) readily provides an update of displacement and slip, with which we update the friction coefficient according to (8) and (9). To update the shear fault traction T T we solve (21) together with the friction conditions (5)- (7). This can be efficiently done on a node-by-node basis, because the matrix Z is diagonal. The solution is: We then compute the relative stress on the fault, τ = T − T o , reinsert it in eq. (18) to update accelerations and finally update velocities (19). Efficient damping of spurious oscillations is achieved by a thin layer of Kelvin-Voigt elements surrounding the fault, with thickness of 1 or 2 elements on each side of the fault. Kelvin-Voigt viscosity introduces a frequency-dependent quality factor Q −1 ( f ) = 2πηf, where f is frequency. We find that the longest period of spurious oscillations is several times the critical time step t fault that a cubic spectral element with same edge length as a fault surface element would have in a medium without viscosity. The viscous time η is set uniformly within the layer to a fraction (10 to 30 per cent) of t fault to achieve significant attenuation of spurious oscillations after a few cycles or less. Larger values damp high frequencies more aggressively but also affect lower frequencies, rupture speed and peak slip velocities. Viscosity modifies the numerical stability of the Newmark time stepping algorithm: the critical time step in a simulation with Kelvin-Voigt damping, t KV c , is smaller than that in a purely elastic simulation, t c . It can be shown, from eq. (9.4.28) of Hughes (2000), that Geometries with narrow angles lead to very small and highly deformed elements, for which the time step needs to be decreased significantly to satisfy the stability condition. This inherent limitation of the explicit time stepping algorithm usually adopted in SEM and FEM can certainly lead to computationally impractical scenarios. This is aggravated by the more stringent stability conditions imposed by the Kelvin-Voigt viscosity included near the fault. Indeed, highly deformed elements tend to require t c t fault , which may lead to t c η and, from eq. (27), to t KV c / t c ∼ t c /2η 1. Strategies to potentially mitigate this issue include local time stepping (Peter et al. 2013), implicit time integration in the critical elements, usage of tetrahedral and prismatic spectral elements in critical regions of the model (Komatitsch et al. 2001;Pasquetti & Rapetti 2006) and optimal damping by flux-based solver procedures borrowed from the DGM (De la Puente et al. 2009;Pelties et al. 2012).

PA R A L L E L I Z AT I O N A N D S C A L A B I L I T Y
In dynamic rupture simulations the computational mesh needs to be dense enough to resolve the breakdown zone at the rupture front, whose size is controlled by the rupture speed, frictional strength drop and slip-weakening distance (Day et al. 2005). The simulation of large earthquakes typically requires a node spacing less than a few hundred metres (Harris et al. 2009). For a total domain size of a few hundred kilometres and a number of GLL nodes per element edge NGLL = 5, the total number of spectral elements needed is on the order of tens of millions. To satisfy the contact and friction conditions the elements carrying the fault interface need to be treated differently than the elements in the rest of the bulk. One approach is to assign during domain decomposition all the spectral elements that are in contact with fault surfaces to a single processor. We initially adopted this strategy (as did Kaneko et al. 2008) for the simplicity of its implementation. However, for large simulations this approach leads to a major load imbalance, with a bottleneck waiting for the processor that contains all the fault elements.
To achieve load balancing, we parallelized the fault solver as well. During domain decomposition, we assigned matching pairs of elements on both sides of the fault to the same processor, the one with lowest rank of the pair. This simplifies the implementation and avoids solver communications across the fault. The fault normal vector (n) and fault boundary matrix (B) were pre-assembled across MPI interfaces along the fault, and internal forces are globally assembled before passing them to the fault solver. Hence, no additional assembly operation (no additional interprocessor communication) is performed by the fault solver. This strategy is expected to generate a minimal impact on the overall cost of computations, which should remain dominated by the bulk wave propagation solver. The original SPECFEM3D code has been shown to have good scaling properties for wave propagation problems (Komatitsch et al. 2009). We demonstrate here that our implementation of fault dynamics does not affect its parallel scalability.
We illustrate the strong scaling of the code in the communitybased SCEC dynamic rupture benchmark problem TPV5. The problem comprises a fault 30 km long and 15 km deep. We placed absorbing boundaries 15 km away along strike from the fault tips, 25 km below the bottom edge of the fault and 30 km away in the fault-normal direction. We adopted a spectral element size of 400 m with 5 GLL nodes per element edge, corresponding to the maximum recommended average grid size of 100 m (Harris et al. 2009). This resulted in 2 265 000 spectral elements. We ran the simulation at the Swiss National Supercomputing Center (CSCS) on Monte Rosa, a Cray XE6 system with 1496 compute nodes consisting of two 16-core AMD Opteron 6272 2.1 GHz CPUs and 32 GB of memory, and with high-performance networking through a Gemini 3-D torus interconnect. The theoretical peak performance of Rosa is 402 Tflops. We choose numbers of processors in powers of 2 ranging from 64 to 8192. We suppressed intermediate outputs, as our focus was on verifying the scaling of the combined bulk-fault solver. Fig. 2(a) shows the total wall clock time taken by the solver (bulk and fault) to complete one TPV5 simulation. The 'ideal' line corresponds to the perfect scaling: wall clock time inversely proportional to number of processors. The code scales well within the range of number of processors we tested. We also tested the scaling of the original SPECFEM3D code without fault implementation. For this purpose, we considered the same domain size and element size as that of our TPV5 simulations, but without the split-node fault surface, and we prescribed an explosion point source at the centre of the domain. We repeated the scalability test in the same system and for the same set of processors as those previously used. The results, shown in Fig. 2, demonstrate that the fault solver does not cause any significant load imbalance and does not affect the overall performance of the code.
While TPV5 was used to analyse strong scaling, for weak scaling we consider a different version of the same benchmark, TPV205. Essentially, a TPV5 (100 m grid size) run on 256 processors is compared with the same benchmark problem solved with 200 m grid size on 16 processors and 50 m grid size on 4096 processors. These three sets of simulations have, in principle, the same load per processor: the total number of operations for fixed domain size and duration scales as 1 per (grid size) 4 . Fig. 2(b) shows the weak scalability results for SPECFEM3D with our fault implementation. Wall clock time is normalized with respect to that of the 50 m grid size simulation. The weak scalability is overall satisfactory. The minor (2 per cent) deviation in weak scaling could be attributed to the fact that the number of spectral elements on the fault plane in the finest resolution mesh is not exactly four times that in the lowest resolution mesh, due to constraints imposed by the fixed fault dimensions.

N U M E R I C A L T E S T S : A S S E S S M E N T O F N U M E R I C A L S O L U T I O N S
To verify our implementation of the dynamic rupture boundary conditions in SPECFEM3D, we have reproduced several 3-D test problems from the SCEC dynamic rupture code validation project (e.g. Harris et al. 2009Harris et al. , 2011) and compared our results to those of other published methods. Here we report only our verification results for the test problems TPV24 and TPV25, which are representative of the non-planar fault geometries that our method can handle. We first summarize the setting of these test problems, then we present our results and qualitative comparisons to other methods. The accuracy of the method in a thrust fault model with shallow dip angle is demonstrated in the Appendix and summarized at the end of Section 5.2.

Description of the rupture problem on a branched fault
The SCEC test problems TPV24 and TPV25 consist of spontaneous rupture propagation on a branched fault system comprising two segments, a main fault and a branch fault, embedded in a uniform elastic isotropic half-space (Fig. 3). The two fault segments are vertical, planar, strike-slip faults that reach the Earth's surface. In TPV24 the faults are right-lateral, while in TPV25 they are leftlateral. The main fault is 28 km long and 15 km deep, and the branch fault is 12 km long and 15 km deep. The branch fault splays off the main fault at an angle of 30 • , at 12 km from the right edge of the main fault. It is assumed that the slip on the branch fault tapers smoothly to zero at the junction with the main fault. The S-wave velocity is 3463 m s −1 , the P-wave velocity is 6000 m s −1 and the density is 2670 m s −1 . The hypocentre is located on the main fault at 8 km to the left of the junction point and at 10 km depth. Rupture nucleation is achieved by prescribing time-weakening over a region that grows with smoothly variable rupture speed.
We employ a semi-spherical mesh with gradual increase of the element size as a function of radial distance (Fig. 3b). Coarsening is an efficient approach to increase the domain size, which improves the accuracy of the static field. The spherical shape of the outer boundary allows the angle of incidence of waves on the absorbing boundaries to be closer to normal, which reduces spurious reflections. A smooth mesh coarsening is needed in the region surrounding the fault to avoid artificial wave reflections.
The time step in our simulations is t = 0.5 ms. To attenuate the spurious high-frequency oscillations, we set the Kelvin-Voigt viscosity as η = 0.3 t fault on the main fault and η = 0.2 t fault on the branch, where t fault = 5.7 ms is the elastic critical time step corresponding to the element size on the fault.

Results and comparison to other numerical methods
The complexity of rupture path selection in branched fault systems has been previously studied by Bhat et al. (2004) and DeDontney et al. (2012). The rupture propagates on the main fault passing through the junction point and, depending on the initial stress conditions, it may jump onto the branch fault. This is illustrated in Fig. 4, which shows a series of snapshots of slip velocity for the TPV24 and TPV25 test problems. The rupture propagation paths obtained in these two cases are remarkably different. In the right-lateral case the rupture abandons the main fault and continues into the fault branch. The rupture successfully continues in the branch fault because it is located in the extensional side, where the dominantly tensile normal stress changes tend to reduce the frictional strength. In the left-lateral case, the rupture mainly propagates into the main fault and continues only a short distance into the branch fault. In both scenarios, rupture on one fault segment past the junction casts a stress shadow on the other segment that inhibits its activation (e.g. Harris & Simpson 1998).
We compare our results to two independent methods, the MAFE code by Ma & Liu (2006) and FaultMod code by Barall (2009). Both are finite element codes with split nodes. We consider solutions on a grid with size 100 m (this is the average GLL node spacing). Fig. 5(a) shows the comparison of rupture times. The evolution of the rupture front from the three methods is in very good agreement on both fault segments and at all distances and directions from the hypocentre. A small discrepancy is observed at shallow depth on the branch fault in TPV24 when the rupture reaches the free surface. The time histories of slip rate produced by the three methods at selected fault locations, shown in Fig. 5(b), are also in qualitatively good agreement.
We found general agreement between the three methods, including the details of rupture initiation, propagation and arrest. These and other SCEC test results not shown here (but available through the SCEC repository) verify our implementation of dynamic rupture in the SPECFEM3D-SESAME code. The software is now suitable to solve complex problems of dynamic rupture with irregular fault geometry, while retaining the existing capabilities of the code for problems of wave propagation with complex media and irregular surface topography.
It is interesting to compare the accuracy of the SEM to that of lower order FEMs that, unlike FDM and BIM, can also treat rupture problems in complicated geometries, heterogeneous media and non-linear bulk rheologies. By design, high-order methods are more accurate than low-order methods for smooth problems. However, this property is not guaranteed a priori for dynamic rupture problems that are typically non-smooth. A quantitative comparison between results of different codes for the SCEC TPV5 benchmark problem conducted at different resolutions (also known as TPV205) is available from the SCEC/USGS Code Verification Web Server (http://scecdata.usc.edu/cvws/ ) based on metrics defined by Barall & Harris (2014). For this particular problem, it is found that SPECFEM3D achieves better agreement with the finest resolution simulation, the solution computed with the fourth-order FDM of Dalguer & Day (2007) with grid spacing of 12.5 m, than do four low-order finite element codes. At given mesh resolution (100 m or 50 m), the rms error of SPECFEM3D is at least twice smaller for rupture time (http://scecdata.usc.edu/cvws/metric_cvv1_u1/tpv5/ metric_cvv1_tpv5_ac_0.html) and about twice smaller for slip rate (http://scecdata.usc.edu/cvws/metric_cvv1_u1/tpv5/metric_cvv1_ tpv5_ar_0.html ).

Background and modeling scope
On 2011 March 11, a M w 9.0 earthquake stroke Japan and triggered a devastating tsunami, causing severe damage in cities and nuclear facilities along the Japanese coast. A combination of seismological, geodetic, bathymetric and tsunami observations revealed a remarkable depth dependency of the frequency content of the source. We define two frequency bands: an LF band from 0.01 to 0.125 Hz and an HF band from 0.5 to 1 Hz. Large slip (∼50 m) close to the trench was inferred by kinematic source inversions of seismic data in the LF band (e.g. Suzuki et al. 2011;Koketsu et al. 2011;Lee et al. 2011) and from differential multibeam bathymetry surveys (Fujiwara et al. 2011). Radiation in the HF band was identified in the deep areas of the plate interface by backprojection of teleseismic data (e.g. Meng et al. 2011). Downdip of the hypocentre the HF radiation was interspersed within a relatively slow rupture process.
Few dynamic rupture models have been proposed to investigate the physical features of this event. Duan (2012)   3-D dynamic models and concluded that the rupture around the hypocentre was enhanced by the stress accumulation due to the preceding M7-class earthquakes and triggered thermal pressurization of pore fluids in the near-trench area causing large slip, which promoted propagation of the rupture over a wide region. Ide et al. (2011) considered that an additional push to the earthquake rupture (slip reactivation) comes from the rupture front backpropagating from the free-surface after rupturing the trench, a phenomena usually observed in dynamic rupture simulations of dipping faults (Dalguer et al. 2001). Goto et al. (2012) used a 2-D inplane rupture model and showed that slip reactivation can result from heterogeneous stress distribution. Kozdon & Dunham (2013) proposed a 2-D model that accounts for depth-dependent material properties, curved fault geometry and seafloor geometry, and showed that despite velocity-strengthening properties at shallow depth, rupture can reach the trench. Huang et al. (2013Huang et al. ( , 2012 also used 2-D inplane dynamic rupture models to provide a physical interpretation of the depth-dependent frequency content of seismic radiation, the variations of rupture speed and slip distribution. They also find that the structure of the subduction wedge contributes significantly to the up-dip rupture propagation and the resulting large slip at shallow depth. Here we propose a minimalistic 3-D dynamic rupture model consistent with this depth-dependent frequency content of slip, where the shallow part radiates coherent energy at low frequency and the deep part at high frequency. The deep HF radiation is interpreted as the rupture of asperities in the bottom part of the seismogenic zone of the megathrust (e.g. Huang et al. 2012;Lay et al. 2012). We set the model parameters by trial and error, taking as a starting point the 2-D dynamic rupture models developed by Huang et al. (2012). The model presented here should be considered as preliminary; a refined model quantitatively constrained by geophysical observations will be presented elsewhere.
Our simulation also serves to illustrate the capability of the SEM to handle non-planar fault geometries and narrow subduction wedges. The model accounts for the free surface, the slope of the outer wedge and the curved geometry of the subduction interface and the trench (Fig. 6). The latter is based on a fault geometry adapted from Simons et al. (2011), which includes constraints from bathymetry, seismic reflection surveys and the Wadati-Benioff zone delineated by seismicity (e.g. Iwasaki et al. 2002;Miura et al. 2003). The software CUBIT generates high-quality hexahedral meshes even in wedges with small dipping angles close to the trench.

Model setup
We consider a homogeneous elastic medium with S-wave velocity 3470 m s −1 , P-wave velocity 5800 m s −1 and density 2700 kg m −3 . We assume the linear slip-weakening friction law and a distribution of asperities defined by heterogeneities of initial stress and critical slip distance D c (Fig. 7). We set an elliptical patch in the region of large slip and a collection of small circular asperities in deeper regions, mainly from 25 to 55 km depth. The number, size and separation distance of the small asperities are set by trial and error to achieve a moderate average rupture speed of 2 km s −1 downdip of the hypocentre.
On the main asperity the stress drop is set to 9 MPa and on the small asperities to 12 MPa. Null stress drop and a high strength excess (24 MPa) are prescribed in the background. In dynamic rupture models constrained by statistical observations, surfacerupturing earthquakes are characterized by a large area of negative stress drop that enhances energy absorption close to the free surface (e.g. Dalguer et al. 2008;Pitarka et al. 2009). Subduction zones with large accretionary wedges exhibit an upper stable sliding region due to the presence of unconsolidated gouge and clays (Marone & Scholz 1988). Therefore we imposed a negative stress drop (average −2.5 MPa) in the shallow part of the fault interface (Fig. 7b). Inspired by the hierarchical patch model developed by Ide & Aochi (2005) and Aochi & Ide (2009) where D c varies with the asperity size, we prescribe D c = 3 m on the large asperities, D c = 1 m on the small ones and D c = 6 m in the rest of the fault (Fig. 7c). The distributions of static (μ s ) and dynamic (μ d ) friction coefficients and normal stress over the fault are shown in Fig. 8. Rupture initiates by reducing the static friction coefficient in the nucleation area of radius 15 km (green circle in Fig. 8a), so that the initial static yielding stress (μ s σ n ) is slightly below the initial stress. This procedure does not alter the stress drop distribution shown in the Fig. 8(b), that is, no overstress has been applied on the nucleation patch.
We generated an unstructured mesh with element size of 2 km on the fault surface. The mesh accounts for the non-planar megathrust and seafloor geometries (Fig. 6) and contains elements close to the trench that are highly deformed. Fig. 12.8 of Cohen (2002) shows that the dispersion properties of spectral elements for wave propagation are not severely affected by the deformation of the elements; the main consequence is a significantly decreased time step. To verify the accuracy of our dynamic rupture simulations in a subduction wedge geometry of shallow angle, we constructed a reduced version of our Tohoku 3-D model, consisting of a narrow vertical cross-section passing through the hypocentre, and performed simulations with 2 km, 1 km and 500 m element sizes. Our results are described in the Appendix and show remarkable agreement of rupture times and peak slip velocity of the primary and secondary fronts at these three resolutions (differences are smaller than a few percent). This indicates that the features of our numerical solution at 2 km resolution discussed here are adequately resolved. In addition, the results in the Appendix illustrate how sensitive are the results to the geometrical assumptions in the immediate vicinity of the trench.

Results and analysis
An overview of the rupture history produced by our model is given in Figs 9-11, which show the spatial distributions of slip velocity at selected times, rupture time and rupture speed, respectively. In particular, Fig. 9 shows slip velocity in three frequency bands: LF 0-0.125 Hz, IF 0.125-0.5 Hz and HF 0.5-1 Hz. In the initial 40 s the rupture propagates mainly updip, starting slowly (about 1 km s −1 ) and gradually accelerating, while the downdip rupture front remains slow and weak. At t ≈ 40 s the updip rupture front reaches the shallow region of negative stress drop and its peak HF slip velocity decreases, while the downdip rupture starts breaking the deep asperities and generating intermittent HF radiation. At t = 55.2 s the rupture has reached the trench and has bounced back downdip. In the updip region, we obtain an average rupture speed of about 3 km s −1 (Fig. 10a), in agreement with finite source inversion models that resolve multiple fronts (e.g. Lee et al. 2011;Koketsu et al. 2011). The average updip speed varies widely among other kinematic source inversion models (Ide et al. 2011;Wei et al. 2012; Yagi & Fukahata 2011). In a narrow band close to the trench our dynamic model produces supershear rupture (Fig. 10b). As elaborated by Kaneko & Lapusta (2010), this is expected locally, at shallow depth in elastic media due to the phase conversion of SV to P-diffracted waves at the free surface. The efficiency of this supershear rupture mechanism can be reduced by low-velocity layers, velocity-strengthening friction (Kaneko et al. 2008) or inelastic or poro-plastic deformation (Ma 2012) at shallow depth. These additional ingredients are the subject of future developments of our model. As shown in Fig. 10(b), near t = 48.8 s a secondary downdip rupture front emerges at the trench, disconnected from the main updip front, and at t = 55.2 s both fronts have coalesced. At t = 65 s rupture of the deep asperities continues downdip of the hypocentre, with rupture speed of 3 km s −1 within the small asperities, 1 km s −1 in their surroundings and an average of about 2 km s −1 . By t = 81.2 s the rupture has started propagating southwards, at speeds of 2.8-3.3 km s −1 . At t = 107.2 s the rupture has broken the southern large asperity and has started propagating up to the trench.   Meng et al. (2011). Our simulation is in agreement with several first-order observations: the slow downdip rupture propagation and the more rapid migration of the rupture towards the South. Differences in rupture timing remain in some areas, which can be reduced by further tuning of the dynamic model parameters, in particular the strength, size and distribution of asperities (e.g. Huang et al. 2012). Fig. 12 summarizes the spatial distribution of LF and HF peak slip rate. The simulation reproduces the general pattern of the observations: LF and IF radiation occurs mainly in the shallow part of the plate interface, from 0 to 25 km depth, where slip is larger than 40 m, whereas HF radiation occurs essentially in the deep small asperities, below 30 km depth, and extends over 300 km along strike. This distinct behaviour is further illustrated in Fig. 12, which shows the slip-weakening curve, slip rate, slip and slip rate spectrum at two locations inside a shallow and a deep asperity, respectively. The deep asperity has a sharp slip rate peak of 10 m s −1 and relatively small slip. The shallow asperity has smoother slip rate with peak of 5.5 m s −1 but larger slip reaching 55 m. The slip rate spectra confirm the different frequency content of slip in these two asperities. We also find that the HF slip rate is localized near the rupture front whereas the LF slip rate lags behind it. Fig. 13 shows seafloor displacements from our simulation and from ocean bottom geodetic measurements at five locations (Sato et al. 2011). The agreement is fair on the vertical components close to the hypocentre and good in the horizontal components at all stations. Close to the trench the seafloor displacements reach 8 m vertically, consistent with the generation of a large tsunami, and 30-40 m horizontally. Even though our model is only conceptual, we obtained a peak value of moment rate comparable to that in the kinematic model of Lee et al. (2011), although it occurs earlier in our dynamic model (Fig. 14a). A second moment rate peak happens at about 120 s in both our dynamic model and Lee's kinematic model, although with different amplitude. Closely reproducing the moment release history is beyond the scope of our preliminary model. In addition, the macroscale representation of the source of the dynamic model is consistent with the kinematic model in the frequency domain as shown in the moment rate spectrum plotted in Fig. 14(b).

C O N C L U S I O N S
We added the capability to model spontaneous earthquake rupture dynamics in the unstructured version of the spectral element code SPECFEM3D. We compared our results of 3-D test problems of the SCEC/USGS dynamic rupture code validation project to other FEMs and found that rupture times and time histories of slip rate are in good agreement. To asses the efficiency of our implementation we performed a strong scaling analysis. The results demonstrate that the dynamic rupture implementation does not cause any significant load imbalance and does not affect the overall performance of the code. A weak scalability also gave satisfactory results. The unstructured SEM coupled with our dynamic rupture implementation makes use of a versatile mesh generation tool (CUBIT) that enables dynamic rupture simulations on complex fault systems, for instance a non-planar faults with branches. Meshing complex fault geometries is possible even with very deformed elements, for instance close to the trench in the subduction wedge of a shallowly dipping megathrust, although at the expense of requiring a small time step. We presented a dynamic rupture simulation of the 2011 M w 9.0 Tohoku earthquake, a complex megathrust event in a non-planar fault with small dip angle close to the trench, which illustrates the versatility and stability of the method. Our dynamic model includes fault heterogeneity and reproduces two important observed features of the Tohoku earthquake: high-frequency radiation in the deep areas of the plate interface and low-frequency radiation and large slip (∼50 m) at shallow depth close to the trench. Overall, our dynamic rupture implementation offers a great potential to simulate more realistic earthquakes in complex fault systems.

A C K N O W L E D G E M E N T S
This study was supported by the QUEST project (Quantitative Estimation of Earth's Seismic Sources and Structure) funded by the 7th Framework Programm of the European Commission, the ASCETE Project (Advanced Simulation of Coupled Earthquake and Tsunami Events) funded by the Volkswagen Foundation within the program 'New Conceptual Approaches to Modeling and Simulation of Complex Systems', by the US National Science Foundation (CAREER award EAR-1151926) and by the Southern California Earthquake Center (based on NSF Cooperative Agreement EAR-0529922 and USGS Cooperative Agreement 07HQAC0026). Simulations were done at the Swiss National Supercomputing Center (CSCS), under the production projects 'Development of Dynamic Rupture Models to Study the Physics of Earthquakes and Near-Source Ground Motion' and 'Development of a Database of Physics-Based Synthetic Earthquakes for Ground Motion Prediction'. We thank Yihe Huang for discussions on observations and modelling of the Tohoku earthquake and for sharing her results of a 2-D convergence test.

A P P E N D I X A : C O N V E RG E N C E T E S T F O R RU P T U R E S I M U L AT I O N I N A S H A L L O W-D I P P I N G S U B D U C T I O N W E D G E
To verify the accuracy of our simulations in a shallow-dipping subduction-wedge geometry with highly deformed elements near the trench, we construct a narrow (2 km wide) vertical cross-section passing trough the hypocentre of the Tohoku's 3-D model (see Fig. A1) and perform simulations with 2 km, 1 km and 500 m element size with spatial accuracy of fourth order (five interpolation points, NGLL = 5). The boundary conditions on the North and South faces of the domain are free stress. The value of the Kelvin-Voigt relaxation timescale η follows the mesh-dependent scaling described in Section 2.3. Fig. A2 shows the rupture time along the whole fault and the slip rate at a distance of 2 km from the trench for the three meshes.  The rupture times (Fig. A2-a) and peak slip velocity of primary and secondary fronts ( Fig. A2-b) show very good agreement, with differences smaller than a few percent between the 0.5 and 2 km resolution results. Moreover, we find convergent behaviour with increasing mesh refinement, that is, the differences between the 0.5 and 1 km resolution results are overall smaller than those between the 1 and 2 km resolution ones. In summary, our simulation at 2 km resolution adequately resolves the features of the model discussed in the main text.
We do not present here a more formal convergence analysis (rms difference as a function of element size, e.g. Day et al. 2005), but note that smooth convergence may not be expected in our unstructured mesh simulations unless even finer resolutions are considered, because of imperfect scaling of the artificial viscosity region upon mesh refinement. Our three unstructured meshes are not obtained by hierarchical refinement, that is, we did not split each element of the 2 km mesh into eight elements to obtain the 1 km mesh. Instead, they were generated independently, although based on exactly the same seafloor and megathrust geometry (Fig. A1, bottom panels). As a result, the shape of the Kelvin-Voigt damping zone, which by design encompasses all elements in contact with the fault surface, is not perfectly re-scaled through mesh refinement, especially near the trench.
In addition, we have found in a different set of simulations that differences in the geometry of the megathrust surface within 2 km of the trench can lead to significant differences in slip rate. For example, Fig. A3 shows the slip rate in a simulation with 500 m resolution in which a concavity (a slope break) in the megathrust was artificially introduced near the trench to avoid too small time steps. This minor change of the model geometry did not change the rupture time of the primary front but provoked quantitative differences of slip velocity larger than 10 per cent. These results illustrate the sensitivity of certain aspects of dynamic rupture to second-order geometrical features of the model in the immediate vicinity of the trench. Figure A3. Effect of model geometry near the trench on slip velocity. (a) Detail of two meshes with element size of 500 m. The mesh shown in the foreground (green) is the same as in the convergence study (Fig. A1). The mesh in the background (black) has a slope break in the megathrust close to the trench leading to a small concavity. (b) Slip velocity histories at a distance of 2 km from the trench (fault node 'S') obtained with the two meshes.