Research ArticleGEOPHYSICS

Mechanism of the 2015 volcanic tsunami earthquake near Torishima, Japan

See allHide authors and affiliations

Science Advances  04 Apr 2018:
Vol. 4, no. 4, eaao0219
DOI: 10.1126/sciadv.aao0219


Tsunami earthquakes are a group of enigmatic earthquakes generating disproportionally large tsunamis relative to seismic magnitude. These events occur most typically near deep-sea trenches. Tsunami earthquakes occurring approximately every 10 years near Torishima on the Izu-Bonin arc are another example. Seismic and tsunami waves from the 2015 event [Mw (moment magnitude) = 5.7] were recorded by an offshore seafloor array of 10 pressure gauges, ~100 km away from the epicenter. We made an array analysis of dispersive tsunamis to locate the tsunami source within the submarine Smith Caldera. The tsunami simulation from a large caldera-floor uplift of ~1.5 m with a small peripheral depression yielded waveforms remarkably similar to the observations. The estimated central uplift, 1.5 m, is ~20 times larger than that inferred from the seismologically determined non–double-couple source. Thus, the tsunami observation is not compatible with the published seismic source model taken at face value. However, given the indeterminacy of Mzx, Mzy, and M{tensile} of a shallow moment tensor source, it may be possible to find a source mechanism with efficient tsunami but inefficient seismic radiation that can satisfactorily explain both the tsunami and seismic observations, but this question remains unresolved.


Most typical tsunami earthquakes (1) are known to occur near deep-sea trenches (25), although their mechanisms are still not fully understood. Submarine landslides triggered by earthquakes can cause unexpected tsunamis (68), although these events may not be regarded as tsunami earthquakes. Another type of tsunami earthquake has been reported to occur near the volcanic island of Torishima along the Izu-Bonin arc. The occurrence of such earthquakes is quasi-regular, approximately every 10 years, always accompanying anomalously large tsunamis (913). The focal mechanism is always a type of non–double couple known as compensated linear vector dipole (CLVD) (14) with the vertical T axis. The moment magnitude (Mw) is always 5.6 to 5.7 (fig. S1). Most CLVD earthquakes with a vertical T (or P) axis have been reported at active volcanoes (1519) on land or on an island but not under the sea. The Torishima earthquakes are unique in that they are submarine CLVD events accompanying anomalously large tsunamis. Tsunamis from the 1984 Torishima earthquake were analyzed (12), but the available data were limited to tide gauge records inside bays more than 200 km from the epicenter; thus, it was challenging to retrieve quantitative information about the anomalous tsunami source relative to the seismic source.

We deployed a pressure gauge array at water depths of 1470 to 2240 m off east of Aogashima Island about 100 km from the epicenter toward ~N10°E from May 2014 to May 2015, which timely recorded seismic and tsunami waves from the most recent event (Mw = 5.7, 2 May 2015) (Fig. 1A). The array consists of 10 stations (table S1) forming equilateral triangles with a station interval of 10 km. Figure 1B shows a comparison of the array records between a volcanic earthquake and an ordinary shallow thrust earthquake (Mw = 5.9, 10 May 2015), which occurred near the Izu-Bonin trench ~200 km to the southeast of the array (Fig. 1A). Amplitudes of seismic waves are similar between the two earthquakes, yet tsunamis are visible only in records of the volcanic event, clearly demonstrating its tsunami earthquake nature.

Fig. 1 Observation outline.

(A) Bathymetry map showing the locations and mechanisms of two events, the 2 May 2015 volcanic tsunami earthquake (Mw = 5.7, depth = 10 km) and the 10 May 2015 ordinary thrust earthquake (Mw = 5.9, depth = 6 km). Ten small triangles indicate stations of the ocean bottom pressure gauge array. (B) Comparison of band-pass–filtered (1 to 80 mHz) records between the 2 May event and the 10 May event. The vertical axis is not scaled so that the amplitudes can be visually compared. Amplitudes of seismic waves are similar between the two events, yet tsunamis are visible only on the 2 May records. The arrows indicate the transition from a very weak downswing to the major upswing of tsunami motion.


To locate the tsunami source, we first performed a waveform analysis of dispersive tsunamis recorded at the array. The tsunami traces are narrowly band-pass–filtered around targeted frequencies (see fig. S2A for an example of 4 mHz). Assuming plane waves, the tsunami peak arrival times for each station are least squares–fitted by a straight line that gives the phase slowness vector and the reference phase arrival time (fig. S2B). The slowness vector is given as the absolute slowness value and wave propagation direction (table S2). The arrival time is reduced to the travel time assuming that the tsunami was generated instantaneously at the earthquake origin time and that the peak of a narrowly band-pass–filtered wave train came from the center of the tsunami source. The obtained absolute slowness value is consistent with the theoretically predicted value (fig. S2C), which depends on the frequency and water depth of the reference station (20, 21). The frequency-dependent travel time (2224) and propagation direction at the reference station are sensitive to the tsunami source location (Fig. 2; see also fig. S3). We place five hypothetical tsunami sources, one at the caldera center and four along the rim (Fig. 2A). The travel time is sensitive to the source location in the north-south direction (Fig. 2C), whereas the propagation direction is sensitive to the source location in the east-west direction (Fig. 2D). The tsunami source within the caldera explains the observations better than those along the caldera rim. Note that the epicenters located by the Japan Meteorological Agency (JMA), the U.S. Geological Survey (USGS), and the Global Centroid Moment Tensor (GCMT) are scattered more widely than these tested tsunami sources. We define the misfit of the ray theoretical calculation (25, 26) to the observation as a function of the assigned source location (Fig. 2B). The minimum of the misfit values lies well within the Smith Caldera, confirming tsunami generation within the caldera. In the above analysis, the phase advance due to the initial phase of cylindrically spreading waves (27, 28) has been assumed to be insignificant.

Fig. 2 Estimate of tsunami source location.

(A) Five tsunami source locations tested in this study (stars). One is within the caldera and four are along the rim. The CMT diagrams obtained by GCMT, USGS, and JMA are shown at their epicenters. (B) Spatial distribution of the misfit of the calculated phase travel time and propagating direction to the observed values. The minimum misfit is located within the caldera. The small circle with a cross indicates the estimated location of source periphery with its uncertainty. (C) Test result showing the sensitivity of tsunami travel time to the shift of source location from north to south through C (A). Black line with error bars shows the observed travel times. (D) Test result showing the sensitivity of tsunami propagating direction to the shift of source location from east to west through C (A). Black line with error bars shows the observed propagating directions.

The initial polarity of the tsunami wave recorded at the array is down, although its amplitude is very small (Fig. 1B). To confirm this, the observed 10 waveforms are mutually displaced along the time axis so that their major upswings start at a common zero-cross time. The average waveform after this mutual time shift clearly shows a weak downswing before the major upswing (Fig. 3). Such a precursory first motion with reversed polarity is often observed for trans-oceanic tsunamis due to dispersion of tsunamis gravitationally and elastically coupled with Earth at long periods (29). However, the tsunami travel distance to the array is too short to cause such dispersion; hence, this small precursory downswing must be due to the tsunami source. It is most simply explained to be caused by a peripheral subsidence associated with the central uplift. The trace from small downswing to major upswing can be regarded as the signal from the boundary between the central uplift and the peripheral small subsidence. We also determine the slowness vector and travel time of the zero-cross onset by the least-squares method. The result (table S2) indicates that the zero-cross signal propagates through the array with the long-wave speed (30), in contrast to the frequency-dependent speed of band-passed peaks. We accordingly back-project the zero-cross signal to the earthquake origin time along the long-wave ray path and locate the source at the boundary of the central uplift near the northern caldera rim (Fig. 2B). This source implies that the tsunami was mainly generated by the uplifting of the entire caldera floor.

Fig. 3 Average trace of 10 tsunami records mutually time-shifted to a common zero-cross time.

The dotted lines indicate the SD of these time-shifted records. This figure emphasizes the weak downswing before the major upswing of the tsunami motion.

To obtain a model of sea-surface uplift that can explain the observed tsunami waveforms, we perform a simulation for the dispersive tsunamis radiating away from the source by solving linear Boussinesq equations using the JAGURS code (31). On the basis of the source location estimated from the ray-tracing analysis (Fig. 2), we assume the main uplift inside the caldera and a small subsiding ring outside. This sea-surface displacement is assumed to occur instantaneously. The tsunami phase arrivals (Fig. 2C), incident azimuths (Fig. 2D), and waveforms (Fig. 4) are all explained well by a sudden uplift of the caldera floor, indicating that the uplift took place within a time much shorter than the escape time (~40 s) of the tsunami peak away from the caldera rim. This is consistent with source times of less than 10 s estimated from seismic waves (see Materials and Methods). These observations seem to suggest that the tsunami and seismic events occurred simultaneously within a time scale of 10 s or less. We assume that the central major uplift and peripheral minor depression of the tsunami source can be approximated by a superposition of two Gaussian curves (fig. S4). We performed a grid search to determine the best combination of amplitude and radius of the main uplift to match the observed waveforms. This search leads to a sea-surface uplift of ~1.5 m with a source size of 4.1 km comparable to the caldera size (Fig. 4, B and C). Agreement between the simulated waveforms and the Ocean Bottom Pressure gauge (OBP) records is remarkable (Fig. 4A). Our model also reproduces a large tsunami wave with an amplitude of around 60 cm observed at a tide gauge inside a bay of the Yaene Port on Hachijo Island, 180 km northward from the epicenter (Fig. 4D and fig. S5). Our question now is whether this tsunami-based estimate of sea-surface uplift is compatible with the seafloor uplift expected from the seismologically determined CLVD mechanism.

Fig. 4 Best initial sea-surface displacement model.

(A) Observed (black) and simulated (red) waveforms at OBP stations (see Fig. 1A for locations). The blue horizontal line represents the time window to calculate the normalized root mean square (NRMS) misfit at each station. (B) The sea-surface displacement for a model of A = 1.5 m and R = 4.1 km (see Materials and Methods). Contour lines represent the bathymetry near the Smith Caldera. (C) The cross-sectional profile of the initial sea-surface displacement. (D) Observed (black) and synthetic (red) waveforms at the Yaene Port on Hachijo Island (see fig. S5).

Although many physical mechanisms have been proposed for the vertical CLVD of volcanic earthquakes (18), little consensus has been reached thus far. In view of this, we consider a very phenomenological model that is simple in geometry and efficient for tsunami generation. It is a horizontal square prism undergoing lateral contraction and vertical expansion without volume change (Fig. 5A). This volume source generates relatively uniform seafloor displacement and hence relatively uniform sea-surface displacement (20, 32). We assume the prism to have a thickness (T) of 2 km and an area (L2) of 8 km × 8 km, which is comparable to the caldera size in a uniform half-space. The deviatoric or scalar moment (18) is set at 4 × 1017 Nm (see Materials and Methods), which corresponds to a relative vertical displacement of ~14 cm across T and a relative horizontal displacement of ~28 cm across L. If the prism is placed at a depth of 9.2 km by referring to the USGS depth of 10 km, the resultant vertical surface displacement is only 4 cm at maximum (Fig. 5A), in marked contrast to the tsunami-based estimate of ~1.5 m. The uplift is also broadened well beyond the caldera rim in contrast to the tsunami-based estimate of source size. The peripheral depression estimated from tsunami observation is not reproduced. The CLVD source at a seismologically estimated depth cannot explain the observed large tsunamis. Taking the uncertainty in hypocentral depth determination into account, we shift the source depth upward without changing other parameters (Fig. 5A). A shallower source produces greater seafloor uplift over a narrower area with a more pronounced peripheral depression. If the source is placed at the shallowest depth of 1.2 km, the maximum uplift is 8 cm, which is still 1/20 of the tsunami-based estimate. Thus, even a shallower source than the seismologically estimated depth cannot explain the observed tsunami amplitudes. We have to involve a source that is efficient in generating tsunami waves but inefficient in generating seismic waves.

Fig. 5 Vertical ground deformations due to two phenomenological sources.

(A) Vertical T–CLVD volume source. Top: source geometry. The dimensions are 8 km × 8 km (horizontal) and 2 km (vertical) with a top depth at δ. Arrows show the relative displacement direction. Bottom: Calculated vertical surface deformation profile for different top depths δ. The source area is shaded. The deviatoric moment is fixed at 4 × 1017 Nm so that the relative vertical displacement is 14 cm. (B) Vertical opening of a horizontal tensile fault. Top: source geometry. The horizontal fault surface is 8 km × 8 km at a depth of 0.2 km. Arrows show the dislocation direction. Bottom: calculated vertical seafloor deformation profile (red curve). The sharp edges of this profile should be truncated when the seafloor disturbance is transmitted to the sea surface (20, 32). The dislocation is fixed at 1 m. The deformation profile in an infinite medium is shown for comparison (blue curve). Note that the vertical scale is different between (A) and (B).

Among typical elementary sources, including shear faulting, tensile faulting, and isotropic volume change (27), only the vertical opening of a horizontal tensile crack at a shallow depth is intrinsically tsunamigenic but aseismic (see Materials and Methods for details). Figure 5B shows the static vertical displacement of such a fault with an area of 8 km × 8 km at a depth of 0.2 km in a uniform half-space. Although this shallow depth is chosen as an extreme case, the qualitative behavior remains unchanged as long as the depth is much smaller than the fault size. The vertical dislocation across the fault is set at 1 m. The calculated surface displacement profile almost mimics the given dislocation profile at the source, implying that the displacement field is mostly confined right above the horizontal tensile fault so that seismic waves would be radiated mostly upward from the fault surface. Figure 5B also shows the displacement profile 0.2 km above the fault in an infinite medium. Vertical displacement in this case is symmetric across the fault surface so that seismic waves would be radiated equally from the upper and lower surfaces of the fault. Confinement of displacement on the upper side of the fault is a unique feature of a vertically opened horizontal tensile fault at a shallow depth, which we suggest to be the predominant process of the Torishima tsunami earthquake. The seafloor uplift of ~1.5 m estimated from the tsunami observation is largely due to this process. The isotropic and deviatoric moments (18) of the CLVD source (Fig. 5A) are MICLVD = 0 and MDCLVD = 4, respectively, whereas those of the tensile fault (Fig. 5B, but with a dislocation of 1.5 m) are MITENSILE ≈ 48 and MDTENSILE ≈ 28 in units of 1017 Nm (see Materials and Methods). The seismologically undetected component (tensile fault) was far greater than the seismologically detected component (CLVD).


Although the vertical opening of a tensile fault may be commonly associated with vertical T–CLVD volcanic earthquakes, it would be challenging to detect it by seismological means unless observations were made directly above the source. For example, the Bárdarbunga, Iceland, volcanic earthquake (Mw = 5.6, 29 September 1996) was identified as a vertical T–CLVD event (16) by far-field data and confirmed by using a regional array at epicentral distances of ~100 km (33); the array was still too far from the source to prove the presence or absence of the vertical opening of a tensile fault. The process should be better detected by geodetic observation of surface deformation right above the source. The Sierra Negra, Galápagos (34, 35), volcanic earthquake (Mw = 5.5, 22 October 2005) was a vertical T–CLVD event (19), for which a progressively increasing caldera uplift was observed by the Global Positioning System network. However, this network failed 4.5 hours before the earthquake (36).

We have argued that the Torishima tsunami was generated largely by the vertical opening of a horizontal tensile fault right beneath the Smith Caldera. Long-period seismic waves were poorly excited by this mechanism. The seismologically determined mechanism is vertical T–CLVD, a solution obtained with the assumption of the vanishing isotropic component of a moment tensor (18, 19). We have shown that this mechanism contributed little to the observed tsunami excitation. The mechanism that satisfies both seismic and tsunami data simultaneously remains to be explored. Because this type of event recurs in the Smith Caldera at approximately every 10 years (fig. S1), it may be understood as a thermomechanical response of the sub-caldera layer to the repeatedly occurring magma ascent. Because shallow horizontal tensile faults are inefficient in seismic radiation, it is difficult to detect and study shallow tensile faults with seismological data. In contrast, tsunamis are efficiently excited by the vertical motion of the top side of a tensile fault, and analyses of tsunami data will help resolve details of horizontal tensile faults. This in turn will lead to a better understanding of the physical processes of magmatic injection.


Event and array information

The reported origin time of the 2 May 2015 event is 16:50:42.6 UT (JMA, or 16:50:43.0 (USGS, The centroid time is 16:50:49.7 (GCMT, The reported hypocenter is 31.517°N, 140.355°E, 28 km (JMA) or 31.529°N, 140.213°E, 10 km (USGS). The centroid location is 31.47°N, 139.94°E, 12 km (GCMT). The reported Mw is 5.7 (JMA and GCMT). The reported moment tensor components Mrr, Mtt, Mff, Mrt, Mrf, and Mtf (in units of 1017 Nm) are 4.47, −2.43, −2.04, 0.44, 0.14, and −1.08 (JMA), or 4.44, −2.73, −1.72, −1.06, −2.81, and −0.826 (GCMT), respectively, where the isotropic moment MI = (Mrr + Mtt + Mff)/3 (18, 37) is assumed to be zero. Their principal axes T, N, and P are 4.50, −1.10, and −3.40 (JMA), or 4.99,−1.43, and −3.56 (GCMT), respectively, with which the deviatoric moment MD (18) is (4.50 + 3.40)/2 = 3.90 or (4.99 + 3.56)/2 = 4.45, respectively. This definition of deviatoric moment is different from that by Bowers and Hudson (37). The plunge of the T axis is 86° (JMA) or 79° (GCMT), either of which is greater than 60°. The non–double-couple component parameter ε is 1.10/4.50 = 0.244 (JMA) or 1.43/4.99 = 0.286 (GCMT), either of which is greater than 0.2. This event is classified as a vertical T–CLVD earthquake (18, 19). In the present paper, we idealized this mechanism as a pure vertical T–CLVD with moment tensor components of 16/3, −8/3, −8/3, 0, 0, and 0 and a scalar moment of 4 × 1017 Nm.

The array consisted of 10 stations, each deploying a free-fall/pop-up ocean bottom pressure meter equipped with an absolute pressure gauge (APG). The APG (PARO-8B7000-I-005) used two quartz crystal resonators (one to convert analog force inputs generated by a Bourdon tube with digital outputs and the other to measure internal temperature to correct for its effect on resonant frequency), and it can measure absolute water pressure up to a depth of 7000 m (68.95 MPa) with nano-resolution technology (38). The sampling and cutoff frequencies were set at 4 and 0.7 Hz, respectively. The array was deployed east of Aogashima Island in a 1-year period from May 2014 to May 2015. Deployment and recovery were carried out by research cruises YK14-07 and YK15-08 of R/V Yokosuka (JAMSTEC), respectively. The coordinates of stations A01 to A10 are given in table S1. We refer to station A05 as the reference station, at which the water depth is 1762 m, corresponding to a long-wave speed of 131 m/s. The array formed equilateral triangles with a site interval of 10 km and a maximum length of 30 km (Fig. 1A). The length (30 km) corresponds to the wavelength of the long wave with a period of 230 s passing through the reference station.

Phase analysis

We decimated the 1-day records of 2 May 2015 to a sampling rate of 1 Hz and applied a two-pole Butterworth band-pass filter with corners at 0.0005 and 0.02 Hz forward and backward over the data. The records were then time-windowed between 60,000 and 63,000 s counted from the beginning of the day. Then, the record sections were narrowly band-pass–filtered forward and backward with corners at n − 1 and n + 1 mHz around a central frequency of n millihertz, where n = 2, 3, …, 10. Figure S2A shows the waveforms in the case where n = 4. We read the peak arrival time on each of the narrowly band-pass–filtered wave trains. At higher frequencies, a picked peak may show an inconsistent arrival relative to those at the adjacent stations or those at the lower frequencies. We picked the one-cycle earlier arrival in such cases.

The data set of the phase arrival time Ti at the ith station (i = 1, 2, …, 10) for each of frequency n was used to determine the frequency-dependent phase slowness vector of tsunami propagation through the array under the plane wave approximation. We deployed a local rectangular coordinate system with its origin at the fifth station (reference station) and the x and y axes directing east and north, respectively (fig. S2B and table S1). For the phase propagating with a slowness vector (Sx, Sy) and reaching the origin (x = y = 0) at an arrival time To, we obtained the following 10 observation equations.Embedded Image(1)where (Xi, Yi) are the coordinates of the ith station. These equations were solved for the three unknowns Sx, Sy, and To by the method of least squares, which can be written asEmbedded Image(2)

Once Sx, Sy, and To were obtained, the error associated with measurement Ti was estimated asEmbedded Image(3)

The errors associated with the estimates of Sx, Sy, and To were calculated asEmbedded Image(4)

The phase slowness S and phase propagation direction θ can be calculated from the slowness vector (Sx, Sy) asEmbedded Image(5)

The errors associated with these quantities were estimated asEmbedded Image(6)

Ray tracing

Phase velocity C of a tsunami wave depends on the wave number k and water depth H asEmbedded Image(7)where f is frequency and g is gravitational acceleration (20, 21). Using this dispersion relation, the bathymetry map was converted into frequency-dependent phase velocity maps. On the basis of the phase velocity maps, frequency-dependent ray paths and travel times were calculated by numerically integrating the ray equations of tsunamis, which were identical to those of seismic surface waves (fig. S3) (25, 26).

A grid search was performed for the tsunami source in such a way that the ray theoretical phase arrival time Tray and propagating direction θray best agreed with the observed phase arrival time To ± dTo (Fig. 2C) and propagating direction θo ± dθo (Fig. 2D), where dTo and dθo were taken to be equal to σT and δθ, respectively. We defined the misfit function at each grid point asEmbedded Image(8)where n indicates the frequency in millihertz at which measurement is made. Grid search was performed at an interval of 0.005°. The tsunami was assumed to be generated at the earthquake origin time. The spatial distribution of the M value is shown in Fig. 2B. The point minimizing M was located within the Smith Caldera (28).

Tsunami source modeling

As described in the Results, the tsunami source consisted of a major central uplift and a minor peripheral depression with their boundary in the vicinity of the caldera rim (Fig. 2B). Assuming a vertical T–CLVD–type mechanism of the earthquake, we modeled an axis-symmetrical initial sea-surface displacement source that is expressed byEmbedded Image(9)where r is the distance (in kilometers) from the center of the caldera. This model has a main central uplift of A (in meters) surrounded by a minor depression of −A/10 (in meters), and the radius of the uplift is R (in kilometers) (fig. S4). We varied A from 0.5 to 5.5 m at 0.1-m intervals and varied R from 0.5 to 5.5 km at 0.1-km intervals. For each initial sea-surface displacement model, the linear Boussinesq equations were solved with the JAGURS code (31) to compute tsunami waveforms at the OBP stations. The complex bathymetry near the Smith Caldera was taken into consideration by setting a grid size of 10 arc sec (250 m) near the caldera and 30 arc sec (750 m) in other regions.

We searched for the best parameters A and R that minimize the misfit between observed and synthetic waveforms. We defined the NRMS misfit for the kth station asEmbedded Image(10)where i refers to the sampling time step; N is the sampling number at each station; obsi and simi are observed and synthetic waveforms, respectively; and Embedded Image is the average of the observed record over the time window at each station. The final NRMS misfit value is an average of the misfits over the stationsEmbedded Image(11)where M is the number of stations (10 here). We set the time window to include the precursory downswing signal and the first and second major upswing and downswing pulses (Fig. 4A). The calculated NRMS misfit value was plotted as a function of A and R in fig. S4. The misfit was small for the amplitude A range of 1.3 to 1.8 m and the radius R range of 3.7 to 4.3 km. The minimum misfit occurred for A = 1.5 m and R = 4.1 km. The waveforms simulated with these values were in remarkable agreement with the observations (Fig. 4A). We concluded that a large uplift of ~1.5 m was confined within the caldera, and a minor depression surrounded the main uplift at least on the northeastern part outside the caldera rim (Fig. 4, B and C).

The JMA issued a tsunami advisory when a tsunami height of 60 cm was recorded at a tide gauge inside a bay of the Yaene Port at Hachijo Island, ~180 km northward from the epicenter (fig. S5A). This height is in marked contrast to the tsunami height of ~2 cm observed at the OBP stations ~100 km toward north-northeast from the epicenter. If our tsunami source model derived from the open-sea observation is a reasonable one, it should explain the above observation as well. To calculate the synthetic tide gauge record, we prepared the nested bathymetry grid data to accurately represent the bay shape (fig. S5), with the finest grid near the port at 0.2 arc sec, or 2 m. Because the nonlinear effects on tsunami waves were not negligible near the coast, we included the nonlinear terms to compute waveforms near the port (bathy04 and bathy05 in fig. S5). The simulated waveform well reproduced the observed waveform at the Yaene Port (Fig. 4D), despite the fact that this observed record was not used in source modeling.

Phase advance and delay

We ignored the effects of the initial phase of cylindrically spreading waves and finite source duration on the estimate of source location in the phase analysis (Fig. 2). We discuss here these effects. In the tsunami simulation for our source model (Fig. 4A), it took 43 s for the peak of disturbance to propagate from the source to the edge at a radius of 4.0 km. On the contrary, the long wave with a speed of 89 m/s traveled only 3.8 km within this time duration, assuming a water depth of 800 m inside the Smith Caldera. The higher speed of the peak propagation in the tsunami simulation was an expression of the initial phase effect of cylindrically spreading waves (27, 28). If this effect is taken into account, the source location as estimated in Fig. 2B is shifted southward by ~0.2 km. However, the location after the shift is still inside the caldera.

We assumed an instantaneous seafloor deformation at the earthquake origin time in the phase analysis. According to the GCMT catalog, the half-duration of the 2015 event was 1.8 s, whereas the centroid time shift was 6.7 s. The source duration expected from an empirical scaling with scalar moment for vertical CLVD earthquakes (18) was 7 to 8 s. Because we located the tsunami source by back-projecting the tsunami phase up to the earthquake origin time, we made a correction for the source location using the reported centroid time shift. This correction shifted the source location northward by ~0.6 km. Thus, the phase shift due to the initial phase of cylindrically spreading waves and the phase delay due to a finite source duration acted in opposite directions so that their effects on source location may be canceled out.

Zero-cross time

We decimated the 1-day records of 2 May 2015 to a sampling rate of 0.25 Hz. No filter was applied except for the high-cut filter associated with the decimation. The records were then time-windowed between 60,000 and 63,000 s and detrended. In this time window, we searched for the maximum peak and its preceding zero-cross time. We confirmed for each of the 10 records that the identified zero-cross time represented transition from the weak downswing to the major upswing of tsunami motion. This weak downswing signal was observed regardless of the choice of causal or acausal filters. The zero-cross times were read, and a slowness analysis was made for them by the same method as for the phase arrivals. The result is given in table S2. We displaced the 10 sets of tsunami traces along the time axis so that their major upswings had a common preceding zero-cross time. We then calculated the average and SD of these time-shifted traces to illuminate the downswing stage just before the zero-crossing (Fig. 3).

CLVD model

Our vertical T–CLVD source model was very phenomenological. We assumed such internal deformation that a finite volume underwent lateral contraction and vertical expansion without volume change. We took a rectangular coordinate system (x,y,z) where the z axis directed vertically upward with z = 0 at the surface of the half-space. The CLVD source was a square prism extending vertically from z = −δ − T to −δ, laterally from x = −L/2 to +L/2 and from y = −L/2 to +L/2. The z-parallel relative displacement was DV between z = −δ − T and −δ. The x-parallel (y-parallel) relative displacement was −DH between x = −L/2 and +L/2 (y = −L/2 and +L/2). The CLVD source required that DVL2 = 2DHLT. Its displacement field was calculated by distributing a large number (1024) of z-normal tensile fault planes between z = −δ − T and −δ and a large number (1024) of x-normal (y-normal) tensile fault planes between x = −L/2 and +L/2 (y = −L/2 and +L/2) using the existing codes (39, 40). The centroid moment tensors of the above three decomposed sources wereEmbedded Image(12)respectively, where μ and λ are Lame’s constants. Because DVL2 = 2DHLT, the resultant centroid moment tensor is such CLVD thatEmbedded Image(13)with a scalar moment of (3/2) μDVL2. If we equate this scalar moment to the reported scalar moment of 4 × 1017 Nm by using a μ value at the reported centroid depth, we can recover the potency (41), (3/2) DVL2, and hence the relative vertical displacement, DV. The potency recovered should be less sensitive to the local heterogeneity of rigidity than the seismic moment (42, 43), if the centroid depth is a correct estimate of the source depth. The thickness T and length L were fixed at 2 and 8 km, respectively. The top depth δ was varied. We obtained DV to be 14 cm, using a μ value of 30 GPa corresponding to the shallowest layer of Preliminary Reference Earth Model (PREM) (44). Relative lateral displacement −DH is given by the relation DVL2 = 2DHLT so that DH = 28 cm. Figure 5A shows the vertical surface displacement profile due to a CLVD source in the three cases where δ = 0.2, 4.2, and 8.2 km. The central surface uplift was 4 to 8 cm, whereas the peripheral depression was 0 to 1 cm.

Seismic and tsunami waves from a horizontal tensile fault

The vertical opening of a horizontal tensile crack at a shallowest depth has such a unique feature that it efficiently generates tsunamis but poorly excites seismic waves. This remarkable feature can be explained in a more formal way by using a modal expression of the displacement field due to a horizontal tensile fault whose moment tensor is given byEmbedded Image(14)where UZ is vertical dislocation across the horizontal square fault plane with length L at depth δ. The displacement at point x for a general moment tensor acting stepwise in time at point xs is given byEmbedded Image(15)(27), where * indicates the complex conjugate and ωi, Qi, ui, and eipq denote the eigen angular frequency, quality factor, displacement function, and strain function of the ith normal mode, respectively. The corresponding discrete horizontal wave number is denoted by k. For the horizontal tensile fault, this expression reduces toEmbedded Image(16)where Δi denotes the ith mode dilatation function and Embedded Image represents the ith mode normal stress acting on the (x,y) plane at point xs. The modes included predominantly gravitational spheroidal modes, that is, tsunami modes (4547). The Embedded Image eigen functions of tsunami modes behaved very differently from those of long-period Rayleigh modes at the same period. They rapidly increased from zero at z = 0 to take their maxima at the bottom of the ocean and slowly decreased with increasing depth (48), whereas those of long-period Rayleigh modes remained nearly zero at the bottom relative to their maxima at much greater depths. For example, the Embedded Image value of 200-s Rayleigh wave at the ocean bottom at 3-km depth was only 0.66% of the maximum near 300-km depth according to the PREM model (44). Thus, if the vertical opening of a horizontal tensile fault happened at a shallowest depth below the seafloor, it generated tsunami waves efficiently but long-period Rayleigh waves poorly.

Displacement due to a finite horizontal tensile fault

We next examined the near-field static displacement due to a finite tensile source in a uniform half-space. The fault plane at depth δ (0.2 km) extended horizontally from x = −L/2 to +L/2 and from y = −L/2 to +L/2, where L = 8 km, across which a z-parallel dislocation was UZ. The source was very shallow in the sense that δ/L << 1. The vertical displacement field was calculated using the existing code (39, 40). Figure 5B shows the vertical surface displacement profile in the case of UZ = 1 m. The vertical ground displacement right above the source was almost equal to the given dislocation UZ and was almost zero further outward, implying that displacement was concentrated on the top side of the fault and was poorly partitioned into the bottom side. This displacement field was in marked contrast to that for a very deep fault (δ/L >> 1), where the vertical displacement was equally partitioned into the top and bottom sides of the fault. Such a contrast points to a great efficiency of tsunami generation from the shallowest tensile fault.

The moment tensor for this fault can be decomposed into a combination of the vertical T–CLVD source and the isotropic source.Embedded Image(17)

This decomposition implies that seismic radiation from the vertical T–CLVD and that from isotropic contraction are practically indistinguishable from each other if the source is shallow enough. Such indistinguishability explains the well-known trade-off between the vertical CLVD and isotropic components of a moment tensor at a shallowest depth (18, 49).


Supplementary material for this article is available at

fig. S1. Repeatedly occurring (approximately every 10 years) tsunami earthquakes near the Smith Caldera.

fig. S2. Phase slowness analysis.

fig. S3. Frequency-dependent ray trajectories and wave fronts.

fig. S4. Search for the best initial sea-surface displacement.

fig. S5. Waveform simulation at the Yaene Port on Hachijo Island.

table S1. Array information.

table S2. Result of array analysis of phase and onset arrivals.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank T. Oi for his help with the instrumentation. We also thank the shipboard scientists and technical staff of cruise YK14-07 and cruise YK15-08 of R/V Yokosuka (JAMSTEC). For the tsunami simulation, we used the bathymetry data, JTOPO30, provided by the Marine Information Research Center of the Japan Hydrographic Association. The bathymetry data near shore were provided by the Japan Nautical Chart Web Shop. The topography data near Yaene Port were provided by the Geospatial Information Authority of Japan. The tide gauge data at Yaene Port were provided by JMA. We also thank A. Gusman for help in the processing of the bathymetry data. Funding: This work was supported by Japan Society for the Promotion of Science (JSPS) Grants-in-Aid for Scientific Research (KAKENHI) grant numbers 25247074, JP17J02919, JP17K05646, and JP16H01838. Author contributions: Y.F. organized the research project. H. Sugioka, A.I., and H. Shiobara contributed to the observational part of the project. O.S., S.W., and K.S. contributed to the analysis part. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The supporting data are in part included in tables S1 and S2. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from Y.F. (fukao{at}

Stay Connected to Science Advances

Navigate This Article