1. Introduction

Figure 1(a) shows the tectonic setting of Japan. The Bungo Channel is located between the islands of Shikoku and Kyushu in southwest Japan (Fig. 1(b)). The Philippine Sea plate is thrusting against the continental plate at a rate of 7 cm/yr in the northwest direction near the Bungo Channel (e.g., Miyazaki and Heki, 2001). In 1997 and 2003, southeastward transients occurred as a result of interplate slow slips, which had similar moment magnitudes (Mw 7.1), similar slip areas, and approximately 1-year durations (Hirose et al., 1999; Ozawa et al., 2004). Considering the time interval of approximately 6 years between the 1997 and 2003 events, it has been pointed out that the next Bungo Channel aseismic slip may occur soon.

Fig. 1
figure 1

(a) Tectonic setting of Japan. Solid lines indicate plate boundaries. (b) Enlarged map of the rectangular area in (a). The Bungo Channel separates the islands of Shikoku and Kyushu.

Full size image

Hereafter, we use the phrase ‘long-term slow slip’ to mean aseismic slip that lasts for one year or more. By contrast, we use the phrase ‘short-term slow slip’ to mean aseismic slip that has a moment magnitude of roughly 6 and lasts for only several days as detected using a tiltmeter (Obara et al., 2004).

The continuous GPS network of the Geospatial Information Authority of Japan (GSI) has detected transient crustal deformation in the Bungo Channel area since early 2009. The spatial pattern of the detrended transient motion shows southeastward movement, which suggests that aseismic interplate slip occurs in this area, taking into account the northwestward subduction of the Philippine Sea plate.

In this research, we estimate the spatiotemporal evolution of the aseismic slip along the plate boundary beneath the Bungo Channel area since January 2009 and compare it with previous long-term slow slips.

2. Data and Analytical Procedure

GPS data for 24 h were analyzed on a daily basis with Bernese GPS software (version 5.0). We adopted the F3 solution (Nakagawa et al., 2008), which uses the final orbit and earth rotation parameters of the International GNSS Service (IGS) and provides a higher S/N ratio than the previous F2 position time series (Hatanaka et al., 2003). We transformed the daily F3 coordinates into local displacement time series (east-west, north-south, and up-down components) with respect to the Misumi site (Fig. 1(b)) for 94 GPS sites around the Bungo Channel.

Since the raw time series include annual and trend components, we first estimated the annual components. We used a polynomial function and trigonometric functions to fit the position time series for the period between 1 January, 1996, and 1 January, 2011. The degree of the polynomial function and the overtone of the trigonometric functions were determinedby Akaike’s information criterion (Akaike, 1974). Since the polynomial function whose degree was optimized by this criterion fitted the data sufficiently well to follow long-term slow slips, we consider that the estimated annual components are unaffected by long-term slow slips. After estimating the annual components, we removed them from the raw time series together with the linear trend estimated for the period between 1 January, 2007, and 1 January, 2009. In the detrending process, we assumed that the crustal deformation rate for the period between 1 January, 2007, and 1 January, 2009, was constant. Hereafter, we refer to this procedure as the detrending of the time series; that is, the detrended time series is the deviation from the positional trend for the period between 1 January, 2007, and 1 January, 2009, with the annual components removed.

Using the detrended time series (Figs. 2, 3 and 4), we applied square-root information filtering (Ozawa et al., 2007) by means of the time-dependent inversion technique (McGuire and Segall, 2003) for the period between 1 January, 2009, and 1 January, 2011. Since we used detrended data, the estimated aseismic interplate slip is the deviation from the steady state for the period between 1 January, 2007, and 1 January, 2009. The GPS sites used in the filtering analysis are plotted in Figs. 3(a)–(e). In our inversion analysis, we weighted east-west (EW), north-south (NS), and up-down (UD) movements with a ratio of 1:1:1/3, considering the standard deviations. We used the fault patch shown in Fig. 3(a), which consists of 10 × 10 knots. The interval between nodes is approximately 24 km in the dipping direction (west-northwestward) and 15 km in the strike direction (south-southwestward). The plate boundary model was borrowed from Satake (1993) and Kyushu University (1999). We incorporated the inequality constraint that the aseismic slip opposes the motion of the subducting Philippine Sea plate, following the method of McGuire and Segall (2003). Aseismic slip distributions for the 2003 and 1997 Bungo slow slips were estimated by Ozawa et al. (2004, 2007) by a similar method. We reanalyzed the 2003 and 1997 events in order to compare the 2010 event with the 2003 and 1997 long-term slow slips under the same fault patch and constraints and using the most precise F3 solutions, which have been available since 2009.

Fig. 2
figure 2

Detrended position time series at selected GPS sites. (a)–(d) correspond to the 2010 event, while (e)–(h) and (i)–(l) correspond to the 2003 and 1997 events, respectively. Southeastward transients are clearly shown. The locations of the GPS sites are indicated in Fig. 3(a). Red lines represent values computed using our model.

Full size image
Fig. 3
figure 3

Detrended movements (black arrows) and computed values (white arrows). Ellipses centered on the tips of black arrows indicate 1σ errors in observations. The tails of white arrows show the locations of GPS sites used in our analysis. Southeastward movements indicate interplate aseismic slip. The rectangle shows the fault patch used in estimating interplate aseismic slip. We used a variable scale of vectors to show characteristic features in the time evolution of each event. (a)–(e) correspond to the 2010 event, (f)–(i) correspond to the 2003 event, and (j)–(m) correspond to the 1997 event. (a) 2009/1/1–5/1. (b) 2009/5/1–9/1. (c) 2009/9/1–2010/1/1. (d) 2010/1/1–10/1. (e) 2010/10/1–2011/1/1. (f) 2003/5/1–7/1. (g) 2003/7/1–11/1. (h) 2003/11/1–2004/2/1. (i) 2004/2/1–5/1. (j) 1996/12/10–1997/6/1. (k) 1997/6/1–11/1. (l) 1997/11/1–1998/2/1. (m) 1998/2/1–5/1.

Full size image
Fig. 4
figure 4

(a) Detrended EW and NS position time series at site 0437, whose location is shown in Fig. 3(a). (b) Detrended EW and NS position time series at site 0681 (Fig. 3(a)). In panels (a) and (b), solid lines represent values computed using our model.

Full size image

3. Results

We can clearly see transient movements in Figs. 3(a)–(e), which show southeastward displacements. The transient motion for the period between 1 May and 1 September, 2009, is up to 1 cm/year as shown in Fig. 3(b), which is estimated from a coordinate change between May and September. The total displacement during this period is larger than the error level and consistently directed to the southeast direction (Fig. 3(b)). The standard deviation of the displacements was estimated by Kalman filtering on the assumption that it does not change over time. A value of 1σ typically corresponds to displacements of 0.2 mm for horizontal components and 0.6 mm for vertical components. These values are translated into 1 mm/year for horizontal components and 3 mm/year for vertical components in Fig. 3(b). The spatial pattern in Fig. 3(b) strongly indicates the occurrence of a slow slip event on the Philippine Sea plate. Figures 4(a) and 4(b) show position time series at two GPS sites (0437 and 0681), whose locations are shown in Fig. 3(a). The time series of ground displacements in Figs. 4(a) and 4(b) show a southeastward change between May and September 2009, especially in the NS component. The ground displacements level off from September to October and the southeastward transient continues with increasing speed over time from approximately November 2009 (Figs. 4(a) and 4(b)). We cannot clearly state that this transient between May and September was caused by a short-term slow slip, since its duration appears to be roughly 2 months from June to July (Fig. 4).

Figure 3(c) shows transients for the period between 1 September, 2009, and 1 January, 2010. The spatial pattern in Fig. 3(c) shows southeastward movements in southwest Shikoku and small transients in Kyushu. After January 2010, the transient motion increased in speed and GPS sites in Kyushu started to record southeastward displacements (Figs. 3(d) and 3(e)). Figure 2 shows detrended position time series at selected GPS sites associated with the 1997, 2003, and 2010 events, indicating similarity among the three transients.

On the basis of the above detrended time series, our filtering analysis shows that interplate aseismic slip occurred in the period between May and September 2009 beneath the Bungo Channel and southwest Shikoku, covering a low-frequency earthquake area (Fig. 5(b)) (Katsumata and Kamaya, 2003). The interplate aseismic slip for this period is much larger than 1σ error as shown in Fig. 5(b). The estimated moment magnitude from 1 May to 1 September 2009 is 6.3, with a rigidity of 30 GPa.

Fig. 5
figure 5

Estimated interplate aseismic slip. Green arrows show movements of the continental plate against the Philippine Sea plate. Ellipses centeredon the tips of green arrows indicate 1σ errors. Red circles represent the epicenters of low-frequency earthquakes during the corresponding period. Solid lines indicate isodepth contours with 20-km intervals. Areas enclosed by solid lines show the adopted fault patch. (a)–(e) correspond to the 2010 event, (f)–(i) correspond to the 2003 event, and (j)–(m) correspond to the 1997 event. These periods correspond to those in Fig. 3. The slip area first appeared beneath southwest Shikoku and expanded into the Bungo Channel, over time, during the 2003 and 1997 events. We used a variable scale of vectors to show characteristic features in the time evolution of each event.

Full size image

Our model for the period between May and September 2009 does not reproduce particularly well the observed transients for this period, shown in Figs. 3(b) and 4. This poor fitting is due to the oversmoothed slip distribution both in space and time in our analysis. Because we used a single hyperparameter to control the temporal evolution of the fault slip, oversmoothing occurs in the case of a small event. Thus, the estimated aseismic slip distribution cannot explain the large displacements in the northern part of Cape Ashizuri. However, the crustal deformation from May to September clearly shows a typical pattern induced by a slow slip event in this area. Thus, we consider that a slow slip event occurred on the upper surface of the Philippine Sea plate during this period, although we cannot decide whether it is a part of the following large slow slip, or an isolated event. This hypothesis is supported by the fact that the observed displacements and estimated interplate aseismic slip are much larger than 1σ error together with high activity of low-frequency tremors for this period (Figs. 5(a) and 5(b)).

Furthermore, Between 1 September, 2009, and 1 January, 2010, the aseismic slip area shifted to beneath the southwestern part of Shikoku, adjacent to the area ruptured by the 1946 Nankai earthquake (e.g., Sagiya and Thatcher, 1999), while the slow slip in the west and central parts of the Bungo Channel gradually vanished (Fig. 5(c)). We consider that the backslip in Fig. 5(c) is not significant. Aseismic slip continued beneath southwest Shikoku between January and October 2010 and expanded into the central Bungo Channel area (Figs. 5(c)–(e)). Figure 5 shows the 1σ errors of the estimated slip parameters, indicating that the estimated slip is well above the 1σ errors. The root-mean-square of the residual between the results of observations and computations using the 2010 model is 2 mm, while those obtained using the 2003 and 1997 models are 3 and 4 mm, respectively. The estimated moment magnitude from 1 September, 2009, to 1 January, 2011, is 7.0 (Fig. 6). The time evolution of the moment magnitude is similar to that of the number of low-frequency earthquakes, indicating a close relationship between long-term slow slip events and low-frequency earthquakes or short-term slow slips (Fig. 6).

Fig. 6
figure 6

(a) Time evolution of the magnitude of moment released by the 2010 event. (b) Cumulative number of low-frequency earthquakes in the area surrounded by blue lines in Fig. 5(a).

Full size image

The resolution of the model parameters is shown in Fig. 7 along with the total slip distribution. The resolution matrix R is expressed as

, where H is the coefficient matrix of the observation equation tions, and HT is the transpose of the matrix H (Menke, 1989). The resolution indicates the reliability of the estimated slip. When the estimated resolution is near 1, we can conclude that a slip in this area has high reliability. Our resolution in Figs. 7(d)–(f) shows the high reliability of the estimated parameters within the fault patch.

Fig. 7
figure 7

(a)–(c) Estimated total interplate aseismic slip. Green arrows indicate the slip of the continental plate against the subducting Philippine Sea plate. Red circles represent the epicenters of low-frequency earthquakes for the corresponding period. The nomenclature matches that of Fig. 5. (d)–(f) Contour maps of resolution of slip parameters. The contour interval is 0.1 m.

Full size image

4. Discussion

From the above results, we draw two main conclusions: a small aseismic slip first occurred beneath the Bungo Channel and southwest Shikoku between May and September 2009 and decayed from September to October 2009 (Figs. 3–5). This was followed by a long-term slow slip beneath southwest Shikoku, which expanded into the Bungo Channel from January to October 2010 with increasing slip speed and then decayed over time (Figs. 5(d) and 5(e)). Our approximate model reproduces these observations well, as shown in Figs. 2 and 3, although it does not closely fit the observations in Fig. 4, which is attributed to spatial and temporal oversmoothing for small events.

From 1 May to 1 September 2009, the activity of low-frequency tremors became higher than that for a period between January and May 2009. However, there were no reports of short-term slow slip in the Bungo Channel or southwest Shikoku for this period (National Research Institute for Earth Science and Disaster prevention (NIED), 2010), although our analysis indicates the occurrence of a slow slip there. The reason for this may be that the duration of the event was sufficiently long for tiltmeters, which are unstable over long periods, to miss the signal.

The Bungo channel slow slip events can be divided into three periods: a rise period, a high-speed period, and a deceleration period. The slip gradually increases in magnitude during the rise period, attains its maximum speed during the high-speed period, and gradually decreases its speed during the deceleration period. These phases appear in the time evolution of the estimated moment shown in Fig. 6.

The long-term slow slip beneath southwest Shikoku increased in speed from approximately March 2010, the time at which the activity of low-frequency earthquakes became higher (Figs. 5 and 6). The reason for this may be that the long-term slip in this event originated in a shallow region away from the low-frequency earthquake area and, over time, expanded into a deeper region near the low-frequency earthquake area, making the stress state more favorable for low-frequency earthquakes to occur during the period of high-speed slip. At the time of the 2003 long-term Bungo Channel slow slip, low-frequency earthquakes also became more active when the slip area expanded into a deeper region as the slip accelerated (Ozawa et al., 2007). The rate of stress change appears to have a strong effect on the occurrence of low-frequency earthquakes, which become more active during the high-speed period of a long-term slow slip.

At the time of the 2003 event, a small slip occurred beneath southwest Shikoku during a rise period (Fig. 5(f)). Since the inverted slip area extends to one of the edges of the modeled fault area, we checked whether the characteristic features change if we use a larger fault patch, and we found they did not change significantly. Then, during the high-speed period (Fig. 5(g)), the center of the slip area moved westward by approximately 30 km (relative to its position in Fig. 5(f)). At the time of the 1997 event, the slip area appeared to be centered on southwest Shikoku (Fig. 5(j)) during a rise period. Then, during the high-speed period and deceleration period (Figs. 5(k)–(m)), the center of the slip area moved westward by approximately 15 km (relative to its position in Fig. 5(j)) while the slip magnitude in the central Bungo Channel area increased over time. Thus, the rupture process of the long-term slow slip after September 2009 was similar to those of the previous Bungo Channel slow slips, in that the rupture originated in a shallower part beneath southwest Shikoku and extended into a deeper region beneath the Bungo Channel during the highspeed period. Our 2003 and 1997 models reproduce these observations well, as shown in Figs. 2(e)–(h) and 2(i)–(l), respectively.

The slip distribution of the 2010 slow slip was similar to those of the 2003 and 1997 events, as shown in Fig. 7. The total moment magnitude of the 2010 event since 1 September, 2009, is 7.0 (Fig. 6), which is similar to that of each previous event (Mw 7.1).

Thus, the 2010 event was similar to the previous long-term Bungo Channel slow slips in terms of both extent and magnitude. This suggests a characteristic slow-slip hypothesis that long-term slow slips in the Bungo Channel occur not only with a quasi-periodicity of 6 to 7 years but also with a similar moment magnitude, slip area, and duration.