Uncertainty quantification of tsunami inundation in Kuroshio Town, Kochi Prefecture, Japan using the Nankai-Tonankai megathrust rupture scenarios

. The Nankai-Tonankai megathrust earthquake and tsunami pose significant risks to coastal communities in western and central Japan. Historically, this seismic region hosted many major earthquakes, and the current national tsunami hazard assessments in Japan consider megathrust events having moment magnitudes between 9.0 and 9.1. In responding to the lack of rigorous uncertainty analysis, this study presents an extensive tsunami hazard assessment for the Nankai-Tonankai Trough 15 events, focusing upon the southwestern Pacific region of Japan. A set of 1,000 kinematic earthquake rupture models is generated via stochastic source modelling approaches, and Monte Carlo tsunami simulations are carried out by considering high-resolution grid data of 10-m and coastal defence structures. Significant advantages of the stochastic tsunami simulation methods include the enhanced capabilities to quantify the uncertainty associated with tsunami hazard assessments and to effectively visualize the results in an integrated manner. The results from the stochastic tsunami simulations can inform 20 regional and local tsunami risk reduction actions in light of inevitable uncertainty associated with such tsunami hazard assessments, and complement conventional deterministic tsunami scenarios and their hazard predictions, such as those developed by the Central Disaster Management Council of the Japanese Cabinet Office.


Introduction
The Nankai-Tonankai megathrust earthquake and tsunami pose an imminent threat to people living in western and central 25 Japan. Historically, the Nankai-Tonankai Trough hosted major earthquakes (Ando, 1975;Garrett et al., 2016;Fujino et al., 2018;Fujiwara et al., 2020; see Figure 1). The most recent events in 1944 and 1946 ruptured the eastern (Tonankai) and western (Nankai) parts of the trough, respectively, and similar segmented ruptures occurred during the 1854 Ansei earthquakes on December 23 and 24 with 32 hours apart. The moment magnitudes (M) of these events were in the range between M7.9 and M8.4 (Fujiwara et al., 2020). In contrast, older megathrust events in the Nankai-Tonankai Trough region, 30 such as the 887 Ninna, 1361 Shohei, and 1707 Hoei events, ruptured multiple segments simultaneously based on geological surveys and historical records, resulting in greater moment magnitudes (M8.6 and above; Fujiwara et al., 2020), and caused wider impacts to coastal areas in western-to-central Japan (Figure 1). The average recurrence period of these major events is approximately 150 years (9 events happened over 1336 years since 684 AD). In light of geological evidence (see Garrett et al. [2016] for a comprehensive review) and the 2011 Tohoku earthquake, which was originally thought to rupture in 35 segments but resulted in synchronized rupture of multiple segments (Stein and Okal, 2011), the Central Disaster Management Council (CDMC) of the Japanese Cabinet Office adopted 11 rupture scenarios of M9.0-9.1 that span from the Hyuga-nada to Tokai region (CDMC, 2012; see Figures 1 and 2). These scenarios are used for disaster risk management purposes in Japan, such as developing tsunami inundation maps for municipalities and constructing physical countermeasures. 40 The 2012 CDMC tsunami source models take into account different rupture and asperity patterns across the entire region of the Nankai-Tonankai Trough by expert judgement (Figure 2). The ruptures reach the accretionary wedge of the Nankai-Tonankai Trough, and the largest slip values of circa 52 m are considered. The significant movements of the accretionary wedge at the shallow part of the subduction interface are capable of generating massive tsunamis, as inferred from the 2011 Tohoku earthquake and tsunami (Lotto et al., 2019). The earthquake slip distributions of the 2012 CDMC models are 45 spatially broader than inverted source models for the 2011 Tohoku earthquake and tsunami that show more concentrated slips (Iinuma et al., 2012;Satake et al., 2013), and are consistent with current geodetic estimates of inter-seismic coupling (e.g. Baranes et al., 2018;Watanabe et al., 2018). It is noteworthy that the 2012 CDMC models are intended for representing extreme rupture scenarios at national scale and are not sufficient for encompassing the variability of tsunami hazards at regional and local scales. Such variability is of critical importance when tsunami hazards at municipality/township levels are 50 concerned. In other words, local extreme situations may not be captured by such national-level earthquake rupture scenarios.
These considerations are essential for developing effective tsunami risk reduction and evacuation strategies for local municipalities.
Since one cannot predict how future megathrust events would unfold, evaluating tsunami hazards based on a broad set of possible earthquake rupture scenarios is a viable strategy for better disaster preparedness. Regional and local tsunami risk 55 reduction actions should be informed by rigorous uncertainty analysis (e.g. Fukutani et al., 2015;Mueller et al., 2015;Park et al., 2017). However, there is no study for the Nankai-Tonankai megathrust events that evaluates the variability of tsunami profiles and inundation extents by considering numerous rupture scenarios and high-resolution elevation data. A study by Goda et al. (2018) may be considered as an exception, which employs a stochastic source modelling method of Mai and Beroza (2002) for tsunami hazard assessments of the Nankai-Tonankai megathrust; however, that study is based on a coarse 60 grid resolution of 90 m and 100 source models alone.
This study presents an extensive tsunami hazard assessment for the Nankai-Tonankai Trough events, focusing upon the southwestern Pacific region of Japan. Numerical simulations are conducted for two magnitude ranges M8.9-9.1 and M8.7-8.9, each with 500 rupture models (in total, 1,000 source models). Tsunami simulations are conducted with 90-m and 10-m exposed municipalities to the future Nankai-Tonankai megathrust event. The 10-m grid simulations allow detailed evaluations of tsunami inundation in coastal communities of Kuroshio Town. This is a significant improvement, compared with Goda et al. (2018), whose source parameters were set to imitate the 2012 CDMC models and tsunami simulations were based on 90-m grid resolution. In this study, rupture models for the Nankai-Tonankai earthquake are generated from statistical scaling relationships of earthquake source parameters and stochastic synthesis of earthquake slip distributions 70 (Goda et al., 2016). The newly generated source models are more suitable for capturing the full range of possible earthquake ruptures of the Nankai-Tonankai megathrust, in agreement with a broad set of inverted source models for global subduction earthquakes (Mai and Thingbaijiam, 2014). The tsunami simulations are conducted by considering kinematic rupture processes as well as the effects of defence structures along the coast. In this regard, estimated inundation results are more accurate than other studies that are based on coarse grid resolutions and ignore the coastal defence structures. To benchmark 75 our tsunami simulation results, the 11 rupture models by CDMC (2012) are employed.

Tectonic Characteristics
The Nankai-Tonankai Trough is the primary source of major offshore subduction earthquakes and tsunamis in western and central Japan. The driving tectonic mechanism of these seismic activities is the Philippines Sea Plate subducting underneath 80 the Eurasian Plate with slip rates between 40 and 55 mm/year (Loveless and Maede, 2010;Yokota et al., 2016). A megathrust subduction event occurs when the accumulated slips are released forcefully and is capable of triggering intense ground shaking and a massive tsunami. The potential rupture region of a future event can be inferred based on spatial distributions of slip deficit rates (Yokota et al., 2016;Baranes et al., 2018;Kimura et al., 2019). A large portion of the Nankai-Tonankai Trough has high slip deficit rates (i.e. locked) and thus can be used to determine the future tsunami 85 scenarios (Watanabe et al., 2018).
The Nankai-Tonankai Trough region can be divided into six segments (Garret et al., 2016), as shown in Figure 1. The geological evidence of the Hyuga-nada segment (Z in Figure 1) is limited, and the current seafloor geodetic observations and numerical results indicate relatively low coupling in this segment (Yokota et al., 2016;Kimura et al., 2019). The Nankai segments (A and B in Figure 1) ruptured many times in the past and abundant evidence has shown that they are capable of 90 producing high tsunami waves and large-scale coastal inundations (e.g. Tanigawa et al., 2018). The rupture history of the Tonankai segments (C and D in Figure 1) is also well studied (e.g. Fujino et al., 2018), resulting in widespread shaking and tsunami effects in the central Pacific region of Japan. Seismic imaging surveys conducted by Kodaira et al. (2006) indicated that the shallow fractured portion of the crust off the Kii Peninsula (between the segment B and the segment C) can serve as a physical boundary between the Nankai and Tonankai segments (e.g. 1854 Ansei and1944 and1946 Showa events). On the 95 other hand, the same study also indicated that the deeper portion of the crust underneath the Kii Peninsula is strongly coupled, and thus a large-scale synchronized rupture, such as the 1498 Meio and 1707 Hoei events, is possible, analogously to the case of the 2011 Tohoku event. The Tokai segment (E in Figure 1) ruptured synchronously with the Tonankai segments in the past and generated widespread geological evidence of large seismic events (e.g. Fujiwara et al., 2020).

2012 Central Disaster Management Council Models 100
For tsunami hazard mapping of a future Nankai-Tonankai megathrust event, the CDMC developed 11 tsunami source models by considering that the synchronized rupture over multiple segments is possible and the magnitude of a future Nankai-Tonankai earthquake can be as large as M9.1 (Figure 2). The models are intended to represent extreme scenarios that can be used for coordinating various disaster preparedness actions (i.e. coastal defence structures and evacuation planning).
Since the 2012 CDMC models are relevant to tsunami hazard and risk assessments in western and central Japan and 105 stochastic source models for future Nankai-Tonankai megathrust events are developed in Section 3 by referring to the fault plane geometry of the 2012 CDMC models, their key features are summarized in the following.
The entire fault plane of the 2012 CDMC models is represented by a set of 5,773 sub-faults, each sub-fault having a size of 5 km by 5 km. The fault plane consists of the main subduction interface (5,669 sub-faults) and the Kumano splay fault (104 sub-faults, which are located in the segment C around a depth contour of 5 km), and geometrical parameters (i.e. strike, dip, 110 and rake) of these sub-faults are variable over the curved, steepening fault plane along the dip direction. The total fault-plane area is approximately 1.44×10 5 km 2 . The top-edge depths of sub-faults along the Nankai-Tonankai Trough are set to 0 km (see Figure 1). The shallow parts of the fault plane represent the accretionary wedge of the Nankai-Tonankai Trough. In determining the average slip D over the fault plane via D = (16/7 2/3 )×S 0.5 / (i.e. circular crack model; Madariaga and Ruiz, 2016), where S is the fault-plane area and  is the rock rigidity (= 40.9 GPa), the average stress drop  of 3.0 MPa is 115 used based on the past M8-9 class subduction-type events (i.e. 1944Tonankai, 1946Nankai, 2003Tokachi-oki, 2004Sumatra, 2010Chile, and 2011. The average slips of the 11 source models range between 8.8 and 11.2 m. The moment magnitude is approximately 9.1, with the corresponding seismic moment ranging between 5.3×10 22 and 6.7×10 22 Nm. To determine the slip distributions of the 2012 CDMC source models, two types of asperity, i.e. large and very large, are 120 considered. The large-slip areas take up circa 20% of the entire fault plane with an average slip-value twice as large as the average slip over the fault plane, whereas the very large-slip areas take up circa 5% of the entire fault plane with an average slip-value four times as large as the overall average slip. The large-slip areas are positioned at depths shallower than 20 km, and the very large-slip areas are positioned at depths less than 10 km within the large-slip areas. Among the 11 models ( Figure 2), the models 1 to 5 have a single asperity region (i.e. large and very large-slip areas) with slips concentrated in the 125 segments C-E, B-D, B&C, A&B, and Z-B, respectively. The models 6 and 7 take into account the splay fault in Kumano Sea with the asperity in the segments D-E or in the segment B. The models 8 to 11 consider two asperities in the segments C&E, B&D, A&C, and Z&B, respectively. The rupture nucleation is set to occur near the center of an asperity region at a depth of about 20 km and varies for individual cases. The slip values for individual sub-faults are determined based on spatially-varying convergence rates at the sub-faults given that the total slip amount (or seismic moment) within the large-slip and 130 very large-slip areas is conserved according to the above-mentioned average slip values. In the tsunami source models, a kinematic rupture is represented by a set of earthquake slip distributions (with temporal interval of 10 s) by considering the rupture propagation velocity equal to 2.5 km/s and the rise time of individual sub-faults equal to 60 s. The rupture durations of the 11 models are between 4 and 5 minutes.

Stochastic Tsunami Simulations 135
Stochastic tsunami simulations for the Nankai-Tonankai megathrust events require the generation of multiple earthquake rupture models that are suitable for representing geometrical characteristics and spatial slip distributions for the target region.
The rupture geometry, such as length and width, and average and maximum slips over a fault plane typically scale with earthquake magnitude (e.g. Murotani et al., 2013;Thingbaijiam et al., 2017), whereas the spatial slip distributions can be specified by wavenumber spectral functions (e.g. Mai and Beroza, 2002). To determine the main source parameters, 140 statistical scaling relationships of the eight source parameters (i.e. fault length, fault width, fault area, mean slip, maximum slip, power transformation parameter for the marginal slip values, along-strike correlation length, along-dip correlation length, and Hurst number) that are developed based on a large number of inverted source models of the past major earthquakes can be employed (Goda et al., 2016). Subsequently, stochastic simulation of constrained random slip distributions is performed to generate 500 earthquake rupture models for the two magnitude ranges M8.9-9.1 and M8.7-8.9. 145 The synthesized earthquake rupture models are used to conduct tsunami run-up simulations for coastal communities in Kochi Prefecture and to quantify the anticipated variability of tsunami inundation due to the future tsunamigenic earthquakes. This section outlines numerical methods to perform stochastic tsunami simulations for the Nankai-Tonankai megathrust events.
Because these procedures were explained elsewhere (e.g. Goda et al., 2016Goda et al., , 2018 in detail, explanations are kept succinct.

Stochastic Rupture Model
In generating the stochastic rupture models for the Nankai-Tonankai megathrust events, the 2012 CDMC fault plane ( Figure   1) is adopted as a baseline. Both strike and dip angles of sub-faults, each having a dimension of 5 km by 5 km, are variable and reflect the current tectonic setting in the Nankai-Tonankai Trough region (Hirose et al., 2008). These features are retained in the stochastic rupture models. Specifically, the sub-faults of the 2012 CDMC model are mapped onto a 2D 155 (rectangular) matrix, noting that sub-faults for the Kumano splay fault are excluded. The size of the 2D matrix is 153 (alongstrike) by 53 (down-dip), and its origin is set to the southwestern corner (upper-left) of the fault plane.
Because the fault length L, fault width W, and mean slip D are interrelated for a given seismic moment via Mo = LWD, values of these source parameters and other related earthquake slip distribution parameters, such as maximum slip, power parameter for the marginal slip distribution, along-strike and along-dip correlation lengths, and Hurst number, are generated 160 from the statistical scaling relationships developed by Goda et al. (2016). These parameters are represented by the multivariate lognormal distribution with a correlation structure that is determined based on the existing inverted source models of the major earthquakes (Mai and Thingbaijiam, 2014). We do not allow extreme realizations of these source parameters; for this purpose, if any of the eight source parameters falling outside of an interval between mean plus/minus two standard deviations, values of the source parameters are resampled. Moreover, additional checks of the sampled source parameters are 165 made by ensuring that the aspect ratios of the fault plane are adequate for the Nankai-Tonankai Trough (i.e. falling between 0.75 and 4.0) and the ratios of correlation lengths to fault dimensions are consistent with empirical ranges of 0.15 and 0.6 for the strike direction and 0.15 and 0.45 for the dip direction. Once a suitable fault geometry is determined, the fault plane is placed by randomly selecting a location of the upper-left corner of the simulated fault plane such that it fits within the overall fault plane of the 153-by-53 rectangular matrix. 170 For a given fault plane geometry, the earthquake slip distribution is determined based on marginal slip statistics and spatial slip distribution parameters. A candidate slip distribution is first simulated from an anisotropic 2D von Kármán wavenumber spectrum with its amplitude spectrum being parametrized by along-strike correlation length, along-dip correlation length, and Hurst number and its phase being randomly distributed between 0 and 2 (Mai and Besoza, 2002). The simulated slip distribution is modified via Box-Cox power transformation to achieve a desirable right-skewed feature of the slip marginal 175 distribution (Goda et al., 2016) To ensure that the simulated earthquake slip distribution has realistic characteristics for the Nankai-Tonankai megathrust events, several checks on the simulated slip distribution are conducted. The maximum slip of the simulated earthquake slip distribution is required to be less than 77 m and 57 m for M8.9-9.1 and M8.7-8.9 scenario magnitudes, respectively. The above threshold values of the maximum slip correspond to the mean plus two standard deviations of the predicted maximum slip using the scaling relationship by Goda et al. (2016) for the central magnitude values (i.e. M9.0 and M8.8). These 185 maximum slip values are large but are still consistent with those reported for the 2011 Tohoku earthquake (Iinuma et al., 2012;Satake et al., 2013). In addition, we ensure that the maximum slip does not occur at sub-faults deeper than 20 km, and at least 60% of the total earthquake slip occurs within sub-faults shallower than 10 km ( Figure 1). In other words, the major asperities are constrained to occur in the shallow part of the subduction interface. This is consistent with how asperities are assigned in the 2012 CDMC models (Section 2.2) and is deemed suitable for generating stochastic rupture models for the 190 megathrust subduction events. If the candidate slip distribution does not meet all criteria, this model is discarded, and another slip distribution is generated. This process is repeated until an acceptable source model is obtained.
In modelling the Nankai-Tonankai megathrust events, a kinematic rupture is taken into account (Goda et al., 2018). The probability density functions for hypocenter locations are defined based on the statistical models developed by Mai et al. and the mean/maximum slip ratios. Subsequently, further constraints are placed to exclude unlikely hypocenter locations for a given slip distribution, using empirical findings by Mai et al. (2005). By combining the preliminary probability density functions and the constraints, the final probability density function for hypocenter locations is obtained. This is used to sample the location of hypocenters. Using the randomly generated rupture propagation velocity and rise time, the kinematic rupture process of the synthesized slip distribution can be simulated. In this kinematic rupture modelling, the rupture 200 propagation velocity is characterized by a truncated normal variable with mean equal to 2. Finally, by repeating the above-mentioned procedure for earthquake rupture modelling 500 times for the two magnitude ranges (M8.9-9.1 and M8.7-8.9), a set of 1,000 earthquake rupture models is generated. These rupture models are used to 210 perform Monte Carlo tsunami inundation simulations.

Tsunami Inundation Simulation
For a given earthquake rupture model, the following calculation steps are implemented to evaluate various tsunami hazard characteristics and parameters, such as wave profiles at offshore locations, maximum tsunami heights along coastal lines, and maximum inundation depths at onshore locations (see Section 4). The numerical tsunami inundation analysis follows a 215 standard procedure, namely, computing a vertical dislocation of seawater triggered by an earthquake rupture (Okada, 1985;Tanioka and Satake, 1996) and solving nonlinear shallow water equations for tsunami wave propagation and run-up (Goto et al., 1997).
To represent a computational domain of the tsunami inundation simulations accurately, a complete dataset of bathymetry/elevation, coastal/riverside structures (e.g. breakwater and levees), and surface roughness is obtained from the 220 Japanese Cabinet Office (same as the 2012 CDMC models). The data are provided in the form of six nested grids following a 1/3 nesting ratio rule, i.e. 2430-m -810-m -270-m -90-m -30-m -10-m. The coarse grid regions of 2430-m and 810-m cover the geographical regions of the Nankai-Tonankai Trough and western-to-central Japan (Figure 1), whereas the fine grid regions of 30-m and 10-m cover low-lying land areas along the coast. In this study, two grid levels are focused upon.
For regional offshore tsunami modelling in Kochi Prefecture, nested grid systems up to 90-m resolution are adopted, 225 whereas for local onshore inundation modelling in Kuroshio Town, nested grid systems up to 10-m resolution are considered.
The ocean-floor topography data are based on the 1:50,000 bathymetric charts and JTOPO30 database developed by the Japan Hydrographic Association and the nautical charts developed by the Japan Coastal Guard. The land elevation data are based on the 5-m grid digital terrain model, derived from airborne laser surveys and aerial photographic surveys by the 230 Geospatial Information Authority of Japan. The terrain data are associated with measurement errors of standard deviations less than 1.0 m horizontally and of 0.3 m to 0.7 m vertically. The reference elevation of the bathymetry and terrain data is the standard altitude in Japan (Tokyo Peil), and no variation of sea levels is taken into account in the tsunami simulations. Tokyo Peil is close to the mean sea level and is typically within 0.2 m differences, depending on location.
The elevation data of the coastal/riverside structures are provided by municipalities, supplemented by the national coastline 235 database. In the coastal/riverside structural dataset, structures having dimensions less than 10 m are represented, whereas those having dimensions greater than 10 m are included in the elevation data. In the tsunami simulation, the coastal/riverside structures are represented by vertical walls at one or two sides of the computational cells. To evaluate the volume of water that overpasses these walls, Honma's (1940)  The vertical dislocation profile of seawater due to an earthquake rupture is computed using Okada equations (Okada, 1985), evaluated at 810-m grid resolution. To account for the effects of horizontal movements of steep seafloor on the vertical 245 dislocation of seawater, a method proposed by Tanioka and Satake (1996) is implemented. To alleviate the abrupt changes of the vertical dislocation of seawater, a spatial smoothing filter of 9-by-9 cells is employed, similarly to the 2012 CDMC models. To represent the kinematic rupture process in a tsunami simulation, these vertical dislocation profiles are evaluated at every 10 s and are used as input in the simulation.
The tsunami modelling is carried out using a well-tested TUNAMI code (Goto et al., 1997) that solves the nonlinear shallow 250 water equations using a leap-frog staggered-grid finite difference scheme and is capable of generating offshore tsunami propagation and onshore run-up. The run-up calculation is based on a moving boundary approach (JSCE, 2002), where a dry/wet condition of a computational cell is determined based on total water depth relative to its elevation. The numerical tsunami calculation is performed for a 3-hour duration which is sufficient to model the most critical phase of tsunami waves for the Nankai-Tonankai scenarios. The multi-domain nesting from coarse-to-fine resolution is conducted to consider large-255 to-small scale tsunami waves, depending on water depth. The time-stepping intervals for the regional (90-m) and local (10m) tsunami simulations are set to 0.5 s and 0.1 s, respectively, to satisfy the Courant-Friedrichs-Lewy condition. Moreover, the effects of ground deformation are taken into account by adjusting the elevation data prior to the tsunami simulation run.
However, the effects of tidal variations (e.g. mean high water level) are not considered in the simulations. This section presents uncertainty quantification of tsunami inundation assessments for Shikoku/Kochi Prefecture (regional focus) and for Kuroshio Town (local focus) subjected to stochastic rupture events originated from the Nankai-Tonankai Trough. The magnitude ranges of the stochastic rupture models are M8.9-9.1 and M8.7-8.9 (500 source models for each range) and are consistent with the synchronized ruptures of the historical events and the 2012 CDMC models. Figure 4 shows elevation maps of Shikoku Island (Kochi is a prefecture facing the Pacific) and Ogata/Saga districts (which are local 265 communities in Kuroshio Town; see also Figure 1). In the local maps of the Ogata and Saga districts, the locations of residential buildings are included to indicate the population exposure in these local communities. There are vertical evacuation towers in the Ogata and Saga districts, and their photos are shown in Figures 4d and 4e.
In Section 4.1, offshore wave profiles and maximum tsunami heights of the stochastic tsunami simulations are discussed, in comparison with the counterparts based on the 2012 CDMC models. In Section 4.2, tsunami inundation areas in Shikoku and 270 Kuroshio Town are adopted as regional and local tsunami hazard metrics, respectively, and their relationships with earthquake source characteristics, such as moment magnitude and regional concentration of earthquake slip, are investigated.
In Section 4.3, tsunami inundation characteristics in the Ogata and Saga districts are focused upon. In particular, the severities of local tsunami hazards in these districts are examined by evaluating critical tsunami scenarios.

Offshore Tsunami Profiles and Maximum Tsunami Heights 275
To evaluate the regional tsunami hazard characteristics for Shikoku Island (mainly Kochi Prefecture), offshore tsunami wave profiles at three locations, indicated as P1, P2, and P3 in Figure 4a The largest first tsunami heights reach 16.8 m, 9.5 m, and 16.6 m at P1, P2, and P3, respectively, for the M8.9-9.1 cases. The maximum first tsunami heights for the M8.7-8.9 cases are smaller than those of the M8.9-9.1 cases (approximately 80%), but these waves are still very large and destructive. The simulated wave profiles based on the 2012 CDMC models (blue curves) 285 fall within the ranges of the stochastic simulation results for the M8.9-9.1 cases (grey curves). As expected, the former is closer to the maximum-minimum envelope of the latter. Moreover, the similarity of the stochastic simulation results and the 2012 CDMC model results depends on the locations. For instance, at P1 and P2, these two sets of the modelled wave profiles based on different approaches are similar, whereas they differ significantly at P3 (i.e. the 2012 CDMC simulation results are less than a half of the most extreme case of the stochastic simulation results). Hence, for these offshore locations, the 290 stochastic tsunami simulation results can accommodate those based on the deterministic tsunami scenarios defined by the the 10th and 90th percentile wave profiles of the stochastic simulations (broken red curves) are generally smaller than those based on the 2012 CDMC models (blue curves) (Figure 5), it should not be concluded that the results based on the 2012 CDMC models are more extreme than top 10% of the stochastic simulation results. This is because the phases of the tsunami 295 wave profiles vary and depend on the earthquake rupture characteristics.
Subsequently, to compare the maximum tsunami heights at near-shore locations along the coast of Shikoku Island based on the stochastic tsunami models and the 2012 CDMC models, 154 observation points are set up starting from Sukumo in Kochi Prefecture to Naruto in Tokushima Prefecture (Figure 6a). The extracted near-shore maximum tsunami height results shown in Figure 6 are based on 90-m resolution. The water depths at these observation points are at sea depths of approximately 5 300 to 10 m. Figures 6b and 6c plot the maximum tsunami heights along the coast for the two magnitude ranges M8.9-9.1 and M8.7-8.9, respectively. The plots of the maximum tsunami heights along the coast indicate that the stochastic tsunami

Tsunami Inundation Areas and Earthquake Source Characteristics
The extent of onshore tsunami inundation is a useful tsunami hazard metric for tsunami damage and loss (e.g. Goda et al., 2019). It is thus essential to investigate the characteristics of regional and local inundation area metrics and to relate these to the corresponding earthquake characteristics. In this section, we study the relationships between inundation area parameters 320 and three types of earthquake source parameters. The first source parameter is the moment magnitude, which captures the macroscopic feature of released energy and is proportional to the total slip over a fault plane. The second source parameter is the slip ratio, which is calculated as the summed slip within a specified segment divided by the total slip over the entire fault plane. For instance, the slip ratio can be calculated for individual or combined segments of the Nankai-Tonankai fault plane model (Figure 1). The third source parameter is the tsunami potential energy E T and is defined as (Melgar et al., 2019): 325 where η is the initial sea surface dislocation due to an earthquake rupture, g is the gravitational acceleration, ρ is the density of seawater, and x and y corresponds to the two horizontal coordinates. E T is based on the sum of the squared sea surface dislocation values over the entire computational domain. Therefore, this parameter accounts for the combined effects of various source parameters, such as slip, depth, strike, dip, and rake, on the sea surface dislocation (i.e. tsunami simulation 330 input). Figure 8 shows histograms of tsunami inundation area in Shikoku and Kuroshio Town for the two magnitude ranges M8.9-9.1 and M8.7-8.9. The inundation areas for Shikoku and Kuroshio Town are computed based on 90-m and 10-m resolution grids, respectively, and when grids are flooded by tsunamis greater than 0.1 m depth, these are included in the calculation of inundation area. From regional perspectives (Figures 8a and 8c), inundation areas in Shikoku based on the 2012 CDMC 335 models correspond to higher percentiles of the stochastic tsunami simulation results for the M8.9-9.1 and M8.7-8.9 cases. On the other hand, inundation areas in Kuroshio Town based on the 2012 CDMC models are more consistent with average inundation areas of the stochastic tsunami simulation results.
To examine the relationships between regional/local inundation areas and earthquake slip distributions in segments, tsunami inundation areas in Shikoku and Kuroshio Town are plotted against slip ratios in segments Z (Hyuga-nada), A&B (Nankai), 340 and C-E (Tonankai-Tokai) in Figure 9. It is noted that for a given earthquake source model, the slip ratios in segments Z, A&B, and C-E sum to 1. From the scatter plots, positive and negative trends can be recognized for the segments A&B and C-E, respectively, whereas no clear trends can be identified for the segment Z (see also the Pearson's linear correlation coefficient r indicated in the figure). Although some of the identified trends between the inundation areas and slip rations are weak (e.g. Figures 9a and 9b), the slope coefficients of these relationships are found to be significant (i.e. non-zero) based on 345 the p-values. These observations are consistent with intuitions that when more earthquake slips are concentrated near the target region/location (i.e. Nankai segments A&B for this case), inundation areas become larger. These results demonstrate the importance of taking into account the asperity characteristics of an earthquake rupture in predicting the tsunami impact.
For investigating the effects of spatial earthquake slip distribution on the severity of tsunami hazard, the tsunami potential energy (Equation 1) is an ideal metric because it consolidates the effects of various features of earthquake sources into a 350 single parameter. Figure 10 presents scatter plots of moment magnitude versus tsunami inundation area as well as those of tsunami potential energy versus tsunami inundation area. Figure 10a and 10b exhibit moderate positive relationships between moment magnitude and inundation areas (which can also be inspected in Figure 8). In contrast, Figures 10c and 10d show stronger positive relationships (i.e. linear correlation coefficients of 0.6 to 0.7) between logarithmic tsunami potential energy (with base 10) and inundation areas. The positive correlation is more pronounced for the regional inundation area in Shikoku 355 than for the local inundation in Kuroshio Town because there is a greater chance that elevated tsunami potential energy will affect the target region/location.
Although not presented in this study, similar correlation analyses are also carried out for other earthquake source parameters, such as fault length/width and maximum slip. These results indicate that tsunami potential energy shows the high degrees of correlation with regional and local inundation areas, and thus can be used as an effective tsunami source predictor of regional 360 and local inundation extents.

Tsunami Inundation in Ogata and Saga Districts
Assessing tsunami inundation hazard and quantifying its uncertainty at local community levels require an integrated modelling and analysis of stochastic earthquake sources and tsunami run-up simulations. This section focuses upon: (1) Figure 11 for the two magnitude ranges M8.9-9.1 and M8.7-8.9. In Figure 11, the inundation results based on the 11 CDMC source models are also presented, noting that the vertical evacuation tower in the Saga district is inundated for 3 cases out of the 11 models (i.e. model 4, 5, and 10; see Figure 2) and thus only three square markers are shown in Figures 11b and 11d. Out of all 1,000 scenarios, there are 5 cases and 1 case where the maximum inundation depths at the vertical evacuation towers in the Ogata 380 and Saga districts exceed the critical water depths. The chances of such exceedance are low, and these scenarios can be regarded as very extreme. From tsunami hazard viewpoints, the existing two vertical evacuation towers are judged to be satisfactory. Nevertheless, there are other aspects, such as evacuation times of the local population, and these need to be evaluated to conclude that these towers are sufficient from broader perspectives. This is beyond the scope of the current study. 385 Next, to derive critical tsunami hazard scenarios based on local inundation areas, inundation areas in the Ogata and Saga districts are investigated, and histograms of these inundation parameters are shown in Figure 12. Note that the areas that are considered for these local communities are smaller than those considered for Kuroshio Town (i.e. Figures 8b and 8d). Based on the histograms shown in Figure 12 and by taking the 50th and 90th percentiles as critical scenario levels (note: other percentiles can be adopted), the critical inundation areas are obtained as 2.20 and 3.71 km 2 , respectively, for the M8.9-9.1 390 cases and as 1.80 and 3.34 km 2 , respectively, for the M8.7-8.9 cases. For comparison, the 2012 CDMC model 5 is adopted as the most critical scenario among the 11 models shown in Figure 2. The inundation area in the Ogata and Saga districts is obtained as 3.33 km 2 .
Once the earthquake rupture models that correspond to the identified critical inundation areas, various stochastic tsunami simulation results, such as regional maximum tsunami heights and local maximum inundation depths, can be extracted. For 395 such purposes, the earthquake slip distribution, the maximum tsunami height distribution, the maximum inundation depth in the Ogata district, and the maximum inundation depth in the Saga district for the 50th and 90th percentile inundation levels are displayed in Figure 13 for the M8.9-9.1 cases and in Figure 14 for the M8.7-8.9 cases. As a benchmark for comparison, the counterpart results based on the 2012 CDMC model 5 are included in both figures. The results shown in Figure 13 indicate that the extents of local tsunami inundation in the Ogata and Saga districts for the 2012 CDMC model 5 are between 400 those for the 50th and 90th percentile scenarios of the M8.9-9.1 cases. On the other hand, for the cases of M8.7-8.9, the results based on the 2012 CDMC model 5 exceed or close to those based on the 90th percentile critical scenario. By inspecting the corresponding earthquake slip distribution as well as maximum tsunami height distribution (i.e. 1st and 2nd rows of Figures 13 and 14), the local inundation results can be better related to earthquake source and asperity characteristics and regional tsunami wave characteristics. For instance, the 50th percentile scenario is less critical for the Ogata and Saga 405 districts because it is a two-asperity earthquake slip distribution in the segments Z&A and D (Figure 13b), whereas the 90th percentile scenario causes severe consequences in the Ogata and Saga districts because its asperity is concentrated in the segments A&B (Figure 13c), similarly to the 2012 CDMC model 5 (Figure 13a), although detailed slip patterns shown in Figures 13a and 13c are different.
Significant advantages of the stochastic tsunami simulation methods and their use in deriving critical tsunami scenarios and 410 related tsunami hazard maps are the enhanced capabilities to quantify the uncertainty associated with tsunami hazard assessments and to visualize the results in an integrated manner more effectively. It is also important to point out that the deterministic tsunami scenarios, such as the 2012 CDMC models, and the stochastic tsunami scenarios, such as the 1,000 rupture models, are complementary. Typically, deterministic scenarios are derived based on current geodetic and available geological data, whereas stochastic scenarios are more inclined to statistical features of earthquake source models of the past 415 events. The systematic comparison of these results will improve the understanding of different modelling approaches and their key assumptions and will allow hazard modellers to gain confidence in the derived results.

Conclusions
Tsunami hazard assessments based on stochastic tsunami simulations offer valuable insights into the degree of uncertainty associated with such investigations. To evaluate tsunami inundation hazard and quantify its variability at local community 420 levels, this study presented an integrated modelling and analysis of stochastic earthquake sources and tsunami run-up simulations for the future Nankai-Tonankai megathrust events. For this purpose, 1,000 kinematic tsunami rupture models were generated for the moment magnitudes ranging from 8.7 to 9.1, and tsunami simulations were performed with highresolution grid data of 10-m by taking into account the effects of existing coastal defence structures in Japan. To benchmark the results from the stochastic tsunami simulations, a set of tsunami source models developed by the CDMC of the Japanese 425 Cabinet Office was employed. Our motivations in comparing the stochastic tsunami inundation results with the deterministic 2012 CDMC models were to quantify the variability of tsunami hazards at municipality/township levels and to account for local extreme situations. To enable a consistent comparison with the CDMC tsunami source models (M9.0 to M9.1) in terms of released seismic moment, the lower magnitude limit of the stochastic sources was chosen as M8.7.
The numerical investigations focused upon the southwestern Pacific region of Japan, i.e. Shikoku Island and Kuroshio Town 430 from regional and local viewpoints. The comparisons of the offshore wave profiles and maximum tsunami heights along the coast of Shikoku Island based the stochastic tsunami simulations and the 2012 CDMC models indicated that the 2012 CDMC results are consistent with the typical stochastic simulation results for the same magnitude range but are unable to capture extreme scenarios of local tsunami hazards and their variability.
To relate the regional and local inundation extents to earthquake source characteristics, correlations between inundation area 435 metrics and moment magnitude, slip ratio in segments, or tsunami potential energy were examined. The results indicated that tsunami potential energy could be used as an effective tsunami source predictor of regional and local inundation extents.
Moreover, to evaluate the sufficiency of the two existing vertical evacuation towers in the local communities of Ogata and Saga in Kuroshio Town, the critical inundation depths of these towers were compared to the stochastic tsunami simulation results. Since the exceedance of the critical inundation depths was rare (i.e. 5 and 1 out of 1,000 cases for the towers in the 440 Ogata and Saga districts, respectively), the existing two vertical evacuation towers were judged to be satisfactory. Finally, critical tsunami scenarios for the local communities in Ogata and Saga were identified based on the local inundation metrics and their tsunami inundation extents were compared with those of the 2012 CDMC models. The use of stochastic tsunami simulation methods improved the quantification and visualization of uncertain tsunami hazard assessments and was able to complement the tsunami hazard predictions based on conventional deterministic tsunami scenarios. Future investigations 445 should extend the hazard assessments into tsunami evacuation problems as well as into multi-hazard risk assessments.          Shikoku and Kuroshio Town for the two magnitude ranges M8.9-9.1 (c) and M8.7-8.9 (d). Figure  4) for the two magnitude ranges M8.9-9.1 (a,b) and M8.7-8.9 (c,d). The sum of the vertical bin heights is 1. Figure 12: Histograms of inundation areas in the Ogata and Saga districts for the two magnitude ranges M8.9-9.1 (a) and M8.7-8.9 (b). In each figure panel, inundation areas that are used to define local critical tsunami scenarios are indicated by the 50th and 90th percentile lines. The sum of the vertical bin heights is 1.