Probabilistic tsunami inundation assessment of Kuroshio Town,Kochi Prefecture, Japan considering the Nankai-Tonankai megathrust rupture scenarios

Abstract. 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 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 defense 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 regional and local tsunami risk reduction actions in light of inevitable uncertainty associated with such probabilistic 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 M8.1 and M8.4. In contrast, older megathrust events in the Nankai-Tonankai Trough region, such as the 887 Ninna, 30 https://doi.org/10.5194/nhess-2020-169 Preprint. Discussion started: 8 June 2020 c Author(s) 2020. CC BY 4.0 License. 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), and caused wider impacts to coastal areas in western-tocentral 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 segments but resulted in synchronized rupture of 35 multiple segments, 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 Figure 1 and Figure 2). These scenarios are used for disaster risk management purposes in Japan, such as developing tsunami inundation maps for municipalities and constructing physical countermeasures.
The 2012 CDMC tsunami source models take into account different rupture and asperity patterns across the entire region of 40 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 spatially broader than inverted source models for the 2011 Tohoku earthquake and tsunami that show more concentrated 45 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 concerned. In other words, local extreme situations may not be captured by such national-level earthquake rupture scenarios. 50 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 reduction actions should be informed by rigorous uncertainty analysis (e.g. Fukutani et al., 2015;Mueller et al., 2015;Park 55 et al., 2017). However, there is no study for the Nankai-Tonankai megathrust events that evaluates the variability of tsunami profiles and inundation extents events 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 grid resolution of 90 m and 100 source models alone. 60 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.7-8.9 and M8.9-9.1, each with 500 rupture models (in total, 1,000 source models). Tsunami simulations are conducted with 90-m and 10-m grid resolutions for Kochi Prefecture and Kuroshio Town, respectively ( Figure 1). Kuroshio Town is one of the most https://doi.org/10.5194/nhess-2020-169 Preprint. Discussion started: 8 June 2020 c Author(s) 2020. CC BY 4.0 License. exposed municipalities to the future Nankai-Tonankai megathrust event. The 10-m grid simulations allow detailed 65 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 probabilistic scaling relationships of earthquake source parameters and stochastic synthesis of earthquake slip distributions (Goda et al., 2016). The newly generated source models are more suitable for capturing the full range of possible earthquake 70 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 defense structures. To benchmark our tsunami simulation results, the 11 rupture models by CDMC (2012) are employed. 75 2 Nankai-Tonankai Megathrust Earthquake and Tsunami

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 the Eurasian Plate with slip rates between 40 and 55 mm/year (Loveless and Maede, 2010;Yokota et al., 2016). A 80 megathrust subduction event occurs when the accumulated slips are released forcefully and triggers 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 scenarios (Watanabe et al., 2018). 85 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 producing high tsunami waves and large-scale coastal inundations (e.g. Tanigawa et al., 2018). The rupture history of the 90 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 offshore 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 other hand, the same study also indicated that the deeper portion of the crust underneath the Kii Peninsula is 95 strongly coupled, and thus a large-scale synchronized rupture, such as the 1498 Meio and 1707 Hoei events, is possible, https://doi.org/10.5194/nhess-2020-169 Preprint. Discussion started: 8 June 2020 c Author(s) 2020. CC BY 4.0 License.
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 defense 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 asperities, 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-https://doi.org/10.5194/nhess-2020-169 Preprint. Discussion started: 8 June 2020 c Author(s) 2020. CC BY 4.0 License.
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 fault 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 probabilistic 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 synthesis 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 54 (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 https://doi.org/10.5194/nhess-2020-169 Preprint. Discussion started: 8 June 2020 c Author(s) 2020. CC BY 4.0 License. from the probabilistic 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). At this stage, a taper function is applied to deeper segments of the fault plane, similar to the 2012 CDMC models. In particular, slip values for sub-faults having top-edge depths of 25-27 km, 27-29 km, 29-31 km, 31-33 km, and 33-35 km are multiplied by the reduction factors of 0.9, 0.7, 0.5, 0.3, and 0.1, respectively. Subsequently, slip values of the sub-faults that are not within the original 2012 CDMC models are set to zero, and all eligible slip values are scaled to match the mean slip. 180 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. https://doi.org/10.5194/nhess-2020-169 Preprint. Discussion started: 8 June 2020 c Author(s) 2020. CC BY 4.0 License. (2005). The preliminary probability density functions for hypocenter locations are specified based on the fault dimensions 195 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.5 km/s, standard deviation equal to Finally, by repeating the above-mentioned procedure for earthquake rupture modelling 500 times for the two magnitude ranges (M8.7-8.9 and M8.9-9.1), a set of 1,000 earthquake rupture models are generated. These rupture models are used to perform Monte Carlo tsunami inundation simulations.

Tsunami Inundation Simulation 210
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 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 215 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 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 220 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. 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 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 245 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 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 250 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 largeto-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. 255

Results
This section presents probabilistic tsunami inundation assessment results for Shikoku/Kochi Prefecture (regional focus) and for Kuroshio Town (local focus) subjected to stochastic rupture events originated from the Nankai-Tonankai Trough. The moment magnitude ranges of the stochastic rupture models are M8.9-9.1 and M8.7-8.9 (500 source models for each range) https://doi.org/10.5194/nhess-2020-169 Preprint. Discussion started: 8 June 2020 c Author(s) 2020. CC BY 4.0 License. and are consistent with the synchronized ruptures of the historical events and the 2012 CDMC models. Figure 4 shows 260 elevation maps of Shikoku Island (Kochi is a prefecture facing the Pacific) and Ogata/Saga districts (which are local communities in Kuroshio Town). 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 265 comparison with the counterparts based on the 2012 CDMC models. In Section 4.2, tsunami inundation areas in Shikoku and 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. 270

Offshore Tsunami Profiles and Maximum Tsunami Heights
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, are examined in Figure 5. These locations are near Shimanto City, Kochi City, and Cape Muroto and are at sea depths of 36 m, 45 m, and 70 m, respectively. For each of the two magnitude ranges M8.9-9.1 and M8.7-8.9, tsunami wave profiles based on 500 stochastic models are shown with grey 275 colors, the 50th percentile curve is shown with a solid red line, the 10th and 90th percentile curves are shown with broken red lines, whereas the simulated wave profiles based on the 2012 CDMC models are shown with blue curves.
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) 280 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 285 stochastic tsunami simulation results can accommodate those based on the deterministic tsunami scenarios defined by the CDMC. This observation is generally applicable to other offshore locations that are inspected as part of this study. Although 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 290 wave profiles vary and depend on the earthquake rupture characteristics.

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 315 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 within a specified portion of the fault plane with respect to the total slip. 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 ET and is defined as (Melgar et al., 2019): 320 https://doi.org/10.5194/nhess-2020-169 Preprint. Discussion started: 8 June 2020 c Author(s) 2020. CC BY 4.0 License.
where η is the initial sea surface dislocation due to an earthquake rupture, g the gravitational acceleration, ρ the density of seawater, and x and y corresponds to the two horizontal coordinates. ET 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 325 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 grids based on 90-m and 10m resolution, 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 330 2012 CDMC 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. Furthermore, reasonable degrees of similarity can be observed between the maximum tsunami heights along the coast of Kochi and Tokushima Prefectures (Figures 6 and 7) and the histograms of tsunami inundation area in Shikoku and Kuroshio Town (Figure 8). 335 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), 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 340 coefficient r indicated in the figure). 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 345 energy (Equation 1) is an ideal metric because it consolidates the effects of various features of earthquake sources into a 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, Figure 10c and 10d show stronger positive relationships (i.e. linear correlation coefficients of 0.6 to 0.7) between logarithmic tsunami potential energy 350 (with base 10) and inundation areas. The positive correlation is more pronounced for the regional inundation area in Shikoku 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 355 https://doi.org/10.5194/nhess-2020-169 Preprint. Discussion started: 8 June 2020 c Author(s) 2020. CC BY 4.0 License. correlation with regional and local inundation areas, and thus can be used as an effective tsunami source predictor of regional 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)  evaluating the sufficiency of these two towers as vertical refugee facilities. To carry out such assessments, histograms of maximum inundation depth at the vertical evacuation towers in the Ogata and Saga districts are shown in Figure 11 for the two magnitude ranges M8.9-9.1 and M8.7-8.9. 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 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, 375 the exiting 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.
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 380 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 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 385 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 https://doi.org/10.5194/nhess-2020-169 Preprint. Discussion started: 8 June 2020 c Author(s) 2020. CC BY 4.0 License. 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 390 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 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 395 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 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 400 segments A&B (Figure 13c Significant advantages of the stochastic tsunami simulation methods and their use in deriving critical tsunami scenarios and 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 405 deterministic tsunami scenarios, such as the 2012 CDMC models, and the probabilistic tsunami scenarios, such as the 1,000 stochastic rupture models, are complementary. Typically, deterministic scenarios are derived based on current geodetic and available geological data, whereas probabilistic scenarios are more inclined to statistical features of earthquake source models of the past 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. 410

Conclusions
Probabilistic 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 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 415 models were generated for the moment magnitudes ranging from 8.7 to 9.1, and tsunami simulations were performed with high-resolution 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 Cabinet Office was employed.
The numerical investigations focused upon the southwestern Pacific region of Japan, i.e. Shikoku Island and Kuroshio Town 420 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 425 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 430 Ogata and Saga districts, respectively), the exiting 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 435 should extend the hazard assessments into tsunami evacuation problems as well as into multi-hazard risk assessments.