Estimation of slip rate and fault displacement during shallow earthquake rupture in the Nankai subduction zone

Abstract

Enormous earthquakes repeatedly occur in subduction zones, and the slips along megathrusts, in particular those propagating to the toe of the forearc wedge, generate ruinous tsunamis. Quantitative evaluation of slip parameters (i.e., slip velocity, rise time and slip distance) of past slip events at shallow, tsunamigenic part of the fault is critical to characterize such earthquakes. Here, we attempt to quantify these parameters of slips that may have occurred along the shallow megasplay fault and the plate boundary décollement in the Nankai Trough, off southwest Japan. We apply a kinetic modeling to vitrinite reflectance profiles on the two fault rock samples obtained from Integrated Ocean Drilling Program (IODP). This approach constitutes two calculation procedures: heat generation and numerical profile fitting of vitrinite reflectance data. For the purpose of obtaining optimal slip parameters, residue calculation is implemented to estimate fitting accuracy. As the result, the measured distribution of vitrinite reflectance is reasonably fitted with heat generation rate $\left(\dot{Q}\right)$ and slip duration (t _r ) of 16,600 J/s/m² and 6,250 s, respectively, for the megasplay and 23,200 J/s/m² and 2,350 s, respectively, for the frontal décollement, implying slow and long-term slips. The estimated slip parameters are then compared with previous reports. The maximum temperature, Tmax, for the Nankai megasplay fault is consistent with the temperature constraint suggested by a previous work. Slow slip velocity, long-term rise time, and large displacement are recognized in these fault zones (both of the megasplay, the frontal décollement). These parameters are longer and slower than typical coseismic slip, but are rather consistent with rapid afterslip.

Background

Coseismic slip behavior of megathrusts shallower than seismogenic zone is a key to understanding slip propagation and the generation of tsunamis. The shallow parts of the thrusts are generally considered to be aseismic with stable slip that inhibits earthquake nucleation and coseismic rupture propagation (Scholz 1998). However, during the 2011 M_w 9.0 Tohoku-Oki earthquake, Japan, coseismic rupture on the megathrust reached the seafloor, producing a huge slip displacement (approximately 80 m) (Ito et al. 2011; Ide et al. 2011). In addition, recent high-velocity (1.3 m/s) friction experiments on fault zone material in subduction zones showed that an earthquake rupture from deeper portions propagates easily through the updip boundary by high-velocity weakening associated with fluid pressurization and frictional heat (Ujiie and Tsutsumi 2010; Ujiie et al. 2013). On the other hand, onland seismic network observations have recorded different types of seismic events called very low frequency (VLF) earthquakes within accretionary prisms (Obara and Ito 2005; Ito and Obara 2006), and recent observations using ocean bottom seismographs have shown a high correlation between VLF events and reverse faulting in a shallow accretionary prism (Obana and Kodaira 2009). VLF events are characterized by a dominant frequency range of 2 to 8 Hz and a duration ranging from tens of seconds to a few minutes (Obana and Kodaira 2009). Thus, slip event in a shallow part of subduction zone shows great variation, from fast slip to long-duration slow slip.

The Nankai Trough is a subduction zone where the Philippine Sea Plate is subducting beneath the Eurasian Plate at a rate of 4.1 to 6.5 cm/year (Heki and Miyazaki 2001) (Figure 1). The age of subducting plate is relatively young (27 to 15 Ma; Okino et al. 1999), and a constant supply of hemipelagic sediments and trench-fill deposits on the oceanic plate is now building an accretionary prism at the margin. Along the subducting plate boundary, destructive earthquakes and tsunamis repeatedly occur (Ando 1975). Especially, a major splay fault, called ‘megasplay’ (Tobin and Kinoshita 2006), is considered as one of the primary coseismic faults that may have contributed to generating large historic earthquakes and tsunamis (Kikuchi et al. 2003; Baba et al. 2006; Tanioka and Satake 2001; Ichinose et al. 2003; Baba and Cummins 2005; Moore et al. 2007). In the Kumano Basin area off the Kii Peninsula, the megasplay fault branches from the plate boundary décollement at approximately 50 km landward from the trench at depth of approximately 10 km (Park et al. 2002).

Figure 1

Locations of the NanTroSEIZE drilling sites (a, b) and Distribution of vitrinite reflectance overlaid on photographs of slab core samples (c). (a) Regional tectonic setting of the IODP NanTroSEIZE drilling transect and rupture of the 1944 Tonankai earthquake (after Moore et al. 2007 and Kikuchi et al. 2003). Black line shows the location of the cross-section in (b). Abbreviations in inset: PP, Pacific Plate; PSP, Philippine Sea Plate; EP, Eurasian Plate; NAP, North American Plate. (b) Seismic profile along the NanTroSEIZE drilling transect (after Moore et al. 2007). VE, vertical exaggeration. (c) Photographs of the core sample that intersects the megasplay fault at Site C0004 (left) and of the frontal décollement (right) (after Sakaguchi et al. 2011). Black dots on the photographs indicate vitrinite particles measured by Sakaguchi et al. (2011) and color indicates Ro, reflectance of vitrinite. White dotted line indicates the location of a slip zone.

The IODP Nankai Trough Seismogenic Zone Experiment (NanTroSEIZE) Expedition 316 drilled several holes in the shallow parts of the megasplay fault and frontal décollement (Kimura et al. 2008; Kinoshita et al. 2009) (Figure 1). At Site C0004, the megasplay fault was recovered at 271 m CSF (core depth below the seafloor). The fault consists of microbreccia and an approximately 1-cm-thick narrow dark layer that dips approximately 50° from vertical axis of core. The frontal décollement was drilled at the toe of the accretionary prism (Site C0007). The décollement recovered at depth of 438 m CSF exhibits a 3-mm-thick dark layer which sharply separates well-brecciated hemipelagic mudstone above and unbroken hemipelagic mudstone below. An age reversal was found across this dark layer (Kinoshita et al. 2009).

To reveal the slip behavior, Hirono et al. (2009) investigated fault rock properties such as saturation magnetization or trace element concentrations across the megasplay fault and found no marked anomaly at the gouge zone. From the results, they suggested that the temperature in the gouge zone did not reach 400°C and that high-temperature (>300°C) water-rock interaction did not occur. Additionally, Hirono et al. (2014) revaluated trace element and isotope compositions of the gouge zone, and they revised the reached temperature downward (>250°C) from the absence of fluid rock interaction. In contrast, Yamaguchi et al. (2011) examined the major element composition and mineralogy of fault rocks along the megasplay fault and reported a progress of transition reaction from smectite to illite, which may be related to a temperature rise during frictional slip. Sakaguchi et al. (2011) measured vitrinite reflectance on polished core samples and detected an approximately 2-cm-thick positive reflectance anomaly along the megasplay fault and the frontal décollement (Figure 1c). They interpreted the anomaly as an evidence of coseismic frictional heat based on a simple assumption that broadening of the anomaly is proportional to the thermal diffusion length. In this study, the vitrinite data was more precisely processed by dynamic heat calculations to clarify the detailed process of frictional slip on the megathrust. Fulton and Harris (2012) proposed an analytical method to determine the slip velocity and distance based on vitrinite reflectance measured by Sakaguchi et al. (2011) in the Nankai Trough. Vitrinite is one of organic materials, whose reflectance has been widely used as a geothermometer. Fulton and Harris (2012) examined the vitrinite reflectance anomaly in and around the slip zone by combining numerical analyses of frictional heating and thermal diffusion. They reported that long slip duration (approximately 1,000 s) and large displacement (10 cm) can explain the measured reflectance anomaly along the megathrusts. Their attempt is appropriate to determine the slip velocity and slip duration (and slip distance). However, their analyses could be inaccurate due to lack of quantitative fitting calculation. In addition, they used the vitrinite reflectance inside of the gouge zone, which is thought to be excited by mechanochemical effect (Kitamura et al. 2012). Therefore, we improve the numerical fitting procedure of vitrinite reflectance data to derive more reliable slip parameters. The fitting is introduced only to the outside of the slip gouge zone to avoid shear influence. This is the first time to quantitatively evaluate the slip parameters on the megathrust in the Nankai Trough. To check the validity of the procedure, we compare the obtained slip parameters with those from other chemical approaches. Finally, we summarize the results and overview a variation of slip behaviors in shallow subduction zones.

Methods

Method to estimate slip parameters based on vitrinite reflectance

Vitrinite is one of the major organic materials in terrigenous sedimentary rocks. Optical reflectance of the vitrinite increases with its thermal maturation due to irreversible decompression of organic material. Therefore, the reflector has been widely applied as a thermometer which records the peak temperature the material experienced (e.g., Waples 1981). The vitrinite reflectance geothermometer has been used to detect the frictional heat both in natural fault zones (Bustin 1983; O’Hara 2004; Suchy et al. 1997; Sakaguchi et al. 2007; Sakaguchi et al. 2011; Maekawa et al., 2014) and in laboratory experiments of shear friction (O’Hara et al. 2006; Kitamura et al. 2012). For example, Sakaguchi et al. (2007) measured the vitrinite reflectance of natural fault rock samples and revealed a marked reflectance anomaly at a slip zone than the background values. O’Hara et al. (2006) conducted high-velocity slip experiments and confirmed an increase in vitrinite reflectance even at short duration of heating on an earthquake timescale.

Kinetic models of the vitrinite maturation and increasing in reflectance have been extensively discussed. Sweeney and Burnham (1990) constructed a kinetic model employing a series of first-order kinetics that describes parallel chemical reactions associated with the progress of vitrinite thermal maturation (Easy%Ro). Suzuki et al. (1993) proposed an improved kinetic model based on single activation energy, which is more suitable for numerical treatment (SIMPLE-Ro). In their model, the vitrinite reflectance R _o is expressed as a function of the fraction of the reacted material, F _c, as follows:

$$ {R}_o= \exp \left[ \ln \left({R}_{oo}\right)+3.7{F}_c\right], $$

(1)

where R _oo is the initial reflectance of vitrinite. F _c can be written as follows:

$$ {F}_C=0.85-0.85 \exp \left[-A \exp \left(\frac{-E}{RT}\varDelta t\right)\right], $$

(2)

where A is a frequency factor (1.0 × 10¹³/s), R is a gas constant, and E is the activation energy which can be represented as:

$$ E=40.7 \ln \left({R}_o\right)+{227}_{\left[\mathrm{kJ}/\mathrm{mol}\right].} $$

(3)

During a slip (0 ≤ t ≤ t _r , where t = time after the beginning of slip and t _r  = slip duration at a given area of fault plane; rise time), the temperature at a distance of x cm from the center of the slip zone, T(x,t), can be calculated from the amount of frictional heat and one-dimension thermal diffusion, as follows:

$$ \frac{\partial T}{\partial t}-\kappa \frac{\partial^2T}{\partial {x}^2}-\frac{\dot{Q}}{C}=0\ \mathrm{f}\mathrm{o}\mathrm{r} - \frac{w}{2}\le x\le \frac{w}{2} $$

(4)

and

$$ \frac{\partial T}{\partial t}-\kappa \frac{\partial^2T}{\partial {x}^2}=0\ \mathrm{f}\mathrm{o}\mathrm{r}\ x<-\frac{w}{2},\frac{w}{2}<x $$

(5)

where $\dot{Q}$ is the rate of heat generation, C is the heat capacity of the slip zone, κ is the thermal diffusivity, and w is the thickness of the heat generation zone in the slip zone. For simplicity, w equals to thickness of the slip zone (gouge zone), and the energy taken up by chemical reactions and heat transfer due to fluid flow are ignored. $\dot{Q}$ equals to the product of frictional coefficient μ, effective normal stress ${\sigma}_n^{\prime }$, and slip velocity v $\left(\dot{Q} = \mu {\sigma}_n^{\prime }v\right)$. C is calculated from bulk heat capacity Cpbulk, bulk density ρ, and w (C = Cpbulk ρ w). Equations 4 and 5 represent the temperature within and outside of the slip zone, respectively. The third term on the left side of Equation 4 is the term of increasing temperature. The temperature after the slip event is also expressed in Equations 4 and 5, by considering that heat generation $\dot{Q}$ is zero. These equations include eight independent parameters: κ, μ, ${\sigma}_n^{\prime }$, v, Cpbulk, ρ, w, and t _r . Among them, seven parameters were determined by physical property measurements and sample observation. Frictional coefficient was determined by a friction experiment. Thermophysical properties were also measured by a calorimeter with laser flash technique, and the width of slip zone was measured by photo scale observation. ${\sigma}_n^{\prime }$ was estimated by considering solid density, porosity (ϕ), and depth of the fault plane (z) under the assumption of hydrostatic pore fluid pressure:

$$ {\sigma}_n^{\prime }=g{\displaystyle \underset{0}{\overset{z}{\int }}{\rho}_{(z)}dz} $$

(6)

where g is the gravitational acceleration. The parameters adopted for calculation are listed in Table 1. Consequently, the unknown parameters are only $\dot{Q}$ and t _r in these equations. v can be estimated by substituting μ and ${\sigma}_n^{\prime }$ into $\dot{Q}$.

Table 1 Parameters for calculating the temperature-time profile and vitrinite maturation

Combining vitrinite kinetic equations, Equations 1, 2 and 3, and temperature equations, Equations 4 and 5, spatial evolution of the vitrinite reflectance driven by frictional heating can be calculated. For numerical calculation, a finite difference method was employed with time and distance steps of 0.01 s and 0.02 cm, respectively. Because vitrinite maturation continues even at a cooling stage, the simulation was continued until the temperature goes down to 150°C. At 150°C, 2 × 10⁶ s is needed to increase the reflectance by 0.1%; thus, this temperature is adopted as a cutoff value of calculation.

The distributions of temperature and vitrinite reflectance both depend only on $\dot{Q}$ and t _r . The slip parameters can be determined by numerical fitting of the measured distribution of the vitrinite reflectance at various $\dot{Q}$ and t _r . However, the data to be fitted should be defined carefully. Based on a high-speed friction experiment and vitrinite reflectance measurement of experimental products, Kitamura et al. (2012) revealed a mechanochemical effect that may also accelerate vitrinite maturation. Accordingly, vitrinite reflectance values in the slip zone might be also influenced by this mechanical effect. Therefore, fitting was applied only to reflectance profile at the outside of the slip zone.

To evaluate the fitting accuracy quantitatively, fraction of residue sum of square and deviation sum of square (RSS) is introduced:

$$ \mathrm{R}\mathrm{S}\mathrm{S}=\frac{{\displaystyle \sum_i^n{\left({y}_i-{f}_i\right)}^2}}{{\displaystyle \sum_i^n{\left({y}_i-\overline{y}\right)}^2}} $$

(7)

where y _i is the measured value of vitrinite reflectance, f _i is the calculated reflectance value from Equations 1, 2, 3, 4 and 5, $\overline{y}$ is an average of measured values, and n is the number of observational points, respectively. The numerator in Equation 7 is a sum of squared errors between the observed and calculated values. The denominator of Equation 7 is a sum of square deviation. Figure 2 shows a schematic illustration of the relation between y _i, $\overline{y}$, and f _i. RSS represents the fraction of squared errors from the measured values. This indicator should be 0 when the measured data are completely reproduced by the calculation. When the residue indicators yield minimum value, the fitting parameters are expected as plausible fault slip parameters. Fitting is performed by grid search method to find optimal combination of $\dot{Q}$ and t _r . RSS is calculated as wide range of parameters set as possible not to converge at local minimum.

Figure 2

Schematic figure of the calculation method to obtain residue. Randomly measured reflectance value (measured Ro) is laid in the direction normal to the fault plane. The reflection values are averaged within each 2-mm interval (black dots), and its error bar represents standard error. y and f represent the value of vitrinite reflectance obtained from measurement and simulation, respectively. y is the averaged value of measured reflectance.

Slip parameters estimation on the Nankai trough megathrusts

The vitrinite reflectance measurement by Sakaguchi et al. (2011) shows that vitrinite particles closer to the slip zone take higher reflectance (Figure 1c). Firstly, the data are arranged as a function of distance from the center of the slip zone. Then, we averaged the reflectance values at every 2-mm interval (Figure 3). The measured reflectance values are scattered (Sakaguchi et al. 2011). This data handling reduces noise of data and provides a clear profile of reflectance distribution. However, it is noted that reflectance inside of the slip zones would be influenced by mechanochemical effects in addition to temperature rise (Kitamura et al. 2012). To eliminate this effect, the fitting is applied to the ranges 0.6 to 3.0 cm below the slip zone of the megasplay fault and 0.4 to 3.0 cm below the frontal décollement. The reason for fitting the data only on the footwall side is that these faults are reverse faults, and thus, it is difficult to determine the exact depth condition when the reflectance on the hanging wall side was acquired. In fact, the vitrinite reflectance in the hanging wall sides indicates higher value and different reflectance peak broadness, for both of the megathrusts. This might be because hanging wall had experienced higher temperature event at the deeper part before they have been lifted up to the current depth by thrust faulting. In any case, the reflectance distribution in the hanging wall would contain some complexity comparing to footwall. Thus, we do not discuss thermal anomalies on the hanging wall side to avoid potential complexities.

Figure 3

Calculated RSSs for given Q and t _r of the megasplay fault (C0004: a, a’) and for the frontal decollement (C0007: b, b’).

The fitting accuracy was estimated by RSS calculations (Figure 3a,b). In these numerical analyses, $\dot{Q}$ and t _r were varied in the range of 10³ to 10⁷ J/s/m² (10^0.04 J/s/m² step), and 10⁰ to 10^5.5 s (10^0.05 s step), respectively. The minimum RSS was obtained when $\dot{Q} = 16,600\ \mathrm{J}/\mathrm{s}/{\mathrm{m}}^2$ and t _r  = 6,250 s for the megasplay and $\dot{Q} = 23,200\ \mathrm{J}/\mathrm{s}/{\mathrm{m}}^2$ and t _r  = 2,350 s for the frontal décollement. Figure 4 shows the calculated temperature distribution (upper panels) and corresponding vitrinite reflectance distribution plotted on the observed data (lower panels) for the megasplay fault and the frontal décollement. Temperature rises with time until t _r ; then, it decreases by thermal diffusion. For both of the megasplay and frontal décollement, estimated t _r reaches as long as 1,000 s, so temperature in the host rock is also increasing even during slipping. The calculated maximum temperature is approximately 340°C for both faults. Considering the normal stress acting on the fault plane $\left({\sigma}_n^{\prime}\right)$ estimated from Equation 6, the slip velocity and total slip distance can be calculated from $\dot{Q}$ and t _r . Under the hydrostatic condition, ${\sigma}_n^{\prime }$ was estimated as 2.16 MPa for the megasplay and 4.20 MPa for the frontal décollement, based on the shipboard porosity data (Kinoshita et al. 2009). Consequently, we obtained the slip velocity and slip distance of 0.91 cm/s and 56.7 m for the megasplay fault and 0.65 cm/s and 15.3 m for the frontal décollement, respectively.

Figure 4

Temperature and vitrinite reflectance measured and calculated by the most appropriate parameters. The shaded area represents the slip zone. Upper panels show the temperature distribution at given times, and color corresponds to the heating time. Lower panels show the distribution of measured (points with error bars) and fitted (line) vitrinite reflectance. Measured points framed by red line were used for RSS calculations. Letters to the left of each panel are calculated parameters for the best fitting.

These results were obtained by the fitting of reflectance data which was averaged at 2-mm intervals. This averaging might make a small reflectance-peak inconspicuous and could broaden the range of reflectance anomaly around the faults. Additionally, several data points show lower reflectance values than the fitted curve, and the differences between the data values and fitted line exceed the ranges of their error bars (Figure 4: two points within approximately 1.5 cm from the localized slip zone in C0004 and one point around approximately 1 cm in C0007). These points below the fitted curve might represent that the heat signature around the localized slip zone would be narrower and steeper than we fitted. They might also indicate the existence of another reflectance peak at 1.5 to 2.0 cm from both of the faults, representing slip within the brecciated zone. The reflectance anomaly, however, is considered to be one broad peak in our method, which fits one calculated curve to the reflectance anomaly. Therefore, the slip parameters might be overestimated. Although it is difficult to investigate the small peak due to the limited data region or a lack of detailed description of microstructure along the faults, the calculated parameters could be implicated as upper limit of slip estimation on each fault. The detailed structural observation is needed to define specifically the broadness of the reflectance anomaly.

The parameters used in these calculations are well constrained by shipboard measurements or onshore experiments (Table 1), except for the effective normal stress and thickness of heat generation zone, w. In these calculations, shear stress (τ) is evaluated by assuming hydrostatic condition and a constant frictional coefficient obtained by friction experiment. These assumptions probably provide a maximum estimate for τ. Using an apparent thickness of the gouge zone (1.0 and 0.3 cm for the megasplay and the frontal décollement, respectively) is the maximum estimate for w because a heat generation zone in each slip event will be narrower than slip zone when slip localization occurs. Therefore, additional analyses should be performed to assess the sensitivity of the results with respect to τ and w. Figure 5 shows a range of slip velocity (v) and slip distance (D) as a function of shear stress. As mentioned above, heat generation rate is expressed by effective normal stress, frictional coefficient μ and slip velocity v $\left(\dot{Q} = \mu {\sigma}_n^{\prime }v = \tau v\right)$. Because $\dot{Q}$ can be determined by minimum residual analysis, τ and μ have an inverse relation with v (Figure 5a). The slip distance is obtained by the product of v and t _r , implying that the slip distance is also inversely proportional to τ. Figure 5b shows that, in a given range of τ, the slip velocity is cm/s order and that the slip distance is larger than 10 m for both the megasplay and frontal décollement.

Figure 5

Sensitivity of slip velocity and slip distance to the amount of shear stress. Sensitivity of slip velocity (v) and slip distance (D) as a function of shear stress (τ); (a) the megasplay fault and (b) the frontal decollement. Black lines indicate the condition of hydrostatic stress, assuming that horizontal stress is equal to vertical stress.

The thickness of the heat generation zone, w, has a strong influence on the calculated temperature; however, it is difficult to quantify this parameter from a geological observation (Hamada et al. 2009). Here, we assumed plausible values for w (0.1 and 0.5 cm) and assessed the sensitivity of the calculation results for the megasplay fault with respect to w. This is equivalent to examine an effect of localized slip for the vitrinite reflectance distribution. Figure 6 shows the results of sensitivity analyses at constant shear stress, slip velocity, and rise time. The distribution of reflectance is well fitted by the same slip parameters for different values of w, indicating that the reflectance distribution outside the slip zone is almost independent on w. Localized slip (thinner w) results in higher temperature because of a smaller heat capacity than that for thicker w (C = Cpbulk ρ w). However, a large temperature difference between a slip zone and a host rock also entails more rapid heat conduction. Thus, the difference of w has little effect on dynamic temperature change as well as the distribution of vitrinite reflectance outside of the slip zone.

Figure 6

Sensitivity of slip velocity and slip distance to the thickness of the heat generation zone for C0004. Results of a sensitivity analysis of slip zone thickness (w) along the megasplay. Shaded area represents the slip zone. (a) w = 0.5 cm. (b) w = 0.1 cm. Light gray shaded area shows the gouge zone (1 cm).

Another factor to be noticed is a cumulative effect of repetitive slips. Because the megasplay and frontal décollement might have slipped repeatedly, the heat signals should reflect the results of cumulative thermal events. Fulton and Harris (2012) discussed about this repeated heating effect to the vitrinite maturation. They argued that the vitrinite maturation may proceed with approximately 100 slips because of the necessity to satisfy the upper limit of the heating temperature (<300°C; Hirono et al. 2009). According to their sensitivity analyses of slip parameters, matching Ro% distribution, to cumulative effects, the estimated slip parameters in the condition of single slip are different from the ones in the case of 1,000 or 10,000 slips. The estimated slip distance of 1 slip event in the 10,000 slips, however, is the same in the order of magnitude as that of the single slip case. This is because the maturation of vitrinite notably depends on the maximum temperature (kinetic theory, Suzuki et al. 1993, and comparison with other thermometer, Barker and Goldstein 1990). It would imply that one maximum heating event has more significant impact than repetitive small slips and mostly overwrite the record of lower thermal events. In addition, the cumulative distance would be approximately 40 or approximately 400 km with 1,000 or 10,000 slips to explain the measured vitrinite maturation (Fulton and Harris 2012). The large accumulative slip may not be a feasible value. Because, the background temperature has little difference between the hanging wall and footwall side and the fact would show the slip zones have not had large displacement. When the case of ten slips is assumed to match Ro%, the cumulative distance (approximately 400 m) may be acceptable and slip parameters are comparable to ones in the case of one slip. Consequently, although slip parameters should be overestimated in the case of assuming single slip, unconstrained parameters such as shear stress have much greater influence on the slip parameter estimation than cumulative slip effect (Figure 5). Here, we considered that this cumulative effect is not very essential for slip parameter estimation of these megathrusts, and we regard that the calculated slip parameters estimated assuming single slip as upper limit of maximum values related to the largest slip event of the past.

Results and discussion

Based on the vitrinite reflectance data on the two faults, we obtained the slip velocity and distance for each fault. Among the megasplay fault, the slip behaviors were investigated by using other chemical approach (Yamaguchi et al. 2011; Hirono et al. 2009; Hirono et al. 2014). Here, we examine the validity of the results by checking the consistency with other constraints provided by the previous works.

Yamaguchi et al. (2011) performed elemental analysis of core sample along the megasplay. X-ray fluorescence (XRF) core-imaging scanner was employed for mapping of major element. They reported a progress of transition reaction from smectite to illite through X-ray diffraction (XRD) analyses for mineralogical investigation. Though the clay transition could not be converted to temperature quantitatively due to unknown of volume fraction of the minerals, they concluded that the smectite to illite transition would be enhanced by a temperature rise during frictional slip along the megasplay. Hirono et al. (2009) investigated fluid-mobile trace element (Sr, Cs, Rb, and Li) concentrations, Sr isotopes, magnetic properties of minerals, organic/inorganic carbon content, and Raman spectroscopy of carbonaceous material on discrete samples recovered from the splay fault in the Nankai Trough and examined whether the fault has experienced increasing temperature due to frictional heating along the fault zone. Concentrations of fluid-mobile trace elements and Sr isotope ratios are very sensitive to fluid-rock interaction at high temperature, making them strong indicators of frictional heat in the slip zone. This geochemical temperature proxy has been applied to not only the cored fault gouge (Ishikawa et al. 2008; Hirono et al. 2009) but also on-land fossil faults (Hamada et al. 2011 and Ishikawa et al. 2014). Paramagnetic susceptibility, χpara, and saturation magnetization, M _s , are related to abundances of paramagnetic and ferrimagnetic minerals, respectively. Paramagnetic minerals can be decomposed at high temperature, and they would change to ferrimagnetic minerals such as magnetite (Mishima et al. 2009). High temperature also induces thermal decomposition of carbonate minerals. Comparison of carbonate content between the fault zone and the host rock is one of tools to detect frictional heating. These analyses were also applied to the Taiwan Chelungpu fault samples, which was collected by Taiwan Chelungpu fault Drilling Project and have succeeded in detecting heat signals within the fault gouge zone (e.g., Ishikawa et al. 2008; Mishima et al. 2006; Hirono et al. 2007). In addition, Hirono et al. (2014) performed X-ray diffraction analyses to detect clay mineral transition and reexamined trace-element and isotope compositions of the fault gouge in detail. However, all of these measurements failed to detect marked anomaly at the megasplay fault. From the results, Hirono et al. (2009) and Hirono et al. (2014) suggested that the temperature in the slip zone did not reach 400°C and that high-temperature (>250°C) water-rock interaction did not occur. This temperature limitation are estimated by chemical measurements of discrete samples, and the data density is different from map analyses of major elements (Yamaguchi et al. 2011) and vitrinite Ro data (Sakaguchi et al. 2011) on split core surfaces. This temperature threshold, however, would be one of the criteria of temperature estimation.

In our numerical analyses, dynamic change in temperature during slip is also calculated to obtain distribution of vitrinite reflectance. Tmax, (temperature at t _r ), was calculated to be 331°C in average for the slip zone of the megasplay fault (Figure 4a) with 16,600 J/s/m² of $\dot{Q}$ and 6,250 s of t _r . The calculated Tmax is lower than the value (approximately 400°C) estimated by Sakaguchi et al. (2011). Although this temperature condition is higher than the constraint by Hirono et al. (2009) and Hirono et al. (2014), this can be probably explained by lower fluid-rock ratio in the slip zone. Trace element behavior is highly dependent on fluid-rock ratio, and in general, higher water ratio tends to facilitate element transfer between fluid and solid (You et al. 1996). In the Taiwan Chelungpu fault, where element transfer was detected in the fault zone, the fluid to rock ratio was estimated as high as 4 to 10 (Ishikawa et al. 2008). It is supposed that lower fluid-rock ratio (for example, <2) might inhibit efficient fluid-rock interaction in the megasplay fault. This absence of the fluid-rock interaction suggests that infiltration of high-temperature fluid may not be a heat source substitute for the frictional heating. Vitrinite reflectance might show similar distribution to a measured one, when the hot fluid flows at a similar heat generation rate to frictional heating. In this case, however, the hot fluid may migrate with long duration (approximately 1,000 s) inducing fluid-rock interaction and will mark a record in the solid grains inside the fault, such as trace element signals. Thus, the rise in temperature can be considered as provided by frictional heat. Estimated slip velocity by this research is approximately 1 cm/s, which is slow relative to a normal earthquake. This slip velocity is also consistent with the previous statement that high velocity (approximately 1 m/s) slip did not occur upon the megasplay fault (Hirono et al. 2009), and we concluded estimated slip parameters as feasible.

Kinetic approach requires an appropriate kinetic expression in principal. Ideally, such kinetic equation should be experimentally determined by using individual fault rock samples (Hamada et al. 2009). It is necessary to investigate whether the reaction kinetics is applicable to heating for short duration. However, it is practically difficult to completely reproduce the condition where faulting occurs, even though a shallow depth condition (<3 km) in subduction zones (high confining pressure, low-oxygen fugacity, and fluid-saturated condition). O’Hara (2004) introduced the Easy%R, to calculate Ro change in the rapid cooling system (approximately 1.0°C). They compared the simulation results with maturation experiment in laboratory (Bustin 1983) and with other Ro-temperature conversion method based on independent geologic data (Barker and Goldstein 1990). They indicated the consistency between Easy%R and rapid (approximately 1.3°C /s) heating experiment. On the other hand, uncertainty of ±50°C in Ro temperature converting is found by the comparison between simulation and other temperature estimation method. We calculated the maximum temperature as approximately 330°C in the megasplay. If the ±50°C error is applied, this estimated value has a variation from 280°C to 380°C, and it would be still compatible with the suggested temperature limitation (approximately 400°C; Hirono et al. 2009). It is difficult to evaluate quantitatively how the temperature error affects the estimated slip parameters. However, slip parameters would not be influenced very much by this temperature error. Even if the required temperature for vitrinite maturation varied within 50°C, heat duration would be still some thousands of seconds to achieve the temperature (Figure 4). Thus, the temperature error does not have a significant effect to obtain the characteristic of slip, slow slip rate of approximately 1 cm/s and long duration of 1,000 s s. However, chemical kinetics for less maturated vitrinite (Ro < 0.5) in rapid heating has not been verified by any experiments so far. Through experiments specialized in rapid heating, development of the chemical kinetics of vitrinite maturation is necessary for more detailed analysis.

It is also possible that reaction mechanism changes drastically in fault zones because of the shearing effect. If such additional effect other than temperature fosters the vitrinite reaction (higher apparent frequency factor and/or lower activation energy), it would lead to a higher reflectance with smaller amount of heat. In this case, we potentially overestimated the slip parameters. Here, we focused only on data outside of the slip zone in order to eliminate the shearing effect. However, if there is no choice other than using an existing kinetics, one should be careful that the numerical analysis could include a certain value of error.

Conclusions

In this study, we examined slip behaviors on the shallow subduction zone faults (megasplay fault and frontal décollement in the Nankai trough), based on numerical analysis of vitrinite reflectance. The two faults exhibit almost similar slip properties: mm/s to cm/s of slip velocities, approximately 1,000 s of slip duration and meters to several tens of meters of displacements. These characteristics imply that the slip events on these faults are extremely slow, long duration, and large displacement compared to those for normal large earthquake in subduction zone.

The slip parameters obtained in this study can give a constraint on other numerical analyses such as tsunami simulation. This estimation also predicts that slow slip velocity, long risetime, and very large total displacement would be common features of slip behaviors in shallow subduction zones. The obtained slip parameters might be consistent with rapid afterslip, but less likely a shallow slow slip earthquake considering their large displacement. Such anomalous slip has never been recorded by seismic observations, but is probably not special at least in a geologic timescale. As mentioned above, the thermal indicator reflects the most dramatic slip event occurred on the relevant fault. A possibility of such slip (cm/s order slip rate and 10,00s s order duration) should not be dismissed just because of lack of observation. The slip zone core samples we focused, however, were small fractions of each huge megathrust, and it is possible that estimated slip parameters do not represent behavior of the whole slip area. Our estimation also would contain some uncertainty owing to data resolution, frictional parameters (shear stress and friction coefficient), and to cumulative effect by repeated slip. The estimated parameters would correspond to values of the upper limit of the most dramatic slip as mentioned in the Slip parameters estimation on the Nankai trough megathrusts and Results and discussion sections. Therefore, similar approaches for on-land fault outcrops with detailed structural observation would confirm that such slow and long-duration slip universally occurs. If such slips truly existed, how do they occur and propagate? To answer this question, the numerical analyses method is needed to be applied to faults in various parts of subduction zone, and to faults whose slip parameters are deduced by seismic observations, one candidate is a deeper part of the megasplay fault in the Nankai trough, which is planned to be drilled by the forthcoming NanTroSEIZE.