Discrete Element Modeling of a Subduction Zone with a Seafloor Irregularity and its Impact on the Seismic Cycle

Seafloor irregularities influence rupture behavior along the subducting slab and in the overriding plate, thus affecting earthquake cycles. Whether seafloor irregularities increase the likelihood of large earthquakes in a subduction zone remains contested, partially due to focus put either on fault development or on rupture pattern. Here, we simulate a subducting slab with a seafloor irregularity and the resulting deformation pattern of the overriding plate using the discrete element method. Our simulations illustrate the rupture along three major fault systems: megathrust, splay and backthrust faults. Our results show different rupture dimensions of earthquake events varying from tens to ca. 140 km. Our results suggest that the recurrence interval of megathrust events with rupture length of ca. 100 km is ca. 140 years, which is overall comparable to the paleoseismic records at the Mentawai area of the Sumatran zone. We further propose the coseismic slip amounts decrease and interseismic slip amounts increase from the surface downwards gradually. We conclude that the presence of seafloor irregularities significantly affects rupture events along the slab as well as fault patterns in the overriding plate.


Introduction
Seafloor irregularities are ubiquitous in subduction systems. They include seamounts, seamount chains, marine ridges, or oceanic plateaus with overall heights up to 6-7 km (e.g., Lallemand et al., 1992;Bilek, 2011, 2014), positioned between perpendicular to parallel to the trench (e.g., the Kyushu-Palau ridge in the Nankai zone and the Cocos ridge in the Costa Rica zone are at ca. 90°; the Nazca ridge in the Peru subduction zone and the Louisville ridge at the Tonga trench are at ca. 45°; the Zenisu ridge is about parallel to the Nankai trench and Loyalty ridge is also almost parallel to the New Hebrides trench). Since decapitating these irregularities during subduction is extremely difficult (Wang and Bilek, 2011), the role of seafloor irregularities on seismic activity in subduction zones remains mostly controversial.
For instance, Singh et al. (2011) conducted a deep seismic reflection survey in Mentawai offshore of the Sumatran subduction zone (Fig. 1a), and identified a seamount (irregularity) of 3-4 km height and 40 km width, which has been subducted to the depth of 30-40 km (Fig. 1b). According to the spatial pattern of seismicity, Singh et al. (2011) suggested that the presence of the seamount reduces the possibility of megathrust earthquakes in this zone. On the other hand, the paleoseismic records (Fig. 1c) of relative sea-level change from coral data in the same area imply that there could be megathrust events with magnitudes of up to 8.8. Moreover, the data covering 700 years of events shows that this area actually ruptures cyclically roughly every 200 years (Fig. 1c) (Sieh et al., 2008). The rupture pattern suggests that the next megathrust earthquake will probably take place within the next several decades, a potential devastating event for the coastal region of central Sumatra. Thus, one line of evidence suggests the likelihood of megathrust events, whereas another downplays that possibility.
The role played by subducting slab irregularities was partly addressed by analyzing fault patterns in the overriding plate, and particularly, fault generation and propagation (Dominguez et al., 1998a(Dominguez et al., , b, 2000, stress field disturbance (Ruh, 2016;Ruh et al., 2016), stress weakening as a function of the seamount dimensions (Ding and Lin, 2016), and fault pattern and dé collement caused by seamount-forearc collision (Morgan and Bangs, 2017). Although some approaches can model earthquake ruptures along the slab (e.g., Scholz and Small, 1997;Yang et al., 2012Yang et al., , 2013Yu et al., 2018), and earthquake cycling ruptures along the smooth slab (e.g., Van Dinther et al., 2013;Petrini et al., 2020), modeling of the associated fault patterns in the overriding plate as well as the slip rupture behavior along the slab with a seamount have not been achieved yet.
We study both fault generation and slip rupture behavior in a subduction zone using numerical experiments. More specifically, we apply the approach to the Sumatran subduction zone where a seafloor irregularity is present on the subducting slab (Singh et al., 2011) in order to analyze the earthquake recurrence interval and rupture dimension along the megathrust interface. We systematically compare our numerical results with seismic observations, seismic profiles, and paleoseismic records. We also discuss the deformation and fault structure of the overriding plate when subjected to subduction considering slip accumulation and rupture recurrence. Finally, we show that the presence of the seamount irregularity could be the source of a seismic cycle in the Sumatran subduction zone.

Methodology
To investigate the effect of a seafloor irregularity (e.g., seamount) on the behavior of a subduction zone and on rupture events, we performed numerical simulations based on the Discrete Element Method (DEM).
In this study, we simulate the deformation of the overriding plate as well as the interface behavior between the overriding plate and the subducting slab of the Sumatran subduction zone by adapting the model proposed by Jiao et al. (2018). The model was implemented in the Yade DEM software (see Donzé, 2008, 2009;Šmilauer, 2015 for details), whose flexibility allowed us to define a specific model relevant to study both faulting and slip along faults.

Model formulation
As in every DEM model, the simulated medium is discretized as an assembly of rigid particles interacting M8.4 earthquakes, respectively; (b) the interpreted cross-section of the CGGV040 shows the geometry of the subducting slab with a seamount (in a purple shape) at the depth of ca. 30-40 km; (c) earthquake cycles from the paleoseismic records at Bulasat (Sieh et al., 2008) show the big rupture events repeat every ca. 200 years and result the vertical slips of ca. 0.7-2.5 m. In addition, this paleoseismic model infers that the big rupture event at the Mentawai zone will repeat in the future. through specific interaction laws. The medium as a whole can deform and fracture based on the elastic-plastic behaviors described by these interaction laws. The overall behavior of the material is governed by the motion of its constitutive particles ruled by Newton's second law. The computing cycle can be decomposed into four main steps related, respectively, to (1) the determination of the positions of the constitutive elements, (2) the determination of their potential interaction, (3) the computation of the forces applied to each of them according to the predefined interaction laws, and (4) the calculation of their updated positions through the integration of the equations of motion. The calculation cycle is repeated iteratively until the simulation stops. Because of the dynamic formulation of the method (explicit time-domain integration), a non-viscous damping is used to dissipate kinetic energy and facilitate convergence towards quasi-static equilibrium. This damping directly acts on the interaction forces before their numerical integration in the equations of motion, so that the displacements are calculated from the damped forces. This is a convenient numerical tool to ensure the convergence of the simulations (see Duriez et al., 2016 for details).
The inter-particle behavior of our DEM model can be decomposed into the normal and tangential directions of the contact plane. The normal contact model accounts for both divergence and convergence (Fig. 2).
In the convergence regime (compression of the contact/ bond), the normal force F n is computed as: F n = K n · U n (1) where U n is the normal component of the relative displacement between particles A and B, and K n is the normal stiffness derived from the properties assigned to the particles, such that: where R A and R B are the radii of the particles and E A and E B , their respective elastic moduli, are directly related to the bulk modulus of the simulated medium.
In the divergence regime (extension of the contact/ bond), the normal force is computed with the same stiffness as that defined for the convergence regime. The inter-particle distance can increase up to U n tensile , for which the maximum admissible tensile force F n max is reached: with t the tensile strength of the interparticle bond and A int = π·{min(R A , R B )} 2 the interacting surface area between A and B. When F n max is reached, the force is not set to zero immediately as it is usually for modeling brittle rocks (e.g., Scholtè s and Donzé , 2013). Instead, F n gradually decreases, describing a softening behavior at the particle scale between U n tensile < U n < U n rupture , according to: where s is a weakening coefficient that needs to be defined. If the inter-particle distance continues to increase, the inter-particle bond breaks when U n > U n rupture and all forces are set to zero. A crack is then defined at the location of the bond breakage.
As in classic DEM formulations (Hart et al., 1988), the tangential force F s at the current time step t is computed incrementally as: with F n t−∆t the force computed at the previous time step, ∆U S the incremental tangential displacement between A and B, and K s the tangential stiffness, defined as K s = a· K n with a, a coefficient related to the Poisson's ratio of the simulated medium.
As for the normal force, a maximum admissible tangential force F n max is defined as: F s max = c· A int (6) with c the inter-particle cohesion (Fig. 2). Once the tangential force reaches this limit, the tangential force stays equal to its maximum value until the maximum admissible normal force is reached.
One crucial aspect of the model is that, unlike in classic DEM approaches where particles behave in a purely frictional way, once the interparticle bonds have broken, every new detected contact are set as cohesive bonds (Fig.  3). These bonds between particles have the same strength as the initial ones. This procedure provides a way to model the healing processes that eventually take place along faults. Since healed contacts are generally located in localized areas, where they tend to weaken the bonds (e.g., fault zones). This healing procedure also permits to control the dilatancy of the medium undergoing failure. (a) Normal behavior, (b) tangential behavior, and (c) failure envelope. Fig. 3. Illustration of kinematics and interparticle bond behavior for particles initially in contact (a), generation of cracks (b) and the new bond (c). The purple particle is moving towards the right side.

Model setup 2.2.1 Geometry and boundary conditions
To address the influence of a seamount on the rupture behavior along a subducting slab, we set up a 2D model corresponding to a cross section perpendicular to a trench. We defined the numerical model geometry based on the Sumatran subduction zone depicted by Singh et al. (2011). The numerical set up is shown in Fig. 4. The overriding plate is set up as a 50-km deep and 230-km wide wedge. The seamount (39 km wide, 6 km high) is located 30 km below the top of the overriding wedge and 160 km away from the trench, similarly to what could be observed for the subduction slab in Mentawai. The rigid slab is 10 km thick, a general thickness for oceanic slabs (Hayes et al., 2012).
All these components were filled in with spherical discrete elements. The mean radius of the elements is ca. 1 km and the porosity of the model is ca. 0.46. The seamount and the slab were modeled as rigid media moving together as one. The subduction rate was set to 6 cm/yr horizontally and to 1.2 cm/yr vertically (Singh et al., 2011). The right boundary of the overriding plate is fixed while the slab continuously subducts, which compares to the creeping patch related to the subducting slab.

Model scaling
We scaled the elastic and strength properties of the model in order to simulate the rupture/seismic cycle behavior in subduction zones.
In order to understand the stress-strain relationship of the material in the overriding plate, the model's parameters we implemented were shown in Table 1, then the numerical medium could show the elasto-plastic behavior (Supp. Fig. S1). Such behavior enables to simulate nucleation and propagation of faults associated to strain softening as well as sliding along faults in the residual state.
To scale the deformation properties, we weaken the elastic properties. We chose elastic stiffnesses one thousandth of the natural stiffness, so that the elastic strain limit of the numerical medium is ca. 10%. In nature, the elastic strain limit of rock (e.g., granite) is ca. 0.01% (e.g., Okubo and Fukui, 1996). Thus, considering the resolution limit of the model, we set that a 1-km deformation in the model is scaled down to a 1-m deformation in nature.
To scale the strength properties, we followed the classical scaling law (Hubbert, 1937(Hubbert, , 1951Peltzer, 1988). σ c represents the average strength of the lithosphere, h c is thickness, g c is the gravitational acceleration, ρ c is density, then, the dimensionless force ratio expresses the balance between the tectonic forces F t and the gravity induced forces F g . If the value of this ratio is similar in nature and in the numerical model, the model is scaled with respect to gravity.
For the continental lithosphere, we consider σ c = 5× 10 5 kPa, h c = 10 5 m, ρ c = 2.8 g/cm 3 , and g c = 9.8 m/s 2 , F t / F g =0.17. On the other hand, we use σ c ′ = 180 kPa, h c ′ = 50 km, ρ c ′ = 2.7 g/cm 3 , and g c ′ = 10 m/s 2 in our model, then we obtain F t ′/F g ′=0.13, similar with 0.17 obtained for the real continent. The strength of the model is from the uniaxial test (Supp. Fig. S1).

Model calibration
Before simulating the rupture behavior along the subduction zone, we calibrated our model by studying the stress-strain responses of the overriding plate and of the interface between the overriding plate and the slab, respectively (Fig. 5).
In order to understand the stress-strain relationship of the material in the overriding plate and of the interface between the overriding plate and slab, we extracted two sub-volumes from the model (Fig. 5a, b) to perform shear tests. The square block was used to characterize the plate behavior (Fig. 5a, c) while the purple block and white slab were used to characterize the megathrust fault behavior (Fig. 5b, d). The stress-strain responses and associated deformation fields obtained from both cases, which indicate the stress accumulation and release behavior along the faults (Fig. 5e, f). These tests confirm that this model is able to simulate the fault generation inside the overriding plate and the stress accumulation-and-release behavior along the generated or pre-existing faults. We use this model to simulate the rupture/seismic behavior (also stress accumulation-and-release behavior) of the megathrust fault (pre-existing fault) and the possible splay fault (generated fault) in the Sumatran subduction zone.

The case without a seamount
In order to address the effect of the seamount in the subduction zone, we present the slip velocities along the slab for the case without seamount. We simulate subduction in a period of ca. 400 years, the same as that in the case with a seamount in the following section. The slip along the interface between the overriding plate and the slab tends to be tremor or creeping events, their recurrence interval is about 10 years and the slip amount for each is about 6 cm (Supp. Fig. S2), since we did not consider other roughness on the slab. It is impossible to identify the different slip behavior between the front and back prisms for the case without the seamount. The effect of a seamount geometry on a slab on slip and fault behavior in the subduction zone will be simulated and discussed in the following model.

Preconditioning
Before simulating the rupture behavior of the subduction zone, the model needed to be preconditioned. The initial state of stress ( Fig. 6) is generated by first letting the prism stabilize under gravity and then, second, by moving the subducting slab until the occurrence of the first megathrust event which occurred after ca. 200 years of deformation in our model (Supp. Fig. S3). This loading sequence corresponds to the first stick-slip sequence that induces a quasi-elastic rebound of the overriding plate. We consider that the model is preconditioned right after the occurrence of this first rebound since further displacement of the subducting slab starts a new stick-slip sequence in the rupture cycle (seismic cycle) that leads to another rebound of the overriding plate. The rebound tends to reset the state of stress within the overriding plate. Thus, we record the further deformation from this stage, 200 years being the actual point of reference of time in our simulations.
This model enables us to investigate the earthquake cycle by estimating the frequency of large and small earthquakes as well as by characterizing the nucleation of splay faults caused by the subducting irregularity. The magnitudes of the events are scaled by their rupture lengths. In addition, we also analyzed the slip behavior occurring along the megathrust and splay faults and discuss their interactions.

Modeling Results and Discussion
To investigate the earthquake recurrence interval and rupture dimension along the megathrust interface, we simulated the subduction process up to the point where the overriding plate travelled ca. 24.5 m, corresponding to the (a) and (b) show the blocks we took from the entire model for the shear tests; (c) shows the simple shear test setting for the generated fault: the blue boundaries are fixed, and the shear loading is exerted on the yellow boundaries. Since we simulate the fault strength on average at depth ca. 15 km, we also apply ca. 400 kPa (hydrostatic pressure ρgh) lateral stress in this test; (d) shows the shear test for the pre-existing fault, which is part of the simulated megathrust fault: the blue boundary constrains the displacement of the purple block, and the loading velocity is ap plied on the white block; (e) and (f) show the stress-strain curves of the generated fault and pre-existing faults respectively. duration of approximately 400 years. The final deformation is shown in Fig. 7. In the following, we analyze how much deformation occurred above the slab and how the overriding plate responded to the constant downward slab motion.

Displacement field of the overriding plate
We calculated the displacement of the overriding plate during the simulated 400 years (Fig. 8). We also tracked the areas where interparticle ruptures occurred (called cracks hereafter) so as to identify fault zones. After the seamount moved with the slab by ca. 24.5 m, the cracks concentrated on three main regions: one along the interface between the slab and the overriding plate (the megathrust fault), and the other two on each side of the seamount; a front thrust fault (also called Splay Thrust, ST) that developed from the landward flank of the seamount towards the seaward side; and a Back Thrust (BT) that originated from the landward base of the seamount towards the landward side. The cracked areas corresponding to these two fault zones are relatively wider than the one along the megathrust fault. The BT is actually thickened with several short conjugate faults visible on Fig. 8a, b.
Beside the crack swarms, the displacement field ( Fig.  8a, b) is conspicuously separated by the ST and BT faults. The summit of the seamount divided the entire subducting interface into two patches. The shallow patch in the frontal accretionary prism is located at the seaward side of the model in front of the ST. This frontal accretionary prism suffered ca. 6 m of horizontal displacement with respect to the fixed boundary (Fig. 8a), and a maximal subsidence of ca. 2 m (Fig. 8b). The slab subducted for ca. 24.5 m during the 400 years. The movement of the upper plate was different from that of the slab because of the rupture events that occurred along the interface between the upper plate and the subducting slab. Due to the ruptures along the ST and the BT faults, the patch at the back of the prism suffered a ca. 3 m uplift (Fig. 8b) and a relatively small horizontal displacement (less than 1 m) compared to that of the shallow patch (Fig. 8a).
Due to the ST, which is related to the seamount, the deformation results show significant horizontal displacement at shallow depths (frontal accretionary prism, seaward side of the model), while maximal subsidence can be observed along the ST (color fields in Fig. 8a, b). At greater depths (back of the prism, landward side of model), horizontal slip is small, while vertical displacements are significant due to ruptures that occurred along the splay fault and the back thrust.
In addition, in order to measure the vertical displacement along the subducting slab, we recorded the vertical motion of a series of points P1-5 during the simulated 400 years (Fig. 8c). Three big events BE1-3 (vertical slip > 0.4 m during events at point P1 close to the trench) occurred in 272, 412, and 551 years, respectively, corresponding to a recurrence interval of 139.5 ± 2.5 (ca. 140) years. Even though the experiment only run for 400 years, the occurrence of these big events clearly seems periodic. The displacement paths of P2 and P3 are similar to that of P1, but these two points, which are close to the splay fault ST, show smaller uplifts compared to P1 during the big events BE1-3. P4 and P5, which are located on the back prism, are not sensitive to these big events.
The vertical displacement record of P1 (Fig. 8c) shows different tendencies before and after the big events. Here, we consider the deformation before and during the big Fig. 6. The initial state (initial balance situation) of the simulation.
The layers with different colors inside the overriding plate help identify how the overriding plate deformed; the key reference locations (points A to B just above the slab and points C to K ca. 5 km above the slab) are marked here for the later analysis.  event BE1 (i.e., in one cycle) as an example. We define the time between 200 to 272 years as the interseismic period before BE1, when the vertical deformation of most part of the frontal prism (e.g., P1) in this period was mostly oriented downwards. We define the time between 272 to 280 years as the seismic period of the BE1, when the vertical deformation of most part of the frontal prism in this period was mostly upwards (e.g., P1 in Fig. 8c). Even though the seismic period in nature is much shorter, about few minutes, we could only observe the fast rupture event occurred in these eight years from the model. Thus, we define this period as the seismic period. Before BE1, the frontal prism moved landwards (Fig. 9a) and downwards (Fig. 9b) with the slab. The back prism was raised by the splay fault and back thrust, delineated by the crack swarms. During BE1, the frontal prism rebounded seawards (Fig. 9c) and upwards (Fig. 9d). The entire overriding plate seemed to behave elastically, and the seaward motion increased gradually from the back to the frontal prism. The overriding plate generally moved upwards, except for the zone located just above the seamount, which moved downwards, affected by the seamount. During this period, there was few slip along either the splay fault or the back thrust fault. In details, the vertical displacement records (Fig. 9e) also show that P1-3, located in the frontal prism, generally moved downwards before BE1 and upwards during BE1. P5, located in the back prism, moved upwards before and during BE1. P4, located in the back prism but above the seamount, moved upwards before BE1 and downwards during BE1. The similar behavior can be identified in the other cycles (e.g., detailed deformation before and during the last event BE3 can be referred in Supp. Fig. S4).
During the coseismic period of the BE1 (i.e., the time between 272 to 280 years), the frontal prism uplifts significantly (vertical displacements at P1-3 in Figs. 9c, 9e). Such displacement could be associated with small ruptures along ST or/and cracks in the frontal prism with occurrences of BEs and immediately after BEs as their aftershocks.

Slip accumulation along the slab
In order to understand the slip cycle along the interface between the slab and overriding plate, we measured the slip deformation that occurred along the slab interface every 0.2 years (Fig. 10). The deformation corresponds here to the relative slip with respect to the slab and is thus different from the absolute displacement presented in Figs. 8, 9. Hence, the elements close to the fixed boundary, located at the bottom-right corner of the upper plate, show a constant slip rate (i.e., homogeneous increasing slip amount along the vertical axis in Fig. 10a) with respect to the subducting slab. Apart from those elements at bottomright corner of the upper plate, the ones located along the interface present various slip rate, implying stress concentration and release along the interface during the subduction process. We then identified the number of rupture events when the accumulated stress releases as well as their corresponding magnitude during the simulated period (400 years). Here, we only measured the slip along the interface as the slip along ST or BT can be interpreted through the displacement field inside the upper plate (Figs. 8, 9). The accumulated slip spectrum (Fig. 10) along the slab also shows that ST separates the frontal and back prism parts. The deformation pattern is dominated by the three big events (BE1-3 in Fig. 10; first labeled in Fig.  8) accompanied with smaller magnitude sub-events, labeled SE. The relative lighter areas (less density of gray Fig. 9. Displacements of the overriding plate during the first cycle (80 years), which is identified from Fig. 8c.  Fig. 8. Grey and red arrows show the movement direction of the overriding plates before and after the big event, respectively. lines) located in the frontal prism correspond to the large cumulated slip that developed during these big events. These ruptures propagated to the trench region, inducing deformation at the surface. Their rupture lengths along the slab interface were ca. 100 km. The SE only ruptured for ca. 50 km (Fig. 10). Some SE ruptures extended to the trench region (such as SE1-4), whereas some did not (such as SE5).
We analyzed the cumulated slip varying with the time and displayed the cumulated slip along the slab in Fig. 10. Such slip history can be further discussed with respect to key reference locations (A to K) located in the corresponding figures.
For further analysis of the deformation process, the segment from points B to G was considered in the main frontal prism. After BE1 and BE3, the accumulated slip amounts from point B to point G have similar features which shows the sub horizontal line from B to G at the top of the lighter areas in Fig. 10. It suggests that these big events may have completely released the accumulated stress along the slab and may have ruptured the interface between the slab and the entire frontal prism from A to G along the slab. The seaward side of point B (ca. 25 km from the trench landwards) including the area surrounding point A was affected by some shallow disturbances. After BE2, stresses were not yet totally released in some parts of the main frontal prism, so BE2 was followed by some smaller rupture events SE2-4, which appeared on the segment close to the trench. During each event, the stress in the frontal prism might be released partially. After these sub-events, the accumulated slips of the main frontal prism were homogeneous distributed over time (the top of the lighter area was not horizontal after the event; Fig. 10). After decades of stress accumulation, the stresses were partially released first during the sub-event SE5, and then entirely released by the next big event BE3.
The cumulated horizontal slip (Fig. 10b) is similar to the accumulated resultant slip, associated with the big events BE1-3 and the sub-events SE1-5, but the cumulated slip distribution history of the back prism (from point H landwards) is slightly different from the resultant one. The cumulated vertical slip field (Fig. 10c) shows more details about the deformation processes for the frontal prism and in particular for the back prism. In the frontal prism, the three big events and five sub-events readily resulted in vertical slips (Fig. 10c) as well as in cumulated resultant slip (Fig. 10a). The back prism suffered much more vertical slip than the frontal prism. The difference between the frontal and back prisms was in fact due to the slip along the splay fault. In addition, from the accumulated slip distribution in the back prism, we could identify smaller events, which ruptured the part from the seamount landwards.

Slip velocity along the slab
As mentioned in the last section, the density of slip amount distribution (gray or white areas) along the slab reflects the slip velocity variety. Different densities of the gray lines in Fig. 10 denote that the slip velocity was temporarily heterogeneous. Some sparser periods are associated to rapid rupture events, whereas some denser periods are associated to slow slip or to smaller events. However, the recurrence of rupture events cannot be clearly identified from the accumulated slip distributions. Hence, we recorded the slip velocities respect to the slab along the slab interface and plotted the resultant and vertical slip velocity spectrums in Fig. 11. Since the horizontal slip distribution (Fig. 10b) along the slab is similar to the resultant one (Fig. 10a), we only dissect the resultant and vertical slip velocity spectrums, shown in Fig. 11a and b, respectively.
In the slip velocity spectrum obtained from the whole simulation (Fig. 11), white regions (no slip velocity) represent totally locked patches, light yellow regions corresponds to partially creeping zones (slip velocity is less than 6.0 cm/yr in Fig. 11a or vertical slip velocity is less than 1.2 cm/yr in Fig. 11b), and yellow regions represent completely creeping zones (slip velocity is ca. 6.0 cm/yr in Fig. 11a or slip velocity is ca. 1.2 cm/yr in Fig. 11b). The dark regions indicate occurrences of rupture slip event (slip velocity greater than 6.0 cm/yr in Fig. 11a or vertical slip velocity greater than 1.2 cm/yr in Fig. 11b). In Fig. 11a, the events associated with large slip velocities, including BE1-3, show long dark lines in the spectrum and thus correspond to large rupture events (also see in Figs. 8c, 10). Besides the three big events (BE1-3), the subevents (SE1-5), shown as the relatively dark color and shorter lines than the big events, happened at years 345, 433, 464, 475, and 505, respectively, corresponding to a recurrence interval of 40 ± 29 (ca. 50) years, which correspond to about one third of the big events recurrence interval (ca. 140 years). The likelihood of these recurrence intervals probably depends on the deviation of slip accumulation or on the distribution of the stress concentration along the slab as discussed in the previous section. Besides these large events, which occurred along the interface between the slab and the frontal prism, the vertical component of the slip velocity spectrum (Fig. 11b) reveals small events along the interface between the slab and the back prism, probably also along the splay fault, e.g., events in seismic group 1 (SG1, taking place in 240-270 years) and seismic group 2 (SG2, taking place in 330-380 years), which occurred at the back prism segment after the first purple dashed line (the seaward base of the seamount).
In Fig. 11b, some smaller events (e.g., events in SG1 and SG2) in the back prism with rupture length less than 50 km, which is smaller than that in the frontal prism. However, from the slip velocity spectrum on the top surface of the overriding plate shown in Fig. 12a, we could also identify these events with longer length, probably because of the dip of the splay fault. Additionally, a low slip velocity region, located between the seaward base (the left purple dashed line in Fig. 11b) and the summit (the gray dashed line) of the seamount, suggests that this segment tended to lock with the seaward flank of the seamount. Such behavior could have favored the formation of the splay fault. This low slip velocity region is less evident in the surface slip velocity spectrum presented in Fig. 12b.
The slip velocity spectrum shows characteristics of every type of earthquake (big and sub events) in a subduction system. The big events (i.e., BE1-3) ruptured the entire frontal prism along the slab to the trench, with rupture length of ca. 100 km and a regular recurrence interval of ca. 140 years. The sub events, associated with slips along different parts of the interface between the frontal prism and slab, ruptured along lengths of ca. 50 km with a recurrence interval of decades. Some of these sub events (e.g., SE1-4) ruptured the segment close to the trench. Some (e.g., SE5) produced slip along moderate length segments of the slab and may have accumulated stress in the shallow part of the prism. Some small ruptures, which produced slips along the interface between the back prism and the slab (i.e., events in SG1-2), reached the surface of the splay fault to generate ca. 50 km length   Surface slip amount also show the big events and sub-events; but the light yellow region (between the left purple and gray dashed lines) is not as clear as that in Fig. 11; the splay fault group rupture longer over the entire seamount, because of the dip of the splay fault (Figs. 8 and 9). The left purple, gray and right purple dashed lines are the same as those in Fig. 11. ruptures. Most of ruptures stopped at the seamount.

Slip sequence of the reference points above the slab
In order to elucidate the slip heterogeneity along the slab, we computed the accumulated slip amounts along the slab at different reference points above the slab (Figs. 13, 14).
The resultant slip curves (Fig. 13a) are used to measure the displacement of the given positions along the slab during events in the frontal prism. The cumulated slip increased gradually from point A to point F. The slip slope generally increased from A to F between rupture events (e.g., before BE1; between BE1 and SE1), implying that the inter-rupture slip rate increased from the trench landwards.
The resultant slips (Fig. 13a) evidence large coseismic slips (ca. 4 m) of the frontal prism during big events, BE1-3, corresponding to the coseismic deformation. Since point A is located at the trench region, it could be affected by shallow disturbance, resulting in different slip behaviors from that of the other points of the frontal prism. The slip curves of points B-F clearly indicate the three big events. Nevertheless, the coseismic slips from C to F seemed to decrease, because the seamount (close to point F) stops the coseismic rupture. Furthermore, we can also discern slip associated to the sub event SE1 between BE1 and BE2, and slips associated to the four sub events SE2-5 between BE2 and BE3, whereas these sub events are difficult to be identified in the vertical slip curves in Fig. 13b.
The slip patterns of all the points located along the slab in the frontal prism are generally similar to each other (points A-F in Fig. 13), but different from those located in the back prism (points G-K in Fig. 14; even though location of point G belongs to the frontal prism, we discuss this point together with points in the back prism). The slip amount accumulated during 400 years changes from ca. 20 m at F to ca. 7.5 m at G, because G is close to the seaward base of the seamount and it is thus strongly affected by the seamount subduction. In Fig. 14, it is difficult to identify the big events contrary to Fig. 13, since those ruptures did not propagate over the seamount. From G to K, the general cumulated slip amount increases (Fig. 14). As mentioned previously, the seaward flank of the seamount (points G and H) seemed to lock with the seamount as evidenced from the vertical component curves (Fig. 14b). Points I-K show similar resultant slip patterns, which can be associated to some inter-rupture slow or small slips (slip amount less than 1 m) that occurred during smaller events (indicated by gray arrows in Fig. 14a), which probably belong to the seismic groups (e.g., SG1 or SG2 in Figs. 11, 12). The vertical slip patterns at points I-K are different due to change of dipping directions when the points went over the summit of the seamount.
From the trench landwards in the frontal prism, the accumulated slip (Fig. 13) somewhat increased during inter-seismic periods whereas it decreased during coseismic periods. Between the rupture events, the accumulated slip increased gradually from the trench landwards, showing that the shallower the segment was, the more locked the interface is during the inter-rupture periods. Thus, between rupture events BEs and SEs, the accumulated slip in the frontal prism at distances larger   (BE1-3), indicate the small corupture slips during smaller events from points I to K curves. Points G and H are close to the seamount, thus co-rupture slip behaviors around them are more complex. The slips generally increase from points G to K. Their slip patterns show more co-rupture slips of frequent small events, which are related to the event slips along the splay fault. The vertical component pattern shows the local seamount geometry effect. The key reference locations (points G to K) are denoted in Fig. 6. than 50 km away from the trench (from point C landwards) also increased gradually from the trench landwards; at shallower level of the interface obtains higher coupling ratio (i.e., seismically locked instead of creeping). The only exception is the region very close to the trench (from Point A to C, ca. 10 km depth), where was unstable due to its proximity to the surface. The region from points C to F (ca. 50-210 km from the trench and ca. 10-45 km depth), gathered the big rupture events (BE1-3). During these rupture events, the coseismic slip could reach ca. 3 m in this region along the slab. In the back prism (Fig. 14) and on the landward of the seamount, the region (ca. 45-50 km depth) can be regarded as a relatively small seismogenic region, since most of events rupture for less than 1 m along either the subducting interface or the splay fault. Since this region is close to the creeping boundary (backstop boundary), some interrupture slips prevailed.

Comparison of Modeling Results to Observations
We discuss here our simulation results with respect to observations related to the Sumatran subduction zone. First, we analyze the effect of the seamount geometry on the deformation of the overriding plate and seismic cycles, and we compare the vertical displacement on the surface with the paleoseismic record in the subduction zone. Then, we compare nucleated splay faulting in the overriding plate with recently seismic observations. Finally, we compare the simulated varies rupture zones with the observed seismic domains along the slab.

Analysis of the effect of the seamount at the subduction zones
The displacements inside the overring plate (Figs. 8,9) and along the slab (Figs. 13, 14), the accumulated slip distribution (Fig. 10), and the slip velocity spectrum (Figs. 11,12) reveal that the seamount divides the deformation of the overriding plate in two distinct zones separated by the splay fault, ST. The seamount acts as a watershed, thus making the displacement paths different in the frontal and back prisms at the surface (Fig. 8c) as well as along the slab (Figs. 13,14). The frontal prism tends to subduct before the big events and rebound during the big evens, whereas the back prism does not respond to big events. Because of the seamount, the big events and sub events (Figs. 13,14), which occurred in the frontal prism, do not affect the back prism zone. The coseismic slip behaviors around the zone close to the seamount are more complex; in particular, the region in contact with the seaward flank of the seamount seems to move along with the slab during the subduction (Figs. 11b, 12b, and 14).
The rupture length and the seismic cycle (Figs. 10,11,and 12) are also different in the two regions separated by the seamount. The rupture length (ca. 100 km) in the frontal prism is longer and the cycle (ca. 140 years) is also longer than that in the back prism (ca. 50 km with few decades interval). The slip patterns (Figs. 11,12) show dominant long slips (ca. 100 km) associated to a few big events in the frontal prism, and more smaller slips (less than 50 km) associated to frequent small events, which also include the small event slips in the back prim, probably along the splay fault.
In addition, the geometry of the seamount significantly influences the vertical slip along the interface between the slab and the overriding plate. The seamount also accounts for the generation and propagation of the splay fault and back fault in the overriding plate. Even though splay faults might be easily generated at the shallow zone with the weak sediments and friction along the slab, it is not easy to be formed by the only friction along the slab at the depth of ca. 30 km. However, the splay fault propagation is observed at Mentawai zone (Wang et al., 2018). Thus, this observed geometry of the seamount might be the key factor in the different motion pattern of the frontal and back prism, as well as the cycle of the big events, and it may also be the key factor in generating and rupturing the splay fault at the depth. In the model, ruptures of the big events also terminate at the position of the seamount, since the seamount is the mean feature to accumulate the stress along the slab, and further controls the big rupture cycle. Thus, the prominent irregularity such as the seamount is the source of the big events at the subduction zones and also the splay fault propagation at depth. Besides, the rupture along the splay fault could interact with seismic events on the subduction slab. In the coseismic period of the BE1 (i.e., the time between 272 to 280 years), the frontal prism continuously uplifts (positive vertical displacements at P1-3 in Figs. 8c, 9e). Such displacement might be associated with ruptures along the splay fault. If the splay fault ruptures together with the interface between the overriding plate and the slab, a larger magnitude is expected that could result in severe ground shaking and maybe tsunami. If the splay fault rupture follows a megathrust event as an aftershock, its occurrence could be simulated through our DEM model or a physics-based model, e.g., Coulomb stress change and rate-and-state friction model (Dieterich, 1994), providing crucial constraint for seismic hazard assessments (Chan et al., 2017).

Comparison of surface displacement and geological observations
As described above, the vertical deformation of most part of the frontal prism (e.g., P1) goes downwards (Fig. 8c) in between big events (e.g., BE1) and generally rebounds up during big events. The surface deformation is thus significantly affected by the seismic cycle, as documented in Sieh et al. (2008). Since the projection of the location of Bulasat (Fig. 8) on the seismic profile is close to P1 (red lines in Figs. 8c, 15), we compared the vertical deformation path of P1 and the paleoseismic record of Bulasat (Figs. 1c,15) (Sieh et al., 2008). Even though the amplitudes of the vertical displacements of P1 and of Bulasat are somewhat different, the overall evolution before and during the events are certainly same. Both in the model and observations, this area (around P1 or the Bulasat) subsides and uplifts before and during big events respectively.
The vertical displacement record of P1 (Figs. 8c, 15), accumulated slip distribution (Fig. 10) and slip velocity spectrum (Fig. 12) shows the recurrence interval of every earthquake (big and sub events) in the subduction system. The vertical slip close to P1 on the surface is ca. 0.5 m, comparable to the paleoseismic observations (0.7-2.5 m) from Sieh et al. (2008). The big events (i.e., BE 1-3) rupture the entire frontal prism along the slab to the trench, corresponding to rupture lengths of ca. 100 km. They also occur regularly every ca. 140 years, which is overall comparable to the ca. 200 years interpreted from paleoseismic data (Sieh et al., 2008). Several reasons can cause the discrepancy between the model and the observation. In the model, the geometry of the interface between the slab and overriding plate is simplified. We focused on the main feature (irregularity/seamount) along the slab, and we thus simplified the geometry of the slab as a planar slab. Furthermore, the properties of the material are also not certain, and we did not take into account any heterogeneity, which most probably exist in such tectonically active zones. The real cycle behavior is controlled by many factors (e.g., the chemistry condition inside the overriding plate or along the slab interface), and so far, there are still some factors, which cannot be well constrained. However, through testing our model, we found that the strain-softening is the key factor in the fault generation and propagation. With strain-softening behavior and scaled stress-strain relationship of the model, our model is able to simulate the seismic cycle and preexisting fault (megathrust) or newborn fault (splay and back fault) propagation.

Comparison of splay fault and seismic observations
During the seamount subduction, ruptures occur along the megathrust fault, splay thrust fault and back thrust fault. We also compared our simulation results with the seismic activity observations in the Mentawai region of the Sumatran subduction zone (Wang et al., 2018). Fig.  16a, b present snapshots showing typical crack swarms occurring during one year at the beginning and at the end of the simulated period, respectively. The locations of the ST lineation denoted by the crack swarms at these two periods (presented in Fig. 16a, b, respectively) are similar. The crack swarms simulated along the ST is comparable to the seismic activity reported by Wang et al. (2018). Using a seismic catalog recorded by a global network, Wang et al. (2018) could identify two earthquake lineations: one along the slab with a low-angle thrust faulting and another along the splay fault with a highangle mechanism. Both seismic groups can be clearly shown in our model. The position of the crack swarms along the ST fault and the slab interface are consistent with the seismic observations. Moreover, one of the frontal backthrust (FBT), main backthrust (MBT) and coastal backthrust (CBT) on the surface could correspond to the point where the splay fault propagates to the surface. The width of the faulted backthrust zone is also similar to what can be observed in Figs. 8, 9. Thus, the depth, position and dip angle of the splay fault formed in the model is in good agreement with the recent seismic observations.

Comparison of the different rupture zones to observations
The shallow, big and small seismogenic regions defined in our model can be compared to the near-trench, central megathrust and downdip domains identified from the seismic observations by Lay et al. (2012).
Along the interface between the frontal prism and slab, the degree of the creeping slip increases from the near trench to the deep regions, whereas the slip during the big events decreases. Here, we can divide the rupture along the slab into three distinguished regions: shallow seismogenic, big seismogenic and small seismogenic regions. The shallow seismogenic region (0-10 km depth) corresponds to the shallowest part of the plate. This region was characterized as the near-trench domain (depth shallower than 15 km) from interpretations based on seismicity records. The downward region (10-45 km depth) concentrates stress and produces large slips during the big rupture events. This region corresponds to the large seismogenic region, which is comparable to the central megathrust domain identified at 15-35 km depth by Lay et al. (2012). In the deeper region at 45-50 km depth, the rupture might be small when it occurs along the slab only, but it can also be great when the rupture develops along the slab and propagates along the splay fault. This region could thus correspond to the modest seismic region described by Lay et al. (2012). If a rupture in this region does not propagate along the splay fault and arrests at the seamount, it could be associated with small events which would result in slow slips and seismic tremors, corresponding to the downdip or transitional domain identified by Lay et al. (2012). Thus, the region division along the slab with seamount subduction rigorously matches the domain definition from seismic observations (Lay et al., 2012).

Conclusions
In this study, we utilized DEM simulations to study the subduction zone at the Mentawai segment in the Sumatra subduction system and inferred the impact of irregularity on the seismic activity in a subduction system. We have simulated the deformation and rupture behavior of the overriding plate during 400 years of the irregularity/ seamount subduction. Our main findings are: (1) The modeling results shows different rupture lengths of events, which varies from tens km to 140 km.
(2) The crack swarms are located in three main regions. One is aligned along the interface between the slab and the overriding plate. The other two regions are located in the splay and backthrust fault regions. The position of the simulated splay fault and its angle to the slab are in good agreement with the seismic observations.
(3) The proposed model is able to simulate fault nucleation and propagation, and permits to evaluate the recurrence of the megathrust ruptures. Along the interface between the slab and the overriding plate, the seaward flank of the seamount, which tends to be fully locked with the slab, separates the rupture behavior along the slab into two main segments. The shallow part ruptures as a result of the big events, with a regular recurrence interval of ca. 140 years during the modeled 400 years. Each big event ruptures ca. 100 km on the interface. Besides these big events, there are five sub events, and the repetition of these events is about decades with significant irregularity. The rupture length of these sub events is ca. 50 km. Some subevents located close to the trench can propagate toward the surface, whereas some located more at depth cannot.
(4) On the landward side of the seamount, the deep part shows more small rupture events. If these events propagate along the splay fault toward the surface, the length of these ruptures might be greater than 50 km. During these events, the rupture along the splay fault tends to uplift the region between the splay and backstop faults.
(5) This is an innovative study using numerical simulation to reproduce and demonstrate the four kinds of seismogenic regions along the slab: shallow, big, modest, and small seismogenic regions, comparable to the domains defined from seismic observations (Lay et al., 2012).
Though this study, we propose that the subduction zones with the prominent irregularity/seamount located along the slab, such as the Mentawai zone, still present the hazard of large earthquakes whose occurrence follows the seismic cycle. In particular, if we have observed the rupture along the splay fault at depth and the seamount geometry along the slab, it could happen a megathrust earthquake in the near future (less than ca. 100 years; one earthquake cycle), which will generate strong ground shaking or/and tsunami.
At the end, we would like to claim that since no model is perfect, the more data constrain, the more reliable the model would be. At the Sumatran subduction zone, the current available data could only provide us limited constrains to our model. Furthermore, our current model did not consider the 3D geometry in the subduction zone, as the models proposed by Furuichi et al. (2018). Due to the lack of information, our model did not consider the temperature, hydraulic pressure effects in the different subduction depth. We expect more observations in the future, and then we could further improve the model (e.g., a 3D model with some more physical factors, such as hydraulic stress and temperature distribution, given by petrini et al., 2020).
grateful to the Asian School of the Environment (ASE) for allowing him to complete this work by extending his contract at NTU for one year and a half. This work is Earth Observatory of Singapore contribution (No. M4430217.B50.706022) and the Ministry of Science and Technology (Grant Nos. MOST 109-2116-M-008-029-MY3, MOST 110-2124-M-002-008, and MOST 110-2634 -F-008-008). We note that there are no data sharing issues in the paper since all of the numerical information has been provided in the table and figures. Jiao Liqing carried out the simulations and prepared the manuscript. Chan Chung-Han contributed to the result interpretation, wrote manuscript, and administrate project. Luc Scholtè s assisted in the development of the methodology and revised the manuscript. Auré lia Hubert-Ferrari contributed ideas and edited the manuscript. Fré dé ric-Victor Donzé contributed on the development of the methodology. Paul Tapponnier leaded the funding acquisition.

Code/Data Availability
We implemented the Yade DEM software (see Kozicki and Donzé, 2008;2009;Šmilauer, 2015 for details), whose flexibility allowed us to define a specific model relevant to study both faulting and slip along faults.