Excitation of back-arc tsunamis from megathrust ruptures: The underdog hazard in the Sea of Japan

The 2011 Tohoku earthquake created a moderate tsunami in the back-arc Sea of Japan basin. This tsunami went largely unnoticed due to its small size and the significant coverage of the large fore-arc waves. We present a physical dislocation model for the excitation of back-arc tsunamis and identify fault dip as the main geometrical contributor to the propagation of back-arc tsunamis. Using numerical simulations and data from the 2011 event, we show that a combination of near- to intermediate-field horizontal and vertical dislocations as well as transient surface waves is necessary to reconstruct the back-arc propagation. We then simulate potential future earthquakes in the Japan trench and Nankai trough to investigate the back-arc tsunami hazard in the Sea of Japan. Our results show that the coseismic excitation of back-arc tsunamis can result in considerable waves exceeding 1 m from megathrust earthquakes.


Introduction
The 2011 Mw 9.0 Tohoku earthquake ruptured a ~200 km × ~400 km patch of the 75 Pacific-Okhotsk plate boundary along the Japan trench [Ammon et al, 2011;Ide et al, 2011;Shao et al, 2011;Lay, 2018] with unusually large slip, up to 50-80 m [Fujiwara, 2011;Ito et al, 2011a]. The shallow hypocenter at ~25 km [Dziewonski et al, 1981;Ekström et al, 2012] combined with the large slip resulted in ocean floor uplifts of up to ~5 m near the coastline, as well as ~30 m vertical and ~60 m horizontal displacements at the trench, respectively [Kido et al, 80 2011, Ito et al, 2011b. The large areas of uplift and subsidence at the ocean floor [Grapenthin & Freymueller, 2011] generated a devastating tsunami reaching 40 m along the Sanriku coast [Mori et al, 2011;Shimozono et al, 2012].
However, the coseismic deformation was not limited to the trench side of the Japan island -the fore-arc. In fact, geodetic measurements documented decimetric and centimetric 85 deformations in the western Japanese coastlines [Pollitz et al, 2011]. This so-called back-arc field of deformation would continue into the Sea of Japan and result in a complex pattern of uplift and subsidence at the sea floor. Numerical models also successfully predict centimetric deformation westward across the Sea of Japan reaching Russia [Pollitz et al, 2011;Sato et al, 2011]. Such back-arc deformations are expected from Okada algorithms for a buried, inclined 90 fault plane (Fig. 1).
as recorded by tide gauges around the Sea of Japan. Small waves were also reported by local fishermen at Niigata Prefecture [Murotani, personal. comm.]. Murotani et al [2015] showed that the arrival times of these tsunami waves matched the origin time of earthquake centroid and hence were excited coseismically. They modeled this tsunami using both vertical and horizontal components from a point source using Okada's [1985] algorithm. 100 Here we further investigate the characteristics of back-arc tsunamis in the Sea of Japan by considering the effects of surface waves from the mainshock in addition to static dislocations.
We also evaluate the back-arc tsunami hazard from various megathrust earthquakes in the Japan trench and Nankai Trough. We show that considerable (~1 m) back-arc coastal amplitudes in the Sea of Japan are possible from such events. Finally, we provide a simple physical model for the 105 generation of back-arc tsunamis highlighting the most important source parameters involved.

Tide gauge data 110
To investigate the propagation of the 2011 back-arc tsunami, we use data from three tide gauges along the western coastlines in Rudnaya Pristan, Preobrazheniye (both in Russia), Busan (Korea), as well as five tide gauges from Geospatial Information Authority of Japan (GSI) on the Japanese coast (Fig. 2a). These stations provide a reasonable coverage of the tsunami around the Sea of Japan. The Busan station is used to analyze the critical interaction of fore-and back-arc 115 arrivals. The time series from Russian and Japanese sites are sampled at 1-minute and 30-second intervals, respectively.

Tsunami Simulation
We apply the algorithm by Mansinha & Smylie [1971] to calculate vertical and 120 horizontal deformations as initial conditions for tsunami simulations. We use the MOST algorithm [Titov et al, 2016] to simulate the tsunamis and truncate the simulations at a shallow depth (10 m) close to the coastline. MOST has been extensively validated and applied in tsunami studies [Synolakis et al, 2008]. Our simulations are carried out using the ETOPO2 bathymetry grid with a resolution of 2 arc-minutes, and over 24-hr time windows with 5 s time steps to 125 satisfy the CFL stability conditions.

Contribution of surface waves
Surface waves from earthquakes can cause tsunamis and seiches for various basin types due to their high amplitudes and wide spectral range [Barberopoulou, 2004;Saito et al, 2019]. 130 Our simulations consider the initial deformation from Rayleigh waves which typically possess higher frequency content than the longer period tsunami waves. Using a normal modes approach, we calculate surface wave eigenfunctions down to ~150 seconds (see supplementary material).
This corresponds to a spatial resolution of ~30 km in the Sea of Japan, which is shallower than 4000 m, and should be sufficient considering the large source dimensions. 135

Simulation Results
We carry out tsunami simulations using (a) vertical deformation, (b) both vertical and horizontal deformations for the CMT solution, and (c) a combination of surface deformation and Rayleigh waves (more details are shown in supplementary material). Fig. 2a shows the 140 simulation result from the superposition of all three deformation sources as initial conditions. Each panel in Fig. 2(b-d) shows the simulated time series compared to recorded data (black curves) as well as their correlation coefficients for 1-hr time intervals (results for a longer simulation window are shown in Figs. S1 & S2). Similar to this study, Murotani et al [2015] were more successful in matching the 145 observed signals from stations at the far side of the Sea of Japan (see supplementary material).
The discrepancy is possibly due to the effect of bathymetry caused by horizontal component. We attribute the practically non-existing coseismic amplitudes and the signal complexity in Busan [similar issues can be seen at Okushiri in the north (Fig. 2)] to its proximity to deformation node

Potential Back-arc Tsunami Hazard in the Sea of Japan
The tectonic convergence at the Japan Trench and the Nankai Trough has resulted in multiple documented megathrust events in the east and south of Japan. Subduction patches capable of hosting similar megathrust events have been studied in detail along both margins [Ando, 1975;Furumura et al, 2011;Satake, 2015;Kusumoto, 2018]. We use the proposed 165 rupture locations and geometries from previous studies as potential future sources that are large enough to produce back-arc tsunamis in the Sea of Japan. Here, we do not include surface waves in the simulations as we are mainly investigating the largest amplitudes, not the first, high frequency portion of tsunamis.

Japan Trench
The subduction zone along the Japan trench can be considered as a set of potential rupture patches in an interwoven pattern [Satake, 2015]. A successive cascade of ruptures on these fault patches can occur during a single large megathrust, similar to the 2011 Tohoku event [Stein & Okal, 2011;Noda & Lapusta, 2013]. As a given paleoseismic window may not capture 175 long-term trends in seismicity [Salditch et al, 2020], such large rupture areas can be keys to understanding the maximum possible earthquake and tsunami hazards.
We construct a composite rupture area by joining the previous MS 7. 1-7.8 events in 1936, 1937, 1938 and 1978 in the Miyagi-oki sequence, spanning a total geographic area of ~300 km × ~150 km (Abe, 1977;Umino, 2006;Satake, 2015), i.e., roughly 50% of the 2011 mainshock 180 rupture area. We replace the 1936 event with its collocated 2005 counterpart because of better data (Table S2). Then using earthquake scaling laws [Geller, 1976] for a constant stress drop and assuming a uniform convergence rate along the trench, we calculate an average slip and magnitude (Mw=8.8) for this earthquake, scaled down from the 2011 Tohoku as reference (the Tohoku event has a Kagan aperture [Kagan, 1991] of ~23º from the total mechanism obtained 185 from the events in the sequence). We also make another composite source with similar geometry,

Parameters Controlling the Back-arc Tsunami
Here we investigate the effects of bathymetry, surface waves and earthquake source 205 parameters in the excitation and propagation of back-arc tsunamis.

Bathymetry & Surface Waves
As discussed in section 4, bathymetry has a significant role in the higher frequency content of back-arc tsunamis. This effect is more noticeable in the presence of smaller deformation, e.g., in the case of Busan where the 2011 tsunami arrived more than one hour later 210 than the Russian sites (Figs. 2b-d). Ray-tracing experiments (supplementary material) confirm that the Japanese sites as well as Busan experience considerably more complex wavefronts than those on the far side of the Sea of Japan, due to the entrapment of tsunami energetics in the complex southern bathymetry (Fig. S7).

Our results in section 2.3 show that the contribution of Rayleigh waves to the 2011 215
Tohoku back-arc tsunami at near-to intermediate-field locations may be important, as their addition slightly improves fits to the observed signal (~30% increase in CC) at higher frequencies (during the first few hours of the records after origin time) in Rudnaya Pristan but not in Preobrazheniye (see between 6-7 hrs in Fig. 2b). While these stations are approximately in the same radiation lobe of Rayleigh waves, the inconsistency may be due to rupture directivity. 220

Earthquake source
The extent to which each earthquake source parameter affects the back-arc tsunami is yet to be discussed. Fig. 1 is a simple cross-section of what we consider as the back-arc deformation model. This model is based on the general Okada solution for a centroid near the trench in the fore-arc. Such models are routinely used in tsunami simulations due to their 225 simplicity and robustness [Titov & Synolakis, 1998;Okal & Synolakis, 2004]. In Fig. 1, the earthquake epicenter and mechanism shown by a beachball is located on a fault plane depicted as an inclined, solid black line with a dip angle of , turning into a dashed line to imply the uncertainty in the extent of down-dip rupture propagation.
In the absence of bathymetry, the deformed ocean floor in Fig. 1 (shown in orange) 230 displaces a body of water from static coseismic deformation in the back-arc (shown in red). It is the geographic span and volume of this amount of water that results in back-arc tsunamis, and consequently the source parameters contributing the most to back-arc tsunamis are those increasing this volume. Naturally, earthquakes with higher seismic moment (i.e., larger ruptures) or closer to 235 land would create larger back-arc tsunamis as they would displace more water on the opposite side of the continent. However, as large (M>8.5), shallow (H<30 km) thrust earthquakes capable of generating large back-arc tsunamis are mainly expected to happen at or in the vicinity of trenches, the most loosely constrained source parameter will be source geometry. Megathrusts tend to occur at shallow dips ( <30º) and slip vectors almost perpendicular to the trench 240 (45º<λ<135º).
To investigate the influence of source geometry (focal mechanism) and depth, we design and carry out tsunami simulations in a fore-arc-back-arc system as shown in Fig. 4. In these synthetic scenarios, a deep fore-arc (4000 m) is separated from a relatively shallow back-arc 14 of 4.0×10 29 dyn-cm (Mw=9.0) and are placed in the fore-arc, next to the continent (~75 km away) at a central latitude (green stars in Fig. 4). Each scenario is designed as a set of (H, , λ) (i.e., depth, dip, and slip) with 10 km ≤ H ≤ 40 km, 10º ≤ ≤ 40º, and 40º≤λ≤90º (benefiting from the axial symmetry along latitude) while keeping the rupture strike the same as that of the trench, assumed to be parallel to the continental arc. 250 As shown in Figs. 4c and 4d, the spectra of the tsunami in semi-closed (red) back-arc basin are significantly different from those at the fore-arc (blue) gauge, while the difference between the fore-arc spectra in Figs. 4c and 4d is very small. This shows that though the fore-265 and back-arc signals are excited by common sources they have very different frequency contents and thus are essentially different [Rabinovich, 1997]. This is particularly noteworthy in the dominant period, Td, of back-arc tsunamis.
Our simulations show that contrary to fore-arc tsunamis [Rabinovich, 1997;Abe, 2006], the dominant periods of back-arc waves are not directly determined by fault width which is 270 usually taken as a fraction of the fault length. On the other hand, fault dip seems to determine the dominant period (Table S3). This relationship is a quadratic function of dip angle, leading to often longer period signals (Fig. 4e). By controlling the vertical distance to the buried fault in the back-arc, Hb, i.e., a secondary fault depth (Hb=W.sin , with W as fault width), dip angle controls the excitation depth of these waves, and therefore their dominant period. The 275 relationship is further compounded by the effect of dip angle on the geographic extent and thus the perturbed water volume in the back-arc. In other words, the back-arc part of deformation in the direction perpendicular to the arc, can be defined as an effective length controlling the displaced water volume, and hence the tsunami. We also find that the dominant back-arc frequency is invariant of fault slip angle (Table S3). 280 In contrast to fore-arc tsunamis wherein larger dip angles result in higher tsunami amplitudes (Titov et al, 1999), back-arc tsunami amplitudes follow a complex trend. For smaller sources (Mw<9.0), the variations of amplitude against changes in dip angle is not uniform (Fig.   S8), because the back-arc basin samples mostly the smaller, negative part of vertical deformation from the megathrust, mixed by a part (if any) of positive portion of the field, depending on dip 285 angle and source depth (Fig. S9a,b). For larger sources (Mw≥9.0), however, the buried fault approaches the surface and due to larger values of slip, the back-arc experiences larger vertical components (both area and amplitude), accompanied by insignificant changes to the deformation pattern (Figs. 4e,S8,and S9c,d). In fact, simulations show higher back-arc amplitudes for sources with larger dip angle for such large sources (Figs. 4e). 290

Discussion & Conclusion
Our simulations of back-arc tsunamis in the Sea of Japan form the 2011 Tohoku rupture confirm Murotani et al's [2015] conclusion on the significant contribution of horizontal deformation. In addition, we find that the contribution of Rayleigh waves to back-arc tsunami 295 amplitudes, especially in the case of small static deformations, must be taken into account. At nodal stations, addition of surface waves would potentially improve the quality of simulations, especially at higher (although low amplitude) frequencies. This is best demonstrated in the case of Busan which was located at a deformation node of the 2011 Tohoku event.
Numerical simulations of the back-arc tsunamis from sources at and around the Japan 300 trench reveal the potential of back-arc tsunamis with coastal amplitudes reaching more than 1 m in western Japan and up to ~50 cm on the far side of the Japan Sea. While such amplitudes are very small compared to the fore-arc run-up values (>20 m close to the epicenter) and often ignored in the wake of fore-arc catastrophes, they must be taken seriously in a basin which is otherwise deemed as safe from tsunamis. We emphasize that the reported coastal amplitudes 305 from the simulations in this study are at a distance from the coastlines and are thus smaller than run-up [Synolakis, 1991]. In other words, the on-land tsunami amplitudes will be relatively larger than our calculations while following similar along-coast patterns.
Simulation of tsunamis from potential sources in the Nankai Trough leads to similar results, albeit on a smaller scale. We attribute this deficiency of back-arc tsunami amplitudes 310 (compared to the cases in Japan trench) to the wider continental arc as well as a large shallow continental shelf in the southern back-arc. Note that even the largest Nankai scenario, i.e., Fig. 3c-d (model V in Table S1) only succeeds in creating moderate amplitudes on a local scale.
While such amplitudes barely exceed 15 cm, they may pose hazard to the coastal communities in bays and harbors [Abe, 2006]. Similar to the case of ruptures in Japan Trench, small coseismic 315 amplitudes in Korea are due to the close proximity of deformation nodes.
Megathrust events tend to occur in shallow depths (<40 km) with the majority of accumulated strain released within the crust and close to the trench. Although as a global average, trenches are located ~250 km away from continents, depending on the shape, size and physical properties of the megathrust zone, these ruptures can take place only slightly closer to 320 land [Bilek & Lay, 2018]. Therefore, the magnitudes of megathrusts determine the extent of down-dip faulting and thus that of back-arc deformation (see Fig. 1).
At a constant magnitude (Mw=9.0), our results show that among the geometrical source parameters, dip angle ( ) has the most significant contribution and determines the dominant frequency of back-arc tsunamis. We also find that ruptures with larger dip angle create more 325 widespread back-arc tsunamis -this can readily be seen in our simulation of the potential source near Honshu [see Table S2] whose larger dip angle compared to 2011 Tohoku leads to a larger back-arc tsunami (Fig. S11). Thus, large megathrust earthquakes (Mw≥9.0) with large dip angles (>20°) are expected to create larger back-arc tsunamis. Based on this model, the effective length calculated as the extent of back-arc depression in a direction perpendicular to the trench 330 determines the main properties of back-arc tsunamis. While larger ruptures result in slip fields closer to surface, hence producing more significant surface deformation, the focal depth of such events does not seem to have a significant impact (Fig. S8). This simple model can be readily used in subduction zones with back-arc basins around the world. While some regions such as Java seem to be immune to the back-arc tsunami hazard due to the large trench-to-coast distance 335 (an average of ~450 km), such hazard needs to be assessed and mitigated in various regions including Alaska, the Philippines, Sumatra (see supplementary material) and Oaxaca.