Numerical simulation of afterslips and slow slip events that occurred in the same area in Hyuga-nada of southwest Japan

SUMMARY We propose a model of numerical simulation for the coexistence of afterslip for ∼ M 7 earthquake and slow slip events in the Hyuga-nada region of Japan that incorporates 3-D geometry of the Philippine Sea Plate. Coseismic slip events, recurrence of slow slip events and afterslip are qualitatively reproduced using the composite law, which is a type of rate- and state-dependent friction law with higher cut-off velocity. In addition, characteristic slip distances in the area are larger than those in other seismic source areas. In our simulation, afterslip, which occurred at the velocity-weakening regime, triggered an aseismic transient event. After the termination of this event, spontaneous slow slip events repeatedly occurred in the same area. After a similar event sequence was repeated, another afterslip occurrence triggered a larger coseismic slip in a wide area including that of the slow slip events. Following this coseismic slip, the aseismic slip area was locked until the next afterslip propagation, which triggered an aseismic transient event. These results suggest that detailed observation of spatial and temporal distribution within the area of aseismic slip may indicate the potential of recurring slow slip events and future large earthquakes.


I N T RO D U C T I O N
Several large thrust earthquakes of M = 6.7-7.2 have occurred with recurrence intervals of a few decades in the Hyuga-nada region at the western edge of the Nankai trough, southwest Japan. Although the details of their source parameters are not well understood, Yamashita et al. (2011) suggested that the locations of source regions in 1899 (M7.1), 1931 (M7.1), 1961 (M7.0), 1996 October (M6.9) and 1996 December (M6.7) were nearly the same (Fig. 1). They consider that the second 1996 earthquake occurred in an unruptured area of the first earthquake of the same year. Therefore, the recurrence interval of earthquakes in this area is approximately 30 yr. Additional large earthquake sets occurring in 1941 and 1970 in close proximity to the 1996 sequence ( Fig. 1) indicate the existence of two coseismic slip areas. Between them, we focus on the southwestern area because aseismic slips were observed at its down-dip extension. Moreover, afterslips of the two 1996 earthquakes (Yagi et al. 2001) and slow slip events (SSEs) have also been observed in this region (Yarai & Ozawa 2010).
The aseismic slip area is located in Hyuga-nada, off the east coast of Kyushu Island and south of Bungo Channel, where SSEs occur with recurrence intervals of several years (e.g. Hirose et al. 1999). Although SSEs have previously been observed only at the Bungo Channel near Kyushu Island, three such events have been detected since 2005 in the aseismic slip area with a recurrence interval of approximately 2 yr and durations of 0.5-1 yr (Yarai & Ozawa 2010). Their study was the first to suggest that SSEs occur within the afterslip area. However, the relationship between afterslip and SSE remains unknown. In addition, the studied afterslip/SSE area is located within the tsunami source area of a 1662 thrust earthquake.
To date, the largest earthquake (M7.6) in the Hyuga-nada region occurred in 1662 (Headquarters for Earthquake Research Promotion 2004) with an estimated tsunami source area that expanded 100 km in a north-south direction (Hatori 1985).
In this study, we focus on aseismic slips during interseismic periods of large earthquakes to show the relationship between afterslip and SSE. We numerically simulate both afterslips and SSEs in cycles of large earthquakes by using a model that incorporates realistic 3-D plate geometry.
Numerical results based on rate-and state-dependent friction laws (e.g. Dieterich 1979) show that afterslips can occur under frictional conditions of velocity strengthening (e.g. Marone et al. 1991;Perfettini & Ampuero 2008;Hetland & Simons 2010) and velocity weakening (Helmstetter & Shaw 2009;Hori et al. 2009;Hori & Miyazaki 2011). Although SSEs also can occur under the conditions of velocity weakening and large nucleation size (Kato 2003(Kato , 2004Yoshida & Kato 2003), velocity strengthening with a C 2012 The Authors  (Yagi et al. 1999) and afterslip more than 0.04 m (Yagi et al. 2001), respectively, of the 1996 October earthquake. The blue represents those of the 1996 December earthquake (Yagi et al. 1999(Yagi et al. , 2001. Colored and grey ellipses of long dashed-dotted lines correspond to the main shock and aftershock areas of the 1931and 1961earthquakes, and 1941and 1970earthquakes, respectively (Yamashita et al. 2011. The red star represents the epicentre of the 1899 earthquake. The purple dashed ellipse represents the tsunami source area of 1662 earthquake (Hatori 1985). Grey contours indicate depth to the upper surface of the descending plate. Inset at upper right shows the location of the Hyuga-nada region and study area, indicated by a grey filled ellipse and red rectangle, respectively. one-state variable in the rate-and state-dependent law cannot cause SSE. Among the many numerical models proposed for generating SSEs (e.g. Liu & Rice 2007Shibazaki & Shimamoto 2007;Ariyoshi et al. 2009;Matsuzawa et al. 2010;Segall et al. 2010), our model was developed on the basis of a simple friction law with a small set of frictional parameters. We adopted the composite law by Kato & Tullis (2001) under the condition of velocity weakening to model afterslip and repeating SSEs, which occur at the same place during the cycle of large earthquakes and we discuss a detailed relationship of the spatial-temporal distribution between these aseismic phenomena.

Numerical model and simulation
Interplate earthquakes, afterslips and SSEs are modelled to represent the release of slip deficit accumulated during interseismic periods owing to a slip velocity slower than that of plate convergence (e.g. Rice 1993). The time derivative of shear stress τ due to plate boundary slip, is calculated as where V is the slip rate in the plate convergence direction; V pl , the plate convergence rate; G, rigidity (30 GPa) and β, the shear wave speed (3.27 km s -1 ). In this study, the subducting plate is discretized into small subfaults. The index j indicates the jth small subfault. K ij is the slip response function representing stress on a small subfault i due to unit slip on a small subfault j. The second term in eq. (1) is a radiation-damping term approximating energy radiation as an elastic wave when the slip rate is high (Rice 1993). Parameters determining frictional properties on the plate interface are assumed to represent a laboratory-derived rate-and statedependent friction law (Dieterich 1979;Kato & Tullis 2001). Eq.
(2) represents a fault constitutive law (Nakatani 2001) that determines the slip rate V i for given stress τ i and the value τ si (= τ s * i + τ si ). If we set the reference velocity V * at a level of slip velocity under which slip is sufficiently slow to be negligible, the value τ si is the threshold level of stress τ i to cause significant slip velocity. Hence, the value τ si is analogous to the 'strength as a threshold' (Nakatani 2001) and we call it as strength. Parameter A = aσ controls the slip increase rate at which point stress reaches strength. We do not treat effective normal stress σ directly because its effect is difficult to separate from that of frictional parameter a (in addition to b, as described subsequently) for natural earthquakes. V * is set to V pl in the calculation below. τ s * = μ * σ is the steady state strength with V = V * and τ s is the variation in strength from the steady state.
When we define strength as such, eq. (3) is an evolution law for strength change τ s , which varies depending on the prior slip history. Parameters B = bσ and L control strength recovery and slip weakening. For slip weakening, B and L mainly determine the amplitude of variation in strength τ s and slip weakening distance D c . This equation is a composite law that combines slowness (aging) and slip laws (Kato & Tullis 2001). The first term of eq. (3) represents time-dependent healing and the second represents slip weakening. V c balances these two terms. When V i V c , the first term of eq. (3), the time-dependent healing process is negligible.
As suggested by Ampuero & Rubin (2008), the range of parameter values permitting SSEs for the slip law would be much narrower than that for the slowness law. Further, the qualitative difference in slip pattern between the two laws originates from the difference in healing behaviour (Ampuero & Rubin 2008). These indicate that the time-dependent healing term in the composite law, which is equivalent to that in the slowness law, can play an important role in SSE reproduction. However, the healing term has been neglected in the velocity range of SSEs because of the low cut-off velocity V c = 10 −8 m s -1 in the composite law (Kato & Tullis 2001;Kato 2003Kato , 2004. To adopt the healing process even in the velocity range of SSEs (approximately 1 m yr -1 ), we set V c at 10 −2 m s -1 and the healing term is neglected only in the coseismic velocity range. The result of the different V c value is briefly mentioned in Section 4.
To solve eqs (1)-(3), we removed the time derivative of stress dτ i /dt and obtained differential equations for slip rate V and strength τ s . Assuming the initial values given below, differential equations are solved with an adaptive time-step fifth-order Runge-Kutta algorithm (Press et al. 1996). For the initial condition, slip velocity V was assumed to be uniform at 0.9 V pl and strength was calculated For the boundary condition, slip velocity was assumed to be constant and equal to V pl outside the model region.

Geometry of the plate interface
We used discretized 3-D Philippine Sea Plate interface geometry (Baba et al. 2002) to calculate the slip response function. We first divided the plate interface into small rectangular subfaults, each of which was approximated by three triangles for calculation of angular dislocation (Cominou & Dundurs 1975). Slip response functions are represented by the stress change on a subfault, which is the combined effect of three angular dislocations within a subfault. This stress change is estimated at the centre of a subfault at the direction of slip acceleration. We set small subfaults with lengths of 0.45 km in the x direction and 0.05 km in the vertical direction (x axis is shown in Fig. 2). In this study, the number of subfaults was reduced to shorten calculation time by using large sizes of apparent subfaults at deep and shallow regions in which stable slip occurs and no unstable patch is set. We calculated stress changes by the combined effect of 81 (9 × 9) small subfaults at depths of 6-15 km and 51-60 km and 9 (3 × 3) small subfaults at depths of 15-18 km and 45-51 km.
We set the plate convergence rate as constant (6.5 cm yr -1 ) on the basis of the slip deficit distribution reported by Heki & Miyazaki (2001). The convergence direction of the subducting plate was N55 • W.

Frictional property distribution
In this study, we assumed depth-dependent heterogeneity of A -B. As shown in Fig. 2(a), we set A -B < 0 (velocity weakening in steady state) in the seismogenic zone for a depth of 6-45 km. In addition, we set constant A -B for a depth of 10-40 km and A -B = 0 for a depth of 45 km.
The characteristics of simulated slip on the circular fault patch are controlled by the ratio of the patch radius r to the critical fault radius r c , which is expressed by eq. (4) (Kato 2003(Kato , 2004: When r r c , ordinary earthquakes with high slip rates repeatedly occurred at the fault patch and postseismic aseismic slip followed on the fault outside the patch. When r ∼ r c , episodic aseismic slip events occurred. Stable sliding occurred when r r c . We set L = 12 m on the background (Fig. 2b) and A -B = −0.35 MPa at depths of 10-40 km (Fig. 2a). Then, the value of r c was more than 900 km at that depth. The length of the model region for the strike direction (∼2r) was less than 200 km. Then, r r c , which satisfies a condition of stable sliding (Kato 2003(Kato , 2004. In addition, this value of L is as large as that of the background area described in Hori & Miyazaki (2011). All of the assumed L values in this study are much larger than laboratory values of ∼10 −5 m (e.g. Marone & Kilgore 1993) but are comparable to several recent experiments of ∼1 m (e.g. Chambon et al. 2002).
We applied patch 1 at the location of the 1996 October large earthquake with a radius of 10 km on the basis of the report by Yagi et al. (2001). In this patch, we set a significantly smaller L (8.4 cm) and a larger |A -B| of 0.67 MPa than that of the background. The value of r/r c was 2.9, indicating that unstable slip is expected to occur on patch 1 (Kato 2003(Kato , 2004.
For patch 2, we set the same value of A -B as background and a larger L (26.4 cm) than that at patch 1 at the location of observed aseismic slips with a radius of 20 km. The value of r/r c was 0.97, which indicates conditional instability. Thus, episodic aseismic slip events are expected to occur at patch 2 when this patch is separated from other patches (Kato 2003(Kato , 2004.
Hereafter, among expected aseismic slip events, we regard that spontaneously occurring aseismic events are included in the definition of 'SSEs'. Further, 'afterslip' describes an aseismic slip with a continuously decreasing slip rate following a large earthquake and aseismic slip that is not attributed to an SSE or afterslip is referred to as an 'aseismic transient' event.
At the beginning of the numerical simulation, we set small subfaults with lengths of 0.9 km in the x direction and 0.1 km in the vertical direction at all depth ranges. Further, we conducted approximately 40 simulation runs by varying the spatial distribution of frictional parameters. Next, we set smaller subfaults with lengths of 0.45 km in the x direction and 0.05 km in the vertical direction, C 2012 The Authors, GJI, 190, 1213-1220 Geophysical Journal International C 2012 RAS Downloaded from https://academic.oup.com/gji/article-abstract/190/2/1213/645282 by guest on 30 July 2018 which are described in Section 2.2, to set a smaller value of r c for patch 1 and we conducted approximately 30 simulation runs. In the following section, we describe the result for the parameter set that most quantitatively reproduced observation data among these simulation runs.

Earthquake cycles and recurrence of SSEs
The resultant temporal variations in cumulative slip are shown in Fig. 3. In simulation, the largest earthquakes (EQ0 and EQ4 in Fig. 3; M w = 7.2) occurred at patch 2. During this interseismic period, four SSEs (SSE1-SSE4 in Fig. 3) occurred spontaneously at patch 2 and three large earthquakes (M w = 6.8) ruptured patch 1 (EQ1-EQ3 in Fig. 3). The recurrence interval of the large earthquake at patch 1 was 67-77 yr and that of the largest earthquake at patch 2 was approximately 218 yr. A snapshot of slip velocity for 210 yr from the occurrence of EQ1 is shown in Fig. 4. Moment magnitudes of simulated earthquakes are calculated from the coseismic slip distribution, where coseismic slip is defined as slip with a slip rate equal to or greater than 1.0 cm s -1 with the rigidity of 30 GPa.
In addition, Fig. 4 shows afterslips (afterslip1-3 in Fig. 3) of these earthquakes that distribute around the seismic source and propagate to deeper regions. These afterslips triggered aseismic events (aseismic transient 1 and 2) and a seismic event (EQ4) at patch 2. Spontaneous SSEs (SSE1-SSE4) occurred twice after the termination of an aseismic transient event when the slip rate was slower than the plate convergence rate.

Triggering of aseismic or seismic events by afterslips
Afterslip following EQ1 slowly propagated to the deeper region. Slip at the shallower point within patch 2 (pink line in Fig. 3b) was accelerated within 0.5-1 yr after EQ1. At this time, slip at the central point of patch 2 (red line in Fig. 3b) did not accelerate because the area was locked (Fig. 4b). Slip rate on the periphery of patch 1 (afterslip 1) close to the plate convergence rate and the aseismic slip (aseismic transient 1) at the centre of patch 2 (Fig. 4c) accelerated approximately 4.5 yr after EQ1. The slip velocity of the triggered aseismic transient 1 was much slower than that of afterslip 1 (Fig. 3b).
Temporal variation in slip velocity following EQ1 (Fig. 3b) showed two peaks corresponding to slip at the shallow (pink line) and central (red line) points of patch 2 despite the relatively short distance of ∼15 km between these two points. Slip area of afterslip 1 distributed around the seismic source for 3.75 yr after EQ1 (Fig. 5a) and slip for aseismic transient 1 during the next 3.72 yr was distributed to patch 2 (Fig. 5b).
After the occurrence of EQ2, a locked area was not observed at patch 2 (Fig. 4i). Slip at the central point of patch 2 accelerated for approximately 300 d after EQ2 (Fig. 3c), triggering aseismic transient 2 (Fig. 4k). The slip area of afterslip 2 distributed around the seismic source for 86 d after EQ2 (Fig. 5c) and slip for aseismic transient 2 was distributed at patch 2 during the following year (Fig. 5d). Slip was distributed mainly at patch 2 during afterslip 2 and aseismic transient 2 (Fig. 5e).
Although afterslip following EQ3 (Fig. 4r), EQ1 (Fig. 4b) and EQ2 (Fig. 4j) propagated to deeper regions, a larger earthquake (EQ4), rather than an aseismic transient event, ruptured patch 2, 279 d after EQ3 (Figs 3d, 4s and 4t). The locked area at patch 2 after the occurrence of EQ4 (Fig. 4u) shrank slowly, delaying the acceleration of aseismic slip such as aseismic transient 1. The triggering of either an aseismic transient event (aseismic transient 2) or an earthquake (EQ4) at patch 2 by afterslip following patch 1 is determined by slight differences in stress and strength around the source; however, a detailed discussion is beyond the scope of this paper.

D I S C U S S I O N
As expected from the value of r/r c , we qualitatively reproduced the recurrence of spontaneous SSEs at patch 2, although the recurrence interval of more than 10 yr was longer than the observed value. At patch 1, we reproduced earthquakes that correspond to M7class earthquakes with an observed recurrence interval of a few decades in the Hyuga-nada region. In addition, it was determined that larger seismic events such as EQ0 and EQ4 also occur with longer recurrence intervals. In our study, EQ0 and EQ4 roughly correspond to the 1662 largest recorded earthquake even though the source area of EQ0 and EQ4 is smaller than that of the actual 1662 earthquake. The former may become comparable to the latter if we add an unstable patch representing the 1941 and 1970 earthquake areas in our simulation. However, such alteration was not applied for simplicity.
Regarding spatial distribution, afterslips following EQ1-EQ3 in our results propagated to deeper regions, which are roughly consistent with observations (Yagi et al. 2001;Miyazaki et al. 2003). The slip of afterslip 1 in our simulation (Fig. 5a) concentrated in narrower area than that of the observations and slip area of afterslip 2 (Fig. 5c) did not propagate to patch 2. These results suggest that the propagation time of afterslip to patch 2 in our simulations was longer than that observed. The slip area of afterslip during the 100 d after the second 1996 earthquake reported by Yagi et al. (2001) was similar to that of the superposition of afterslip 2 and aseismic transient 2 (Fig. 5e). Thus, the slip area reported by Yagi et al. (2001) may include an aseismic transient event triggered by afterslip. Miyazaki et al. (2003) showed that aseismic slip occurs after the termination of afterslip, more than 100 d after the second 1996 earthquake. Then, temporal variation of the moment rate for the aseismic slip shows two peaks, as shown in Fig. 3b. These suggest that the aseismic slip is triggered or occurs spontaneously. The former corresponds to aseismic transient 1 and 2 in this study and the latter corresponds to SSE 1 and 3.
We obtained similar results regarding the recurrence of SSEs and earthquakes even when we used a 5 per cent larger or smaller value of L and r for patch 2. However, only small SSEs occurred at patch 2 when we used a 27 per cent larger L for patch 2. On the contrary, when we set an 18 per cent smaller value of L for patch 2, earthquakes ruptured patch 2 continuously without SSEs; however, SSEs occasionally occurred during interseismic periods at patch 2. In both cases of V c = 10 −8 m s -1 and V c = 10 −2 m s -1 , we detected triggered aseismic transient event and spontaneous SSEs. However, the recurrence frequency of spontaneous SSEs and triggered aseismic transient event during interseismic periods at patch 2 decreased to 2 and 1, respectively. In addition, the recurrence interval of earthquakes at patch 2 was 30 per cent shorter than that in the case of V c = 10 −2 m s -1 . A similar result was obtained when we set a 20 per cent larger value of r at patch 2.
Under these conditions, trade-offs exist between parameter values and parameter distribution. Thus, determination of the best-suited parameter sets is difficult. To reduce the trade-offs and obtain a robust model, we should analyse GPS observations with resolution as effectively as possible. Because our model already explains several observations qualitatively, we can constrain the location and size of patches through comparison with fine observation data. We will then choose suitable frictional parameter sets that show greater consistency with the observation without trade-offs between patch geometry and frictional properties.
Several previous studies suggested that spatial heterogeneity in frictional properties on a plate interface is mainly related to the   location and frictional properties of large earthquakes, through comparison of observed data and model results for northeastern Japan (e.g. Kato 2008;Hori & Miyazaki 2010) and the San Andreas Fault (e.g. Johnson et al. 2006). In this study, however, we show the importance of heterogeneity within the area of aseismic slip. Therefore, by detailed observation of spatial and temporal distribution within the area of aseismic slip, we can estimate frictional properties for future large earthquakes and hidden frictional heterogeneity that cause aseismic and seismic events.

CONCLUSION
We calculated earthquake cycles in the Hyuga-nada region by using a model that incorporates realistic plate geometry. We were able to reproduce afterslips, recurrence of spontaneous SSEs and large earthquakes at conditionally unstable patches under the frictional conditions of velocity weakening by considering the time-dependent healing process even in the velocity range of SSEs. In addition, we reproduced aseismic transient events triggered by afterslip following nearby large earthquakes. These results suggest that it is possible to describe frictional heterogeneity on a plate interface by spatial and temporal distribution within the area of aseismic slip.

A C K N O W L E D G M E N T S
We thank two reviewers, Naoyuki Kato and Eric Hetland and editor Jeannot Trampert for constructive comments that improved this manuscript. We also thank Shinichi Miyazaki and Keisuke Ariyoshi for valuable discussions. This work was partly supported by the project "Evaluation and disaster prevention research for the coming Tokai, Tonankai and Nankai earthquakes" of the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. The Earth Simulator was used for all simulations. We thank Dr. William Stuart for providing the subroutines for calculating stress fields due to triangular dislocations, which he coded with Dr. Robert Simpson. We used Generic Mapping Tools (Wessel & Smith 1998) to draw all figures.