1 Introduction

The Hyuga-nada region is located in the southwestern part of the Nankai Trough in the Pacific Ocean (Fig. 1A). M7 class interplate earthquakes repeatedly occur because of the subduction of the Philippine Sea Plate beneath the overriding Eurasian Plate or Amur Plate (e.g., Heki et al., 1999; Seno et al., 1993). In the recent past, thrust-type earthquakes such as the 1968 earthquake (Mw 7.5) and the 1996 earthquakes (Mw 6.8 and Mw 6.7) occurred in this subduction zone (Yamashita et al., 2012). As indicated by the distribution of slow earthquakes in southwest Japan (Obara & Kato, 2016), the shallow and deep parts of the plate boundary are seismically active with respect to slow earthquakes in Hyuga-nada. Broadband seismograms were analyzed to detect shallow very-low-frequency earthquakes (Asano et al., 2015); very-low-frequency events were also observed around the up-dip portion of the seismogenic zone of the plate boundary in this region (Obara & Ito, 2005). Additionally, ocean bottom monitoring for seismicity off the southern Kyusyu region recorded a complete episode of low-frequency tremor lasting for 1 month in 2013, which was associated with very-low-frequency earthquake activity at the shallow plate interface (Yamashita et al., 2015). Similar seismic activities were also seen in New Zealand. From Wallace et al. (2016), seafloor deformation observations revealed direct evidence for slow slip events occurring close to the trench (within 2 km of the seafloor). Long-term slow slip events and frequent down-dip of the Hyuga-nada stable sliding zone were also detected (Takagi et al., 2019). The seismic activity in Hyuga-nada over the years is shown in Fig. 1B.

Fig. 1
figure 1

A Map of Japan and Hyuga-nada subduction zone with depth contour interval of 2000 m. B Seismic activity in Hyuga-nada in southwestern part of Nankai Trough. Blue stars indicate interplate earthquakes greater than M6.5 since 1923 detected by the Japan Meteorological Agency. Orange circle denotes long-term slow slip events (L-SSEs) with cumulative slip of > 40 cm from 1996 to 2017 from Takagi et al. (2019). Green area shows active source of shallow low-frequency tremors from Yamashita et al. (2015). Purple dashed ellipse shows tsunami source from Hatori (1985)

Full size image

The largest earthquake in history in the Hyuga-nada subduction zone was the 1662 Hyuga-nada earthquake (M7.25–7.75), which occurred off Miyazaki Prefecture in the southeastern part of the Kyushu region (Usami, 2003) and generated a tsunami (Watanabe, 1998). Tens of people died, and a few thousand structures were damaged by the earthquake and the tsunami (Watanabe, 1998). The tsunami heights were estimated to be least 4–5 m along the coast in Miyazaki City, as estimated from historical records (Hatori, 1985). The 1662 tsunami was much larger in scale than those usually generated by the repeated M7 class interplate earthquakes that occur in this subduction zone. From the damage to structures that occurred along the coast of Miyazaki Prefecture, it can be inferred that the JMA seismic intensity was high, up to a maximum of 6. Historical documents indicate local subsidence of about 1 m and changes at river currents of Oyodo and Kaeda in the Miyazaki plain (Hatori, 1985). The inundation and runup heights of the tsunami that remained as historical documents were also surveyed and re-estimated by Tsuji et al. (2018).

Tsunami deposit investigation occurred to elucidate the inundation area of the 1662 tsunami by previous studies. However, these deposits are recognized only in Miyazaki plain where much subsidence occurred by the 1662 Hyuga-nada earthquake. Therefore, these geological records are control points but are not enough to determine tsunami inundation area accurately.

In this study, we hypothesize that the unusually large tsunami of 1662 was caused by the coseismic slipping of the focal area of not only the repeated M7 class interplate earthquake but also shallow slow earthquakes. We investigated the tsunami deposits in the lowlands along the coast of Miyazaki Prefecture and compared surveyed tsunami deposits and computed inundation area using a constructed fault model.

2 Fault Model

2.1 Fault Models in Previous Studies

Hatori (1985) estimated the source of the tsunami that followed the 1662 Hyuga-nada earthquake considering tsunami heights of 4–5 m and a subsidence of approximately 1 m at the Miyazaki plain, derived from historical documents. The estimated tsunami source was located offshore near the coast included land area (Fig. 1B). Matsu’ura et al. (2003) estimated the source area of the 1662 Hyuga-nada earthquake, considering seismic intensities in the Kyusyu region. Seismic intensities were re-examined based on the structural damage caused by earthquakes. The earthquake was characterized as an interplate earthquake with estimated magnitude of M7.2–7.5. The local subsidence of about 1 m at the Miyazaki plain was interpreted to be due to coseismic deformation or land subsidence and liquefaction resulting from strong motion. Nakata et al. (2012) deal with a fault model of the 1662 Hyuga-nada earthquake as the interplate earthquake of M7.6 in numerical simulation of the seismic cycle in Hyuga-nada.

Thus, details including the rupture area and type of the 1662 Hyuga-nada earthquake are not clear. A fault model of this earthquake that satisfies the conditions of constraint, namely a subsidence of about 1 m and tsunami heights of 4–5 m at the Miyazaki plain, has not been proposed before.

2.2 Relation Between Slow Earthquake and Great Earthquake of Tohoku-Oki

The relationship between the slow earthquake and the 2011 great earthquake of Tohoku-oki is key to estimating the fault model of the 1662 Hyuga-nada earthquake. Near the Japan Trench, active seismicity corresponding to slow earthquakes in the shallow part of the plate interface is observed. Before the 2011 Tohoku-oki earthquake, slow slip events occurred off the Tohoku region (Ito et al., 2013). Two slow slip events in 2008 (Mw 6.8) and 2011 (Mw 7.0) were recorded using ocean-bottom pressure gauges and an onshore volumetric strainmeter. The estimated slow slip fault from the model corresponding to Mw 7.0 is located within the large coseismic slip area exceeding 30 m of the 2011 Tohoku-oki earthquake. The observed slow slip events before the 2011 Tohoku-oki earthquake possibly played a key role in causing the huge coseismic slip near the Japan Trench. This is because that clear velocity weakening behavior occurred upon the arrival of the positive velocity step during propagation of the nearby large coseismic rupture into the slow slip source area. Episodic slow slip and seismic behavior occurred on the same portions of the megathrust fault, suggesting that the faults that undergo slip in slow slip events can also rupture seismically. Therefore, it can be inferred that slow slip events precede moderate or large interplate earthquakes.

The key to understanding the mechanism of tsunami earthquakes could be in slow slip events inducing a relatively large earthquake while the slow slip event is going on. The effect of slow slip events facilitating weakening could also indicate a mechanism of tsunami earthquakes (Ito et al., 2017). Frictional instability may have been enhanced near the Japan Trench because of frictional slip weakening, leading to the occurrence of tsunami earthquakes.

2.3 Estimation of Fault Model

Considering the relation of slow slip events and large earthquakes, we constructed the fault model including the parameters of the source region where there is active seismicity related to slow earthquakes in the shallow part of the plate boundary in Hyuga-nada. Our fault model of the 1662 Hyuga-nada earthquake was based on recent geophysical findings, and modern understanding is an improvement over the fault model of the interplate earthquake, obtained by changing the length (Fig. 2). Slips of subfault 1 and subfault 2 at the up-dip part of the plate interface are needed to generate the large tsunami. Slip of subfault 3 at the down-dip part of the plate interface is needed to generate the strong ground motion that could explain widespread structural damage along the coast of Miyazaki Prefecture. A length of 80 km and width of 25 km were set to subfault area. A strike of 212° was set in a direction perpendicular to the plate subducting angle calculated from the UNAVCO plate motion calculator (Kreemer et al., 2014). A rake of 90° was set assuming an earthquake of pure thrust type. From a plate boundary model obtained by compiling the structural data from offshore surveys (Nakanishi et al., 2018), fault model parameters were assumed as follows: dip of 8.1°, 9.2°, 16.9°, and the depth of the upper edge at the subfault of 11.4 km, 15.0 km, 19.1 km, corresponding to subfault number of 1, 2, 3, respectively. The locations of the three subfaults described above are shown in Fig. 2, and the fault model parameters are listed in Table 1.

Fig. 2
figure 2

The rupture area of the estimated fault model of 1662 Hyuga-nada earthquake is represented by a red rectangle. Area within green line denotes shallow tremor active area. Field survey point is indicated by pink dot. Depth contour interval is 1000 m

Full size image
Table 1 Estimated fault model parameters of 1662 Hyuga-nada earthquake
Full size table

3 Tsunami Deposit Investigation

As mentioned of introduction, some 1662 tsunami deposits have been reported in the Miyazaki plain (e.g., Kamataki et al., 2013, 2015; Niwa et al., 2020). Kamataki et al. (2015) reported 1662 coarse sandy deposits interbedded with fine sands at ~ 1.5 m depth from ground surface. Niwa et al. (2020) also reported that sedimentary environment changed from terrestrial to shallow marine at the boundary of 1.2–1.5 m below ground surface by changing geochemical proxy even in the area affected by 1662 subsidence, although tsunami deposits themselves were not founded.

Based on their results, we bored to a depth of 2 m using a hand coaler (sample length: 50 cm, diameter: 4 or 5 cm), peat sampler (sample length: 30 cm, diameter: 1.5 cm), and shovel to predictively confirm 1662 deposits at 62 survey sites in the lowlands along the eastern Kyushu (Miyazaki and Kagoshima Prefectures). The study sites, which were located several hundred meters to several kilometers from the present coastline at several meter elevation, were selected with reference to the tsunami inundation areas reported by Crisis Management Course of Miyazaki Prefecture (2013). As a result, we newly found a sandy event deposit interbedded with silt in Komei, Nichinan City (Fig. 2). The locality is found at the approximately 7 m present elevation and 450 m inland from the present coastline. This area was previously used as fields, but is now a wasteland.

Event deposits (bed up to approximately 10 cm thickness) consisted of well-sorted moderately rounded, very fine to fine grain-sized sand and show normal grading and locally parallel lamination at the top of the layer (Fig. 3A–C). Beach sands near the survey area were characterized by well-sorted and well-rounded medium sands, while river sands were mainly pebbles and cobbles with minor poorly sorted angular coarse sand particles (Fig. 3D–F). The beach sands also contain modern shell fragments (Fig. 4A, B), which are not recognized in river sands but in event deposits (Fig. 4C–F), although their species cannot be identified because of fragments. Based on distribution of surrounding basement rocks (Kimura, 1959), the Miocene turbidites (Miyazaki Group) were distributed around the excavation sites. Therefore, it was not considered that shell fragments could have been supplied from the landward side into the study sites. Based on the above clastic composition and roundness, the provenance of event deposits is not river sands, suggesting that event deposits were not formed by foods. It was difficult to identify whether event deposits were formed by tsunamis or storms. However, previous studies reported that storm deposits did not exist even in Miyazaki plain (Kamataki et al., 2013, 2015; Niwa et al., 2020), which is located in a neighboring study site and has higher potential to preserve storm deposits than this study site. Based on the provenance of event deposits, the distance from the coastline to the excavation sites, and elevation, it is reasonable that event deposits were formed by a large-scale tsunami event rather than a storm. According to the previous study (Kamataki et al., 2015), sedimentary rates of fine background sediments after 1662 were estimated as approximately 0.1–0.4 cm per 100 years. Thus, event deposits recognized in this study could have been formed by a tsunami approximately several hundred years ago, suggesting correlation with the 1662 tsunami event. We used this new information as a control point for tsunami simulation in the southern area, where there have been no reports of tsunami deposits.

Fig. 3
figure 3

Photographs of A core samples and B X-radiographs obtained in Komei, Miyazaki Prefecture. C Mean grain size distribution of event deposits and back ground silts. D Mean grain size distribution of beach and river sands. E Photograph of beach area. F Photograph of river area. GL Ground level

Full size image
Fig. 4
figure 4

Photomicrographs of A, B beach sand, C, D event deposit, and E, F river sand. A, C, and E: plane polarized light. B, D, and F: cross-polarized light. S shell fragment, Ns Neogene sedimentary rock, Pl plagioclase

Full size image

4 Numerical Simulation of Tsunami

To estimate the coseismic deformation, the vertical ocean bottom deformation was computed using the equation of Okada (1985) based on the fault model detailed in Sect. 2.3. The vertical maximum uplift of 2.64 m near the trench and minimum subsidence of − 0.46 m was computed. Assuming that the ocean bottom deformation is equal to the vertical ocean surface deformation, ocean surface deformation was used as the initial tsunami. To compute tsunami heights in the large computational area, linear shallow water equations were solved with a staggered leap-frog method (Satake, 1995), while in the small computational area to compute tsunami inundation, non-linear shallow water equations were solved using a finite difference scheme (Goto et al., 1997). In the numerical simulation of tsunami inundation, a moving boundary condition was used (Imamura, 1996). A Manning’s roughness coefficient of 0.025 m−1/3 s was assumed as the bottom friction term in the non-linear shallow water equations. The computed coseismic deformation was also considered at a small computational area in the tsunami inundation simulation. The duration was set to 60 s, computation time to 1 h, and output time to every 1 s in the simulation. The time step was set to 0.05 s in the tsunami inundation simulation to satisfy the Courant-Friedrichs-Lewy (CFL) stability condition. The bathymetry and topography are based on the topography data of the tsunami fault model (Cabinet Office, 2016) published by the Association for Promotion of Infrastructure Geospatial Information Distribution. The bathymetric and topographic grid with a spatial resolution of 3 arc-seconds interpolating the base data was used for computing the tsunami within the coordinates 130–134°E and 30–34°N. The bathymetric and topographic data with a grid interval of 1 arc-second interpolating the base data were used for tsunami inundation simulation at Komei in Nichinan City within 131.44–131.5025°E and 31.65–31.7125°N. In tsunami inundation simulation, artificial structures, such as embankments, have been removed according to the topography of the time.

5 Results

From the rupturing of the source area of slow earthquakes that occurred at the shallow part of the plate boundary in addition to the usual repeated M7 class-defined approximately M7.0–7.5 interplate earthquakes, it can be inferred that the 1662 Hyuga-nada earthquake is the largest of the Hyuga-nada earthquakes in historical records. In the estimated fault model, the rupture areas of subfault 1 and subfault 2 containing the source area of slow earthquakes together were needed to generate the large tsunami; a slip of 8 m at subfault 1 and a slip of 4 m at subfault 2 were needed to explain the tsunami heights along the coast in the Miyazaki plain in historical records and the tsunami inundation area at Komei. A length of 80 km was needed to explain the tsunami inundation area at Komei in the southern part of Miyazaki Prefecture, and a slip of 2 m at subfault 3 was needed to explain the widespread structural damage along the coast of Miyazaki Prefecture indicated in historical records. The total seismic moment as estimated from the fault model was 9.8 × 1020 Nm (Mw 7.9), assuming a rigidity of 3.5 × 1010 N/m2 (The Headquarters for Earthquake Research Promotion, 2022).

Computed tsunami heights along the coast in Miyazaki Prefecture from Nobeoka City to Kushima City were output with 3 arc-second grid intervals, as indicated in Fig. 5. The computed tsunami heights along the coast of the Miyazaki plain agree with tsunami heights of 4–5 m in historical records. Computed tsunami heights in the southern part of Miyazaki Prefecture were larger than those in the northern part. This is due to the location of the tsunami source and variations in topography near the coast with many lowland plains in the northern part of Miyazaki Prefecture and deeply indented coastlines developed in the southern part. Differences in topography between the present and the past may also have an effect.

Fig. 5
figure 5

Computed tsunami heights using fault model of 1662 Hyuga-nada earthquake along coast of Miyazaki Prefecture indicated as red line. Blue dots show tsunami heights in historical records from previous studies (Hatori, 1985; Tsuji et al., 2018) including inundation heights and runup heights. Pink dot shows elevation where event deposits were identified in this study. Yellow dots show elevations where event deposits were identified in previous studies (Ikuta et al., 2014; Kamataki et al., 2013, 2015; Niwa et al., 2020)

Full size image

At Komei in Nichinan City, the first wave of tsunami resulting in the maximum wave arrived 20 min after the earthquake occurred. The flow depths from the results of numerical simulation for Komei and the detailed surveyed points are indicated in Fig. 6. The computed tsunami inundation area at Komei covers locations where event deposits were found by field survey.

Fig. 6
figure 6

A Map of Nichinan coast in Miyazaki Prefecture with contour interval of 10 m. B Computed tsunami inundation area at Komei in Nichinan City using fault model of 1662 Hyuga-nada earthquake with contour interval of 10 m. Surveyed points are described; pink circles are identified event deposits and crosses are unidentified

Full size image

6 Discussion

The estimated fault model of 80 km length can explain the tsunami inundation area at Komei. A 2 m slip at subfault 3 can explain the widespread structural damage along the coast of Miyazaki Prefecture as seen in historical records. In the case of typical repeated M7 class earthquakes in Hyuga-nada, an average slip of approximately 1 m and a maximum seismic intensity of 5 were observed. Considering the seismic intensity of 6 estimated from historical records for the 1662 Hyuga-nada earthquake, a 2 m slip can be considered an appropriate value.

Computed tsunami heights represent heights at points on the sea surface along the coast. On the other hand, tsunami heights of historical records represent inundation heights or runup heights at points on the land. Thus, the compared tsunami heights are not at the same points; different results of computed and historical tsunami heights are indicated in some cases. At some points, it is assumed that tsunami heights in historical records are high because of the effects of the runup tsunami along the river.

Computed tsunami inundation areas reached to all points where event deposits were found by this field survey at Komei and by previous studies at Miyazaki plain. At all six points where event deposit was identified by field survey, distributed over an area of 5 m elevation, the thickness of event sand layers is almost 5–10 cm. Computed flow depths of 1–4 m at these points are a good fit in the tsunami inundation simulation. At the highest point where event deposit was not identified by this survey, located at 15 m elevation, computed tsunami was not reached. Therefore, the estimated fault model can explain the distribution of event deposits that were possibly formed by the 1662 tsunami.

In previous studies, event deposits interpreted by the 1662 tsunami were found at Shimayama in the Miyazaki plain (Kamataki et al., 2015). At this site, local subsidence due to the 1662 event was also identified by Niwa et al. (2020). Event deposits possibly from the 1662 tsunami were also identified at Kakiwara (Kamataki et al., 2013) and at Awakigahara (Ichihara et al., 2015) in Miyazaki plain. Computed tsunami heights along the coast of Miyazaki plain cover elevations where event deposits possibly from the 1662 tsunami were found in previous studies.

7 Conclusions

To determine the tsunami source of the 1662 Hyuga-nada earthquake, we surveyed geological evidence of the tsunami and estimated the fault model of the earthquake. The numerical simulation of the tsunami was carried out using the constructed fault model; the computed tsunami inundation area can explain the distribution of the possible tsunami event deposits at Komei in Nichinan City. The simulation also can explain tsunami deposits found in Miyazaki plain in previous studies. Computed tsunami heights can similarly explain the tsunami heights along the coast of Miyazaki Prefecture in historical records. From the fault model of the interplate earthquake, including the source area of slow earthquakes, the seismic moment of the 1662 Hyuga-nada earthquake was calculated to be 9.8 × 1020 Nm (Mw 7.9), assuming a rigidity of 3.5 × 1010 N/m2. It is possible that the 1662 Hyuga-nada earthquake is the largest in Hyuga-nada in historical records, caused by large slip at the active source area of slow earthquakes.