Self-similar stochastic slip distributions on a non-planar fault for tsunami scenarios for megathrust earthquakes

Megathrust earthquakes that occur repeatedly along the plate interface of subduction zones can cause severe damage due to strong ground motion and the destructive tsunamis they can generate. We developed a set of scenario earthquakes to evaluate tsunami hazards and tsunami early warning systems for such devastating earthquakes. Although it is known that the slip distribution on a fault strongly affects the tsunami height distribution in near-field coastal areas, the slip distribution of future earthquakes cannot be exactly predicted. One way to resolve this difficulty is to create a set of scenario earthquakes in which a set of heterogeneous slip distributions on the source fault is stochastically generated based on a given slip probability density function (SPDF). The slip distributions generated in this manner differ from event to event, but their average over a large ensemble of models converges to a predefined SPDF resembling the long-term average of ruptures on the target fault zone. We created a set of SPDF-based scenario earthquakes for an expected future Mw 8.2 Tonankai earthquake in the eastern half of the Nankai trough, off southwest Japan, and computed the ensuing tsunamis. We found that the estimated peak coastal amplitudes among the ensemble of tsunamis along the near-field coast differed by factors of 3 to 9, and the earliest and latest arrivals at each observation site differed by 400 to 700 s. The variations in both peak tsunami amplitude and arrival time at each site were well approximated by a Gaussian distribution. For cases in which the slip distribution is unknown, the average and standard deviation of these scenario datasets can provide first approximations of forecast tsunami height and arrival time and their uncertainties, respectively. At most coastal observation sites, tsunamis modeled similarly but using a uniform slip distribution underpredicted tsunami amplitudes but gave earlier arrival times than those modeled with a heterogeneous slip distribution. Use of these earlier arrival times may be useful for providing conservative early warnings of tsunami arrivals. Therefore, tsunami computations for both heterogeneous and uniform slip distributions are important for tsunami disaster mitigation.


Introduction
Megathrust earthquakes along the Nankai trough off southwestern Japan, such as the 1944 Tonankai and 1946 Nankai events, have caused severe damage due to strong ground motion and large tsunamis. Devastating earthquakes such as these have occurred at intervals of 100 to 150 years (e.g., Ando 1975;Ishibashi 2004), so we should be prepared for more of them in the future. Earthquake and tsunami hazard analyses in this region have been done based on multiple sets of scenario earthquakes of various magnitudes and rupture areas (e.g., Earthquake Research Committee, the Headquarters for Earthquake Research Promotion in Japan, 2013, 2020). These analyses have also been used to develop tsunami early warning systems to mitigate the disastrous effects of the tsunamis generated (e.g., Yamamoto et al. 2016b). Goda et al. (2018) carried out probabilistic tsunami hazard analyses in this region based on stochastic earthquake scenarios.
Tsunami warning systems generally estimate coastal tsunami heights and arrival times on the basis of earthquake source parameters, crustal deformation, and initial tsunami heights estimated by inversion of seismic, geodetic, or offshore ocean-bottom pressure gauge (OBP) records, respectively (Tatehata 1997;Titof et al. 2005;Blewitt et al. 2009;Tsushima et al. 2009Tsushima et al. , 2011Tsushima et al. , 2012Tsushima et al. , 2014. Tsunami early warning systems of this type are operated by the Japan Meteorological Agency (JMA) (Tatehata 1997, Tsushima et al. 2012) and NOAA (Titof et al. 2005).
In recent decades, to allow earlier and more precise issuance of earthquake and tsunami warnings, dense offshore seismic and OBP networks have been deployed on the ocean floor above the source areas of expected large earthquakes; these include DONET along the Nankai trough (Kaneda et al. 2015;Kawaguchi et al. 2015) and S-net along the Japan Trench Mochizuki et al. 2016;Uehira et al. 2016). Yamamoto et al. (2016b) developed a method to rapidly estimate initial tsunami heights by using dense OBP networks. The deployment of OBP networks in the near field of seismic sources has accelerated the development of other approaches for tsunami early warnings that need to estimate neither earthquake nor tsunami sources. For example, Maeda et al. (2015) developed a method to use real-time OBP data to estimate coastal tsunami heights based on data assimilation.
Simple methods have also been proposed that estimate coastal tsunami heights by applying amplification factors to offshore tsunami heights recorded at OBP sites (Baba et al. 2014a;Igarashi et al. 2016), or by selecting candidate tsunami sources based on spatial correlations of observed tsunami height distributions at OBP sites with those expected from scenario earthquakes (Yamamoto et al. 2016a;Takahashi et al. 2017Takahashi et al. , 2018. These approaches require a pre-computed database of tsunami parameters for a set of scenario earthquakes in the target area of various magnitudes, source locations, and source geometries. After creation of such a tsunami database, the computational cost of estimating coastal heights and arrival times is lower. Because of the low cost and simplicity of implementation of this type of system, without the need for near-real time source estimations, they are used by some local governments and suppliers of lifeline services (e.g., Takahashi et al. 2017Takahashi et al. , 2018Ishibashi et al. 2018).
The accuracy of tsunami heights estimated using scenario earthquakes depends on the ensemble of scenarios providing adequate coverage of the range of possible future earthquakes. To date, most studies have assumed uniform slip distributions on planar faults in their scenario earthquakes, whereas the spatial distribution of slip is in reality heterogeneous, and this heterogeneity affects coastal tsunami wave height distributions in the near field (Geist 2002;Mueller et al. 2015;Li et al. 2016;Melger et al. 2019). Although it is impossible to predict the slip distributions of future earthquakes, the use of large suites of rapidly produced stochastic slip distributions offers a means to account for some of the inherent uncertainties.
In general, stochastic fault slip distributions can be well approximated by the k-square model, where the wavenumber spectrum decays as k −2 beyond the corner wavenumber k c , which is dependent on earthquake magnitude (Herrero and Bernard 1994;Mai and Beroza 2002;Causse et al. 2010). There are two principal methods used to produce stochastic self-similar slip distributions. One is based on describing the slip spectra in Fourier space, which is then converted to the spatial domain using fast Fourier transforms (Herrero and Bernard, 1994). The alternative, which is the method we used, is the composite source model. It involves the placement of a large number of circular subevents (i.e., patches of slip) whose size-frequency distribution on the fault plane follows a power law (Frankel 1991;Zeng et al. 1994;Ruiz et al. 2011), such that the summation of all the subevents produces a self-similar slip distribution. The placement of these subevents is based on a probability density function (PDF) defined spatially across the fault plane, hereafter referred to as the slip probability density function (SPDF). Individual earthquakes generated using a particular SPDF will differ from event to event, but over a large ensemble of earthquakes, the average slip will converge to the predefined SPDF (Murphy and Herrero 2020). Therefore, the choice of the function used to describe the SPDF depends on a priori knowledge or assumptions of the likely first-order spatial distribution of slip. For example, in situations where there is no knowledge of preferential locations for slip, a uniform PDF could be used , whereas if a specific earthquake or a particular rupture feature is assumed to occur (e.g., shallow slip amplification) a more specific function can be applied (Murphy et al. 2016;Scala et al. 2020).
Our aim in this study was to create a family of megathrust earthquakes that adequately represents possible scenarios that might occur in the Nankai trough. The problem that remained was how to set the SPDF. Recent studies have indicated that the coseismic slip distributions of megathrust earthquakes are possibly correlated with plate interface coupling (Hashimoto et al. 2009(Hashimoto et al. , 2012Perfettini et al. 2010;Moreno et al. 2010), but in some instances, stress accumulated during an interseismic period is partially released, resulting in partial rupture of a strongly coupled region (Konca et al. 2008). These observations suggest that the probability of slip during a megathrust earthquake may be proportional to the interplate coupling. Watanabe et al. (2018) modeled tsunamis generated by a Mw 8.9 earthquake along the Nankai trough by assuming that the slip distribution is proportional to the slip deficit rate (SDR) distribution. Based on this hypothesis, for the SPDF, we used the SDR obtained along the Nankai trough by Yokota et al. (2016). Using this technique, we generated 200 scenario earthquakes for an anticipated M w 8.2 Tonankai earthquake, which would rupture the eastern half of the Nankai trough. By computing the ensuing tsunamis, we evaluated the probability distributions of peak tsunami amplitudes and arrival times along the coast facing the source area to evaluate potential inherent errors in early warnings. We also compared these scenarios with tsunamis computed on the basis of a uniform slip distribution on the same fault area. Our aim was to determine the relative importance of tsunami computations based on heterogeneous and uniform slip distributions for effective tsunami early warnings.
2 Methods/experimental 2.1 Fault geometry of the Nankai trough subduction zone The geometry we used for the subducting Philippine Sea plate is based on that of the Japan Integrated Velocity Structure Model (JIVSM; Koketsu et al. 2012) with the deep plate interface placed at the boundary between oceanic and continental crust (Fig. 1). As a rupture on the plate interface approaches the surface, the rupture plane rises from the top of the oceanic crust and cuts through the overlying sedimentary sequence to reach the trough axis at the seafloor (e.g., Moore et al. 2007Moore et al. , 2009. Following the method of Ide et al. (2010), we represented the plate interface by using an adjustable tension, continuous-curvature, surface-gridding algorithm (Smith and Wessel 1990) to define a surface that smoothly connects the top of the oceanic plate and the trough axis. We represented the non-planar fault of the Nankai subduction zone defined above by using a triangular mesh with nodes at 5 km intervals that extended from the seafloor to a depth of 50 km (Fig. 1).

Stochastic scenarios for a future Tonankai earthquake
As a test case, we generated a set of scenario earthquakes for a future Tonankai earthquake that ruptures the northeastern half of the Nankai trough in an area similar to those of the 1854 and 1944 events (Ando 1975). We assumed an earthquake of M w 8.2 in a source area in the northeastern part of the Nankai trough (Earthquake Research Committee, the Headquarters for Earthquake Research Promotion in Japan, 2013) (Fig. 1). We assumed the SPDF to be proportional to the SDR obtained by Yokota et al. (2016) and set the probability to 0 beyond the southwestern boundary of the source area shown in Fig. 1. We used common rupture area to generate tsunamigenic scenario earthquakes so that we could evaluate the uncertainties of tsunami heights and arrival times in that region without knowledge of the slip distribution on the earthquake fault. We also considered the different effects if the tsunamis were generated by an earthquake with uniform slip on the same rupture area.
To stochastically generate heterogeneous slip distributions on the non-planar fault, we used the method of Herrero and Murphy (2018) but adapted it to facilitate surface ruptures on the seafloor, that is, to allow large fault slip near the trough axis . Because an M w 8.2 earthquake would occur in an area smaller than the entire Nankai subduction zone, we needed to subdivide the meshed fault plane into areas of slip (rupture areas) and no slip. This was done by first randomly choosing a single cell from the SPDF generated from the SDR of Yokota et al. (2016). Cells from the mesh were then added to the rupture area one at a time until the rupture area marginally exceeded the empirical areamagnitude scaling relationship of Strasser et al. (2010). Cells for which at least one of its three sides would be common with a cell of the pre-existing rupture area were selected. When the rupture area had been established, families of circular subevents that had a power-law size distribution, each with an Eshelby slip profile (Eshelby 1957;Ruiz et al. 2011), were placed within the rupture area on the fault plane. The distribution of these subevents was based on the original SPDF after it had been re-normalized over the now-defined rupture area.
To generate stochastic slip distributions, we used two parameter sets: (R max , n a ) = (0.2, 60) and (0.25, 20), where R max is the maximum radius of subevents relative to fault width and n a is the total number of subevents used to generate one slip distribution. We generated 100 scenarios for each parameter set and used the resultant 200 scenarios in our analysis. Figure 2 shows examples of individual slip distributions for the two parameter sets and the average slip distribution for all models. Correlations between the stacks of increasing numbers of slip distributions and the SDR for the two parameter sets (Fig. 3) show that the correlations converge when there are more than about 50 slip distributions. This confirms the validity of our assumption that the average of many slip distributions converges to the SDR distribution and that the ensemble of slip distributions is therefore representative of the SDR.
We also generated a uniform slip distribution for the same rupture area, in which we averaged 100 slip distributions randomly generated from (R max , n a ) = (0.1, 50, 000) with a uniform SPDF.

Tsunami computations
We computed the tsunamis generated by each scenario earthquake. Vertical displacements due to dip-slip motion of the triangular fault elements were computed using the dislocation method of Meade (2007), and the total seafloor displacement for each scenario was computed by summing the contributions from all the elements within the rupture area. Tsunamis due to these static seafloor displacements were computed using the JAGURS code (Baba et al. 2014b). We solved the linear shallow-water equations in a finite difference scheme with a nesting algorithm comprising three layers of which the finest grid spacing was 2 arc-seconds. Although the JAGURS code can also solve the non-linear equations, its computation cost is much higher than that solving the linear one, and we need computations for hundreds of scenario earthquakes. We have confirmed that the discussion does not significantly change if we use the non-linear equation after some preliminary runs. We used bathymetry from the M7000 Digital Bathymetric Chart provided by the Japan Hydrographic Association. The JAGURS code considers tsunami inundation based on the method of Jakeman et al. (2010). We assumed a 40 s rise-time for crustal deformation for M w 8.2 events to initiate the tsunamis. Tsunami observation sites were placed at 2.5 km intervals on the 3-m isobath along the coastline facing the source area (Fig. 4). We included vertical crustal deformations in the tsunami heights because we focused on the near-field effect of the tsunamis. We simulated the tsunami propagations for 2 h with a time step of 0.1 s to satisfy the stability conditions solving for the finite-difference methods.

Results and discussion
We measured the peak near-shore tsunami amplitude (PNTA) and arrival time of each of the tsunamis generated by the 200 scenario earthquakes. We defined the arrival time as the time of arrival of the first tsunami for which the wave-height exceeded 0.2 m above sea level. This threshold is based on the JMA "tsunami advisory," whereby a warning is issued if a tsunami height is expected to exceed 0.2 m above sea level at a target site.
The conditional probability of exceedance of PNTA at each site (Fig. 5 a) and the ratio of the maximum and minimum values of PNTA at each site (Fig. 5 b) indicate that the variations in computed tsunami heights are related to the heterogeneous slip distributions. At most sites, PNTA values differed by factors of 3 to 4, depending on the slip distribution, but they differed by as much as a factor of 9 at some sites. After removing the highest and lowest PNTA values to consider the middle 95% of scenarios, these differences were reduced to a factor of about 3 at most sites ( Fig. 5 b), implying that the largest variations occurred at only a few extreme cases. The distribution of the cumulative probability of exceedance of PNTA at each site (Fig. 6 a) is well approximated by a Gaussian distribution, which indicates that the average and standard deviation can be used to make first approximation estimates of tsunami height and its error, respectively, for early warning systems when the slip distribution is unknown.
The differences between the earliest and latest tsunami arrival times at each observation site varied by 400 to 700 s depending on the scenario (Fig. 7) but, as was the case for PNTA, were reduced (to 300 to 600 s) after removal of the most extreme cases. The distribution of the probability of exceedance at each site is again well approximated by a Gaussian distribution (Fig. 6 b). We also observe some deviations from a Gaussian distribution for the PNTA and arrival times especially at sites located outside the rupture area (site numbers smaller than 50 and site numbers larger than 150 in Fig.  6). Their distributions seem asymmetrical with respect to the average at each site. Other flexible distributions as lognormal, Weibull, Pareto, or others might better fit the probability distribution, but we consider that a Gaussian distribution could be used as a first approximation to the distributions at these sites.
To investigate the magnitude of the effect of slip distribution on PNTA, we compared the scenarios that caused the largest and smallest PNTAs near the Shima Peninsula (location in Fig. 4). The slip distributions that generated both the largest and smallest PNTAs are broadly similar (Fig. 8 a and c), with slip concentrated mainly at a large asperity close to the Shima Peninsula, but the patterns of initial seafloor deformation for these two scenarios differ considerably. The areas of seafloor uplift and subsidence (Fig. 8 b and d) are distributed roughly symmetrically and parallel to the coast for the largest tsunami, whereas for the smallest tsunami the uplift area partially wraps around the area of subsidence.
For the largest tsunami, the propagation of the initial wave is relatively simple. The wave packet of greatest amplitude propagated directly toward the Shima  Peninsula with some focusing due to bathymetry (Fig. 4), which led to the largest tsunami reaching the coast of the Shima Peninsula (Fig. 9 a, Movie01). For the smallest tsunami, the largest amplitude wave packet propagated toward the Enshu Nada coast (location in Fig. 4) east of the Shima Peninsula ( Fig. 9 b, Movie02) because the largest asperity was farther east than that of the largest tsunami (Fig. 8 c). Furthermore, the later coalescence of tsunami waves with positive and negative initial amplitudes west and east, respectively, of the Shima Peninsula (t = 120 s in Fig. 9 b) reduced the tsunami amplitude at the peninsula. These results show that the near-field PNTA distribution can be strongly affected by both fault-slip distributions and interference of tsunami waves coming from different directions. Bathymetry can also affect PNTA distributions. This modeling demonstrates how difficult it can be to estimate PNTA distributions for future tsunamigenic earthquakes and shows that for effective tsunami early warning and disaster prevention it is important to consider a wide variety of scenarios.
Comparison of the PNTAs of the scenario with uniform fault slip distribution and those with heterogeneous slip distributions shows that at almost all observation sites the PNTA for the uniform slip distribution is smaller than those of about 50% of the heterogeneous slip scenarios (Fig. 5 a). At several sites, all the PNTAs for the heterogeneous slip scenarios were greater than that of the uniform slip scenario (e.g., in Fig. 5 a the uniform slip PNTA curve is at or below the 100% probability of exceedance for the heterogeneous case at stations 110-170). This result is consistent with previous studies that have shown that compared with scenarios with heterogeneous slip distributions, a uniform slip distribution mostly underpredicts the near-field PNTA (Geist 2002;Mueller et al. 2015;Li et al. 2016;Melger et al. 2019). Additionally, González et al. (2020) observed that stochastic slip distributions based on the SDR better accounted for historical flow depth measurements observed in Iquique compared to cases where only uniform slip distributions or deterministic slip distributions based solely on the SDR were used. This underprediction is the result of asperities concentrating large slip within limited areas in the source region that in turn excite larger amplitude tsunami waves than would be the case with a uniform slip distribution within the same source area. In contrast, the uniform slip distribution scenario produced relatively small but widely distributed slip in the source area, thus causing smaller initial seafloor displacements over a wider area (Fig. 8 f). Consequently, the tsunami wavefield generated was also wider with smaller amplitude, resulting in smaller PNTAs (Fig. 10,  Movie03).
At most observation sites, the tsunami arrival time for the uniform slip scenario was earlier than those for 50% of the heterogeneous slip scenarios, and earlier than all those scenarios at about 20% of sites (Fig. 7 a). This is because the wider distribution of seafloor uplift for the scenario with uniform slip (Fig. 8 f) resulted in shorter propagation distances and therefore earlier tsunami arrivals at most sites. Accordingly, the earliest arrival times for both the heterogeneous and uniform slip scenarios should be considered in examinations of tsunami hazards. Clearly, near-field PNTA and arrival time distributions depend on the scenarios that were used to generate them. We based all of our scenarios on M w 8.2 events on the same source area in the Tonankai earthquake source region. Coastal PNTA distributions will change if scenarios with different magnitudes and source areas are added. Moreover, we used the SDR distribution of Yokota et al. (2016) to obtain the SPDF we used. Use of different SDRs (e.g., Noda et al. 2018;Nishimura et al. 2018) for the Nankai trough will also affect the source area and slip distributions for the scenario earthquakes. Both Noda et al. (2018) and Nishimura et al. (2018) obtained smaller SDRs beneath the Enshu Nada than did Yokota et al. (2016), the use of which may in turn affect PNTAs along the coast. In the future, these (and perhaps other) SDR distributions, with appropriate weightings, should be incorporated in PNTA and travel time calculations. We used the empirical scaling relationship of Strasser et al. (2010) to establish the rupture area for the slip distributions. However, there are several scaling relationships that could also be applied (e.g., Thingbaijam et al. 2017;Allen and Hayes, 2017). The relationship of Thingbaijam et al. (2017) gives about 10% larger rupture area than that of Strasser et al. (2010), while that of Allen and Hayes (2017) gives about 10% smaller one. This difference would also affect the distributions of PNTA and arrival times. The geometry of the plate interface, especially the shallow part, will also affect the patterns of crustal deformation and the tsunamis they generate. The slip on splay faults can also affect the generation of tsunamis (Park et al. 2002). Further careful investigations of these factors are necessary to properly evaluate tsunami hazards and prepare appropriate ensembles of scenario earthquakes for use in providing effective tsunami early warnings.

Conclusions
We developed a procedure to stochastically generate slip distributions on a non-planar fault surface to generate scenario earthquakes for potential use in tsunami early warning applications. We generated 200 scenario earthquakes for an anticipated M w 8.2 Tonankai earthquake in the Nankai trough and found that the estimated nearfield tsunami amplitudes and coastal arrival times were strongly dependent on slip distribution. The cumulative probability distributions of the variations of both coastal tsunami amplitude and arrival time among the scenarios were well approximated by Gaussian distributions, thus indicating that if slip distributions are unknown the average tsunami amplitude could be used as a first approximation forecast of tsunami height and the standard deviation as its uncertainty. For a uniform slip distribution on the fault, at most observation sites the estimates of peak tsunami amplitudes were underestimated and the arrival times were earlier than those for which a heterogeneous slip distribution scenario was used. Our results highlight the importance of the use of heterogeneous slip distributions in estimations of tsunami amplitude, but also showed that the use of uniform slip distributions may be useful for calculating conservative arrival times for tsunami early warnings. The detailed statistical characteristics of tsunami amplitudes and arrival times clearly depend on assumptions about earthquake magnitude, source region, and the probability density function used to generate slip distributions. Therefore, a comprehensive ensemble of scenarios covering a wide range of earthquake magnitudes and source regions are necessary for accurate evaluations of tsunami hazards and effective tsunami early warning systems.