Tsunami simulations of mega-thrust earthquakes in the Nankai–Tonankai Trough (Japan) based on stochastic rupture scenarios

Abstract In this study, earthquake rupture models for future mega-thrust earthquakes in the Nankai–Tonankai subduction zone are developed by incorporating the main characteristics of inverted source models of the 2011 Tohoku earthquake. These scenario ruptures also account for key features of the national tsunami source model for the Nankai–Tonankai earthquake by the Central Disaster Management Council of the Japanese Government. The source models capture a wide range of realistic slip distributions and kinematic rupture processes, reflecting the current best understanding of what may happen due to a future mega-earthquake in the Nankai–Tonankai Trough, and therefore are useful for conducting probabilistic tsunami hazard and risk analysis. A large suite of scenario rupture models is then used to investigate the variability of tsunami effects in coastal areas, such as offshore tsunami wave heights and onshore inundation depths, due to realistic variations in source characteristics. Such investigations are particularly valuable for tsunami hazard mapping and evacuation planning in municipalities along the Nankai–Tonankai coast.

Mega-thrust subduction earthquakes may trigger massive tsunamis that can lead to catastrophic consequences for coastal cities and towns. Recent devastating tsunamis include the 2004 Indian Ocean tsunami (Murata et al. 2010) and the 2011 Great East Japan (Tohoku) tsunami (Fraser et al. 2013). Globally, tsunami exposure is not negligible and more accurate assessment of tsunami hazard due to future mega-thrust subduction earthquakes is necessary to mitigate the potential tsunami risk and to enhance the tsunami resilience of coastal communities (Løvholt et al. 2014). One of the most challenging problems in tsunami preparedness and risk mitigation is to determine the range of possible tsunami scenarios for a given geographical area. Methods, such as probabilistic tsunami hazard analysis (Geist & Parsons 2006;Gonzalez et al. 2009;Horspool et al. 2014), are useful for identifying tsunami source regions and corresponding scenarios that have major impact to a site of interest. In developing such scenarios, many experts with different backgrounds are involved and the scientific views that are expressed by these experts may be diverse. Such expert judgements rarely lead to a consensus model, but rather generate a set of disparate scenarios that need to be weighted in a logic-tree approach. Nevertheless, in face of imminent tsunami risk in major subduction zones, scientists and engineers need to evaluate the risk quantitatively and effectively, and produce reliable and testable hazard/risk maps at national/regional/local levels. In the ideal (future) case, such hazard and risk maps do not depend on difficult-to-quantify expert solicitations, but on algorithmically developed and objective selection criteria. Typically, such testable hazard quantification will require a probabilistic/ ensemble simulation-based approach.
Recent advances in tsunami hazard analysis aim at incorporating the 'anticipated variability' associated with source characteristics of future tsunamigenic earthquakes, such as location, magnitude, geometry and slip distribution (Geist 2002;McCloskey et al. 2008;Løvholt et al. 2012;Goda et al. 2014;Fukutani et al. 2015;Mueller et al. 2015).
Methods for developing stochastic earthquake source models, based on spectral analysis of slip heterogeneity of inverted source models and subsequent spectral synthesis of constrained-random slip fields (Mai & Beroza 2002;Lavallée et al. 2006), facilitate the generation of possible scenarios with different earthquake slips and fault geometry. The approaches also allow for incorporating epistemic uncertainty in the stochastic source characterization by considering one or more reference source models, if available (Goda et al. 2014). Through Monte Carlo tsunami simulation, probabilistic inundation depth maps can be evaluated using numerous synthesized source models (Mueller et al. 2015). Furthermore, by integrating tsunami hazard estimates with tsunami fragility models (Tarbotton et al. 2015), probabilistic tsunami risk maps, as well as tsunami loss curves, at local and regional scales can be developed (Goda & Song 2016). An advantage of such hazard and risk assessments is that the currently known sources of uncertainty related to the tsunami source characteristics are taken into account, whereby such simulation-based approaches also provide a basis to develop testable and easily updatable hazard maps.
Historically, the Nankai-Tonankai Trough hosted many M w 8+ mega-thrust subduction earthquakes that triggered intense ground shaking and massive tsunami across the central to SW part of Japan (Ando 1975). In this region, the Philippines Sea Plate subducts underneath the Eurasian Plate with an average slip rate of 44 mm a 21 (Loveless & Meade 2010). The source areas span over offshore regions of Nankai (Kyushu-Shikoku-Kinki region) and Tonankai (Kinki -Tokai region), as shown in Figure 1. The magnitude of a future Nankai-Tonankai earthquake can be as large as M w 9.0 (if the entire trough region ruptures in a single event), which is similar in size to the 2011 Tohoku earthquake. The Central Disaster Management Council (2012), under the auspices of the Cabinet Office of the Japanese Government, has developed a new tsunami source model for the Nankai-Tonankai earthquake based on expert judgement, and has produced updated estimates of tsunami hazard, risk and economic impact due to such a catastrophic event. Notable features of the Nankai-Tonankai tsunami source model by the Central Disaster Management Council (hereafter, the 2012 CDMC model) include: (i) consideration of 'large-slip areas' and 'very-large-slip areas' (socalled asperities); and (ii) consideration of multiple scenarios (in total, 11 cases) to represent different slip characteristics in terms of the number and locations of very large asperities. The large-slip areas are placed in the shallow part of the source region (depths of less than 20 km) and consist of about 20% of the entire fault plane, whereas the verylarge-slip areas are located within the large-slip areas, consisting of about 5% of the entire fault plane. For the Nankai-Tonankai scenarios, assigned average slip values are about 10 m, whereas the large-slip areas and very-large-slip areas typically show 20 and 40 m fault displacement, respectively. These slip parameters reflect the characteristics of inverted source models for the 2011 Tohoku earthquake (e.g. Yamazaki et al. 2011;Gusman et al. 2012;Satake et al. 2013). Although the locations of asperities are varied among the 11 cases, it is unlikely that the considered slip distributions comprehensively cover possible future scenarios. Consequently, they are probably insufficient to fully quantify the uncertainties associated with tsunami hazard and risk predictions for the Nankai-Tonankai event (note: the 2012 CDMC model was not intended for modelling/evaluating such uncertainties). Nonetheless, such assessments are crucial for improving disaster preparedness by understanding the consequences of different scenarios, and by communicating uncertainties associated with hazard predictions among stakeholders (Goda & Song 2016).
This study develops constrained-stochastic earthquake slip models for the Nankai-Tonankai event. The models are based on the approach of generating stochastic source models for the 2011 Tohoku earthquake by Goda et al. (2014), but in addition reflecting key features of the 2012 CDMC tsunami source model. One of the innovative features of the newly developed stochastic tsunami source models is the incorporation of the kinematic fault rupture evolution: that is, specifying the point of rupture nucleation, the rupture propagation speed and the slip duration (rise time) on the fault. These models can then be used to carry out detailed analyses on tsunami sensitivity due to varying source characteristics for the Nankai-Tonankai event, as well as for the resulting hazard and risk predictions. Such investigations, albeit challenging, are particularly valuable for tsunami hazard mapping and evacuation planning in municipalities along the Nankai -Tonankai coast.
This paper is organized as follows. First, we briefly summarize geological aspects of the Nankai -Tonankai subduction region and the 2012 CDMC tsunami source model. Subsequently, details of the stochastic Nankai-Tonankai tsunami source model are described, focusing on the slip distribution generation and kinematic rupture modelling. Using the constrained-stochastic source models for the Nankai -Tonankai earthquake, we conduct Monte Carlo tsunami simulation. Variability of tsunami wave profiles and tsunami inundation depths due to varying source characteristics is evaluated to discuss the extent of uncertainty associated with tsunami hazard predictions for the Nankai-Tonankai tsunami. Finally, we discuss the main conclusions of this study, as well as future investigations that need to be conducted using our newly developed tool.
Nankai -Tonankai mega-thrust subduction earthquake Nankai -Tonankai subduction region The Nankai-Tonankai Trough accommodates plate movements between the Philippines Sea Plate and the Eurasian Plate. The convergence rate varies along the trough. Recent geodetic measurements in Japan indicate that the slip rate increases from about 36 mm a 21 at the eastern end near Suruga Bay to about 63 mm s 21 at the western end near Kyushu (Loveless & Meade 2010) (Fig. 1). The accumulated strain is released periodically during major subduction earthquakes. Since 1700, five M8+ earthquakes occurred (Ando 1975). The most recent ones were the 1944 M w 8.1 Tonankai earthquake and the 1946 M w 8.4 Nankai earthquake. In 1854, large dual earthquakes (Ansei earthquakes, M w 8.4 and M w 8.5) occurred in sequence, separated by 32 h (note: the eastern Tonankai segment ruptured prior to the western Nankai segment). Further dating back to 1707, the M w 8.7+ Hoei earthquake ruptured the entire Nankai-Tonankai region (where it is possible that the actual event occurred as an earthquake doublet with a short time gap : Ando 1975). By examining the earthquake rupture patterns of past Nankai-Tonankai earthquakes since 684, Ando (1975) suggested two predominant rupture modes: segmented ruptures, in which only a single segment fails; and synchronized (or coupled) rupture modes in which two segments fail in rapid succession. The recurrence interval of large subduction events in the Nankai -Tonankai Trough varies from 90 to 264 years (Ando 1975).
Issues related to fault zone segmentation are crucial for earthquake and tsunami hazard assessments in general, but it appears that this is particularly important for the Nankai-Tonankai events. Seismic imaging surveys conducted by Kodaira et al. (2006) indicated that a physical segmentation boundary exists ( just off the Kii Peninsula: Fig. 1) that separates the Nankai and Tonankai source regions, and thus controls the occurrence of segmented or synchronized rupture. The shallow portion of the crust offshore the Kii Peninsula is fractured, while the deeper portion of the crust is strongly coupled. When the rupture starts near the strongly coupled part, the synchronized rupture may be triggered (Kodaira et al. 2006), affecting a wider geographical region simultaneously. Considering the large-scale synchronized ruptures of the 2004 Sumatra earthquake and the 2011 Tohoku earthquake, it is difficult to predict whether the next mega-thrust earthquake along the Nankai-Tonankai Trough will rupture in segments or simultaneously, as this region has shown events occurring in both rupture modes in the historic past.

Central Disaster Management Council tsunami source model
The 2012 CDMC tsunami source model is developed for the synchronized rupture case. The model is intended to represent one of the more extreme scenarios that may need to be considered from an earthquake disaster mitigation viewpoint. Therefore, it corresponds to a M w 9-class event and its rupture area spans over a vast region (Fig. 1). The key features of the 2012 CDMC model include: † The entire fault plane is represented by a set of 5773 sub-faults (including the Kumano splay fault: Fig. 1); the size of each sub-fault is 5 × 5 km. The main fault plane consists of 5669 sub-faults, while the splay fault plane is composed of 104 sub-faults. The total faultplane area is about 1.4 × 10 5 km 2 . The top-edge depth for sub-faults along the Nankai-Tonankai Trough is set to 0 km or very shallow depths (i.e. the rupture plane reaches or is very close to the ocean bottom). † The average slip (D) over the fault plane is determined based on a scaling law: D ¼ (16/7p 2/3 ) DsS 0.5 /m, where Ds is the average stress drop, S is the fault-plane area and m is the rock rigidity (¼40.9 GPa). The representative average stress drop is determined by analysing stress drop parameters estimated for the past M w 8-9-class subduction-type events (i.e. 1944Tonankai, 1946Nankai, 2003Tokachi-oki, 2004Sumatra, 2010 Chile and 2011 Tohoku). The mean plus 1 SD (standard deviation) of the stress drop parameters from the considered inversion studies is 2.2 MPa (Central Disaster Management Council 2012): whereas, for smaller subduction events (M w 7-8), a representative stress drop value is about 3.0 MPa. In the 2012 CDMC model, Ds is set to 3.0 MPa. † Two types of asperities are defined. The largeslip areas take up about 20% of the entire fault plane, with an average slip-value twice as large as the average slip over the fault plane. The large-slip areas are positioned at depths shallower than 20 km. The very-large-slip areas take up about 5% of the entire fault plane,with an average slip-value four times as large as the average slip over the fault plane. Note that the very-large-slip areas are located within the largeslip areas and are constrained to be at depths less than 10 km. † There are 11 scenarios (see Fig. 2). Cases (1) -(5) consider a single asperity region (i.e. large and very large-slip areas): the asperity is positioned off (1) Suruga Bay to Kii Peninsula, (2) Kii Peninsula, (3) Kii Peninsula to Shikoku, (4) Shikoku, or (5) Shikoku to Kyushu. Cases (6) and (7) take into account local splay fault sources in the Kumano Sea (see Fig. 1); for these cases, asperities are located off the Kii Peninsula (either on the western or eastern side of the Kumano Sea). Cases (8) -(11) consider two asperities; their locations are off: (8) Suruga Bay to Kii Peninsula and Kii Peninsula to Shikoku, (9) Suruga Bay to Kii Peninsula and Shikoku, (10) Kii Peninsula to Shikoku and Shikoku to Kyushu, or (11) Shikoku and Shikoku to Kyushu. Note that rupture nucleation is set to occur near the centre of asperity regions at a depth of about 20 km, and varies for individual cases (therefore, the total rupture duration for the 11 cases also differs). † The slip values for individual sub-faults are determined based on spatially varying convergence rates at the sub-faults given that the total slip amounts (or seismic moments) within the large-slip and very-large-slip areas are conserved according to the above-mentioned average slip values. † The average slip values for the 11 cases range between 8.8 and 11.2 m, whereas the maximum slip values for the 11 cases range between 39.8 and 52.0 m, increasing from east to west (Fig. 2). † The moment magnitude is about 9.1, and the corresponding seismic moment ranges between 5.3 × 10 22 and 6.7 × 10 22 N m (for the 11 cases). † The geometrical parameters (i.e. strike, dip and rake) are variable over the curved fault plane (approximated by the sub-faults). † The kinematic slip distributions are generated at an interval of 10 s. The assumed rupture propagation velocity is 2.5 km s 21 and the rise time of individual sub-faults is set to 60 s.

Methodology
The stochastic source models for the Nankai -Tonankai tsunami combine the 2012 CDMC model and the spectral synthesis method for megathrust subduction earthquakes (Mai & Beroza 2002;Goda et al. 2014). Importantly, our rupture models are not intended for verifying or replacing the 2012 CDMC model, but to constrain our modelling. It is noteworthy that the fault-plane parameters of the 2012 CDMC model (as a whole) are heterogeneous in terms of fault geometry and slip-vector orientation. Both strike and dip angles of individual sub-faults vary (Fig. 1). Because the fault geometry and slip features of the 2012 CDMC model are based on current seismotectonic settings in the Nankai -Tonankai Trough region (e.g. Hirose et al. 2008), the fault-plane information of the 2012 CDMC model is incorporated into our stochastic slip modelling. Specifically, the sub-faults of the 2012 CDMC model are mapped onto a 2D (rectangular) matrix (note: sub-faults for the Kumano splay fault are excluded). The size of the 2D matrix is 54 (downdip) × 153 (along-strike), and its origin is set to the SW corner of the fault plane (i.e. offshore Kyushu: Fig. 1). In our stochastic simulations, slip values for the 54 × 153 matrix are generated, and the synthesized slip values that correspond to the 5669 sub-faults (main fault plane only) are mapped back to the CDMC rupture surface by ignoring slip values outside the main fault plane. The procedure for generating stochastic earthquake slip models for the Nankai-Tonankai tsunami and for simulating tsunami waves is illustrated in Figure 3. The method consists of stochastic synthesis of slip distributions and subsequent temporal rupture specification by considering uncertain rupture start locations. The model parameters for generating stochastic slip distributions include: stochastic synthesis parameters (i.e. correlation lengths in the downdip and along-strike directions, Hurst number for the von Kármán auto-correlation model, and Box-Cox power parameter), target slip parameters (i.e. mean and maximum slips) and kinematic rupture parameters (i.e. rupture nucleation point, rupture propagation velocity and rise time). For future Nankai-Tonankai scenarios, these parameters are uncertain or unknown, and thus are treated as random variables in stochastic slip simulations. In the developed computer code for Monte Carlo tsunami simulation, the probabilistic information of the above-mentioned parameters can be represented by: (i) multiple discrete values (i.e. a set of parameter values and corresponding weights); (ii) uniform random variable (i.e. lower and upper values); or (iii) truncated normal variable (i.e. mean, standard deviation, lower limit and upper limit). In this paper, we adopted the truncated normal distribution to characterize the variability of the source parameters.

Slip distribution modelling
The key slip characteristics for stochastic source modelling are specified in terms of slip statistics, slip distribution parameters and asperity areas. The slip statistics, such as mean and maximum, are used to control the overall features of the resulting slip distribution. The spatial heterogeneity of earthquake slip is modelled in the wavenumber domain by adopting a von Kármán auto-correlation function (Mai & Beroza 2002), but other spectral models could be used. The relevant parameters of the von Kármán function are the correlation lengths (for the downdip and along-strike directions) and the Hurst number. The correlation length determines the absolute level of the power spectrum in the low wavenumber range, and captures the anisotropic spectral features of the slip distribution (when different correlation lengths are specified for downdip and along-strike directions). The Hurst number controls the slope of the power spectral decay in the high wavenumber range, and is theoretically constrained to range between 0 and 1. In addition, realistic tail characteristics of the slip distribution (i.e. slip values often have right-skewed distributions) are modelled via Box-Cox transformation to represent non-normal distribution of earthquake slip (Goda et al. 2014). The asperity areas, which are used for pattern matching for trial synthesized slip distributions, can be specified either based on the 2012 CDMC model (11 rupture cases: Fig. 2) or based on depth limits. For the former, multiple CDMC models can also be considered with different weights. For the latter, the asperity areas are determined by the lower and upper depth limits of the sub-faults (e.g. depth range between 1 and 10 km).
For a given set of stochastic synthesis parameters, slip distributions are generated using a Fourier integral method (Pardo-Iguzquiza & Chica-Olmo 1993). To ensure that synthesized slip distributions have some desirable asperity features (e.g. large asperities are located at shallow depths), multiple trial slip distributions are generated and evaluated by a pattern-matching algorithm that implements selected criteria. In the adopted procedure, two criteria are considered: (i) slip concentrations around the maximum slip are sufficient to form asperity areas; and (ii) the maximum value of the synthesized slip distribution is located within asperity areas of the fault plane. For instance, slip near the main asperity areas (about 10% of the total fault plane) should have more than 30% of the total slip over the fault plane (Mai et al. 2005;Goda et al. 2014).
Subsequently, additional (optional) adjustment of slip values may be needed. A taper function is applied to deeper segments of the fault plane; slip values for sub-faults with top-edge depths of 25 -27, 27-29, 29-31, 31-33 and 33 -35 km are multiplied by the reduction factors of 0.9, 0.7, 0.5, 0.3 and 0.1, respectively (note: the maximum top-edge depth of the 2012 CDMC model is ,35 km, and a diminishing taper is applied in the 2012 CDMC model). The mean slip is adjusted to match the target mean slip (treated as a random variable in our approach). Then unrealistically large slip values, exceeding the target maximum slip (treated as a random variable), are resampled. For resampling of very large slip values, a histogram is constructed based on the synthesized slip distribution with lower and upper limits, for instance, by adopting slip values between two-thirds of the target maximum slip and the target maximum slip (i.e. slip values corresponding to the very large asperities). The simulated slip distributions, after the adjustments, are then used as stochastic slip models for the Nankai -Tonankai tsunami (as instantaneous rupture).

Sampling of model parameters
-Correlation lengths and Hurst number -Box-Cox power parameter -Target mean and maximum slips -Rupture propagation velocity and rise time -Asperity zone Stochastic synthesis of slip distributions -Generation of a random field -Non-linear slip scaling using the Box-Cox --transformation -Pattern matching for asperity characteristics -Tapering of slips for deep segments of the fault -Adjustment for the target mean slip -Adjustment for the target maximum slip

Kinematic rupture modelling
The kinematic rupture process is taken into account following Mai et al. (2005). The probability density functions (PDFs) for hypocentre locations are defined based on the statistical models developed by Mai et al. (2005). The preliminary PDFs for hypocentre locations are specified based on the fault dimensions and the mean/maximum slip ratios. Subsequently, further constraints are placed to exclude unlikely hypocentre locations for a given slip distribution, using empirical findings by Mai et al. (2005). By combining the preliminary PDF and the constraints, the final PDF (or strictly weighting function) for hypocentre locations is obtained. This is used to sample the location of hypocentres. Using the randomly generated rupture propagation velocity and rise time (assumed constant for all sub-faults), the kinematic rupture process of the synthesized slip distributions can be simulated.

Illustration of constrained-stochastic source rupture simulation
The above-mentioned method for stochastic source-rupture modelling for the Nankai-Tonankai tsunami is illustrated in the following. Consider that all model parameters are sampled from assumed probability distributions. For a particular illustration case, the correlation lengths for the downdip and along-strike directions are 60.95 and 104.70 km, respectively; the Hurst number is 0.84; the Box-Cox parameter is 0.24; the target mean and maximum slips are 11.74 and 57.67 m, respectively; the seismic moment and moment magnitude are 6.8 × 10 22 N m and 9.19, respectively (using the constant rock rigidity of 40.9 GPa); the asperity areas are defined based on the depth limits of 1.0 and 12.5 km; the rupture propagation velocity is 2.83 km s 21 ; and the rise time is 60.58 s. Figure 4a shows a synthesized slip distribution prior to the non-linear scaling of slip values based on the Box -Cox transformation. The generated slip values on the fault plane are approximately standard normally distributed (i.e. mean equal to 0 and standard deviation equal to 1). This slip distribution already satisfies the criteria set for pattern matching (thus, it is adopted for further adjustments). The black grid marks the sub-faults, with top-edge depths between 1.0 and 12.5 km (as specified for the asperity areas). The grey square represents the maximum slip over the fault plane, and the grey rectangle is the asperity zone defined for pattern matching (fractional dimensions of the asperity zone are 25 and 40% of the fault width and length, respectively). The threshold for slip concentration is set to 30%. Figure 4b shows the final (instantaneous) slip distribution after non-linear scaling of slip and subsequent adjustments. Owing to the nonlinear transformation, the asperity areas, slip values of which range around 40-55 m, become more prominent. The effects of tapering slip values for deep sub-faults can be seen by comparing slip values near the downdip edges of the fault plane (Fig. 4a, b). The mean-slip value is adjusted to achieve the target mean-slip value of 11.74 m. The standard deviation of the slip values is 12.22 m; the coefficient of variation of the slip values is comparable to those from stochastic models for the 2011 Tohoku earthquake (Goda et al. 2014). The target maximum slip for this case is 57.67 m, whereas several slip values (after the Box-Cox transformation and the meanslip adjustment) exceed this value. To limit the maximum slip value to within a reasonable range, resampling of the slip values for sub-faults with very large slips exceeding the target maximum value is carried out. The effects of resampling slip values for the asperity areas can be seen in Figure 4, which may change somewhat the location of maximum slip. Figure 4c shows the histogram of slip values after the adjustments, exhibiting heavy right-tail features but less than the target maximum value.
Subsequently, a hypocentre location (i.e. rupture starting point) is determined based on the method proposed by Mai et al. (2005). Figure 5 shows: (a) the preliminary PDF for hypocentre locations based on the fault dimensions and the mean/maximum slip ratios; (b) constraints based on large, as well as very large, asperity areas; and (c) the final PDF for hypocentre locations. All values shown in Figure 5 are normalized to 1. The grey square represents the maximum slip of the final simulated slip distribution. The preliminary PDF function indicates that the likely location of hypocentres is near the centre of the fault plane and around sub-faults with relatively large slip. The contour maps for the constraints on the hypocentre locations reflect empirical findings from the statistical analysis of source-inversion models (Mai et al. 2005). The resulting constraint function indicates that the hypocentre is more likely to be close to very large asperity areas (but not too close) and that it is very unlikely to have hypocentres far away from the very large asperity areas. The final PDF combines the preliminary PDF function and the constraint function, and can be used to sample the location of hypocentres. The sampled hypocentre location for the illustrated case is indicated as a grey circle in Figure 5c. Figure 6 displays snapshots of the overall temporal rupture evolution for the simulated slip distribution (between 0 and 180 s), in which the hypocentre location (rupture nucleation) is marked by a grey circle. The slip distributions are generated at an interval of 10 s by considering the rupture propagation velocity of 2.83 km s 21 and the rise time of 60.58 s. The entire rupture process takes 239.10 s (i.e. 178.52 s to reach the farthest sub-fault from the hypocentre plus 60.58 s to complete the subfault rupture process). This information can be used in tsunami simulation to take into account the kinematic process of earthquake rupture.

Monte Carlo tsunami simulation for Nankai-Tonankai earthquake
In this section, we discuss our Monte Carlo tsunami simulation for the Nankai -Tonankai regions. The main purpose of the investigations is to demonstrate that the newly developed method can produce useful results for tsunami hazard mapping and preparedness. In the following, a computational set-up of the tsunami simulation for the Nankai-Tonankai earthquake is described, followed by tsunami simulation results for both offshore and inland areas. A rigorous assessment of uncertainty, as well as a detailed sensitivity analysis of tsunami hazard predictions to model parameters, is beyond the scope of this study, but will be subject of a subsequent investigation.

Tsunami modelling and input data
The tsunami modelling is carried out using a welltested numerical code (Goto et al. 1997) that solves the non-linear shallow-water equations using a leap-frog staggered-grid finite-difference scheme and is capable of generating offshore tsunami propagation and inundation/run-up. The run-up calculation is based on a conventional moving boundary approach, 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 h duration, which is sufficient to model the most critical phases of tsunami waves for the Nankai -Tonankai scenarios. The multi-domain nesting from coarse to fine resolution is conducted to consider large-to small-scale tsunami waves, depending on water depth.
A complete dataset of bathymetry/elevation, coastal/riverside structures (e.g. breakwater and levees) and surface roughness is obtained from the Cabinet Office of the Japanese Government. The data are provided in the form of nested grids (2430 -810 -270 -90 -30 -10 m), covering the geographical regions of the Nankai-Tonankai Tough and western-central Japan (Fig. 1). The lowlying land areas along the coast are covered by 10 m grids for accurate inundation modelling. Because this study focuses on regional-scale effects of variable slip distributions on simulated tsunami waves, the minimum considered grid size for the Monte Carlo tsunami simulation is set to 90 m.
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 raw data are gridded using a triangulated irregular network. The land elevation data are based on the 5 m grid digital elevation model (DEM) developed by the Geospatial Information Authority of Japan. The raw data are based on airborne laser surveys and aerial photographic surveys. These data have measurement errors (in terms of standard deviation) of less than 1.0 m horizontally and of 0.3-0.7 m vertically. The reference elevation of the bathymetry and elevation data is the standard altitude in Japan (Tokyo Peil), and no variation of sea levels is taken into account.
The elevation data of the coastal/riverside structures are primarily provided by municipalities, supplemented also by the national coastline database. In the coastal/riverside structural dataset, only structures with dimensions less than 10 m are represented, noting that those with dimensions greater than 10 m are included in the DEM data. In the tsunami simulation, the coastal/riverside structures are represented by a vertical wall at one or two sides of the computational cells. To evaluate the volume of water that overpasses these walls, the overflowing formulae of Honma (1940) are employed for coastal breakwater modelling at sub-grid scale (Japan Society of Civil Engineers 2002).
In the tsunami simulation, the bottom friction is evaluated using Manning's formula, following the Japan Society of Civil Engineers (2002) standard. Manning's coefficients are assigned to computational cells based on national land use data in Japan: 0.02 m 21/3 s for agricultural land, 0.025 m 21/3 s for ocean/water, 0.03 m 21/3 s for forest vegetation, 0.04 m 21/3 s for low-density residential areas, 0.06 m 21/3 s for moderate-density residential areas and 0.08 m 21/3 s for high-density residential areas. In the dataset provided by the Cabinet Office, no information on the roughness coefficient is provided for grids greater or equal to 90 m. In such cases, a uniform value of 0.025 m 21/3 s (i.e. ocean/water) is adopted, and thus the extent of the inland inundation may be overestimated.

Simulated earthquake slips and slip patterns
The constrained-stochastic source modelling for the Nankai-Tonankai earthquake is applied in generating 100 slip distributions by taking into account variability in the stochastic synthesis parameters, target slip parameters and kinematic rupture parameters. All these parameters are modelled using the truncated normal distribution. The considered values of the parameters are summarized in Table 1. They are determined based on the stochastic source models for the Tohoku earthquake (Goda et al. 2014) and the 2012 CDMC model. The moment magnitudes of the source models are controlled by the mean-slip parameter, which can vary from 6 to 12 m (the average is 9 m). The adopted range of the mean-slip parameter is determined based on the observed variation of this parameter for the 11 inverted source models for the 2011 Tohoku earthquake (Goda et al. 2014). We propagate this variability in our tsunami simulations; therefore, the tsunami simulation results discussed in the following need to be interpreted with this consideration. Note that detailed assessments of the assumed model parameters -and their influence on the tsunami effects -need to be performed in future studies, in which alternative sets of model parameters may also need to be considered.
The 100 simulated source models have various slip (asperity) characteristics. To illustrate the variations of the synthesized slip distributions, six source models are shown in Figure 7. They are selected based on three slip patterns: west, middle and east. The classifications are determined by calculating the ratios of slip amount in the west, middle and east regions to the total slip over the entire fault plane. The boundaries of the west-middle and middle -east regions are set to the western end and eastern end of the Kii Peninsula (Fig. 1), dividing the region of the Nankai-Nonankai Trough into 40, 20 and 40% segments. The synthesized source models are categorized as west/middle/east slip pattern, when slip concentrations in the three regions exceed 60, 30 and 50% of the total slip; otherwise, the source models are left with no classification. This leads to 17 west slip patterns, 35 middle slip patterns, 25 east slip patterns and 23 unclassified patterns that have multiple asperities over two or three regions. Figure 7 presents two source models for the west/middle/east slip pattern. It is noted that these classifications are considered for the sake of presentations of the results only. Figure 7 shows that the location and size of asperities vary significantly for different slip patterns (e.g. Fig. 7b v. Fig. 7d v. Fig. 7f). The maximum slip values also differ significantly (e.g. Fig.  7d v. Fig. 7e). These features are controlled by the stochastic synthesis parameters and the target slip parameters, and have a major influence on the tsunami simulation results (Goda et al. 2014). The developed stochastic source models can generate numerous realistic earthquake scenarios, and thus are useful for exploring the effects of different scenarios on the triggered tsunamis and their impacts on coastal communities. These are investigated in the following.

Offshore tsunami results
For each of the 100 synthesized slip distributions, initial boundary conditions for tsunami modelling (i.e. water surface elevation) are computed using the analytical formulae for elastic dislocation by Okada (1985) together with the equation by Tanioka & Satake (1996). The latter is to take into account the effects of horizontal movements of steep seafloor topography on the vertical water dislocation. Subsequently, the calculated water surface elevations are smoothed using a 9-cell × 9-cell moving average function (same as for the 2012 CDMC model) to avoid very steep initial water surface profiles. The water surface elevation profiles are calculated with an interval of 10 s to account for the temporal evolution of the kinematic earthquake rupture process. In the above-mentioned procedure for calculating the water surface elevation, the hydrodynamic response of seawater is not accounted for. Caution needs to be exercised when modelling tsunamigenic earthquakes with major asperities on a gently dipping fault plane at a shallow depth (applicable to the Nankai -Tonankai earthquake), because large peak vertical deformation is generated along the top edge of the dipping fault (Goda 2015) as a result of the assumption of elastic behaviour of the rock encoded in the Okada (1985) solutions. In reality, strongly non-linear (plastic) deformation occurs, which is typically not accounted for in tsunami simulations: instead, the seabed deformation calculated from the Okada (1985) equations is used directly as an input sea surface condition in the tsunami simulation.
First, offshore tsunami wave characteristics due to the stochastic source models are examined by focusing on the Shikoku and Tokai regions (see Fig. 1). Figures 8 and 9 show the maximum tsunami wave-height (maximum surface elevation -this convention is used throughout) contours for Shikoku and Tokai, respectively, by considering the six source models presented in Figure 7. As expected, slip patterns have a significant impact on the maximum wave heights: if major asperities are near to the locations of interest, larger tsunami wave-amplitudes are generated. By inspecting the source model and the corresponding wave heights, the causal relationship between the source and the impact can be easily understood. The comparison of the results for Figure 8a, b (both are for west slip patterns) also suggests that the local features of the slip patterns have a major influence on the tsunami waves along the coastal line due to the near-shore bathymetry. The complex spatial distribution of tsunami wave heights can be observed in the regions between offshore and coastline (Figs 8a-c & 9b, e, f ). The tsunami wave profiles are significantly affected by refraction due to bathymetry and by overlapping of wave propagation due to complex slip patterns, giving non-uniform incoming tsunami profiles to the land. It highlights the importance of the spatial slip distribution, in conjunction with the local bathymetry and coastline, on the tsunami hazard assessment.
To further illustrate the variability of temporal tsunami wave profiles, two specific sites are selected. The first location is off Kochi in the Shikoku region, and the second one is off Hamamatsu in the Tokai region (see Fig. 7 for the locations on the maps). The water depths at the two locations are about 45 and 60 m for Kochi and Hamamatsu, respectively. The time series of water surface elevation from Monte Carlo tsunami simulation results based on the 100 stochastic source models are shown in Figure 10a, b for Kochi and Hamamatsu, respectively. In the figure panel, individual tsunami wave profiles are shown with light blue lines, whereas the statistics (i.e. median and 10th/90th percentiles) of the tsunami simulation results are shown by thick black (solid/broken) lines to represent major trends of the tsunami wave profiles. It can be observed that the tsunami wave profiles at these locations are highly variable. Although the average wave profiles are not so large, rare tsunami waves (e.g. exceeding 90th percentile waves) can be several times greater than the average tsunami wave levels. For example, at both locations, the maximum tsunami wave heights attained by the severest scenarios are more than twice as large as the 90th percentile wave levels. To further investigate the effects of regional slip patterns on the simulated tsunami profiles, average tsunami wave profiles for the west, middle and east slip patterns are compared with the statistics of the 100 tsunami simulations, and the results are presented in Figure  10c, d for Kochi and Hamamatsu, respectively. As one intuitively expects, tsunami wave profiles are strongly affected by the proximity to the major asperities (i.e. for Kochi, the west slip patterns represent more critical scenarios: while, for Hamamatsu, the east slip patterns have more influence), corroborating the importance of the spatial slip distribution.
Lastly, the overlay of tsunami wave profiles provides useful information on the arrival times of major tsunami waves. For instance, large tsunami waves exceeding 5 m are expected to arrive within 20-30 min at Kochi, whereas such significant waves may arrive earlier (15-20 min) at Hamamatsu (the differences are caused by the proximity to the asperities and the seafloor topography). The difference of 5-10 min has important implications in successfully completing horizontal/vertical evacuations in coastal cities and towns (Federal Emergency Management Agency 2008). Hence, the stochastic tsunami simulations facilitate the quantification and visualization of the tsunami-hazard prediction uncertainty, and are particularly useful for emergency managers who wish to be informed of the range of possible consequences they need to deal with.

Onshore tsunami results
The tsunami wave heights from the source to nearshore regions are significant, as shown in Figures  8 and 9. The spatial variability of the maximum tsunami inundation height along the coastal line in the Shikoku and Tokai regions is shown in Figure 11.   For this analysis, computational cells that are close to zero (i.e. depth between 21.0 and 1.0 m) are defined based on DEM data at a horizontal spacing of approximately 1.0 km. In Figure 11, individual simulation results as well as 10th, 50th, and 90th percentiles of the simulation results at individual locations are shown. Individual results from the simulation exhibit significant variability of the inundation height due to stochastic tsunami scenarios. Tsunami wave heights for average and rare cases are spatially variable, influenced by source, path and topographical effects. Although the mean and extreme tsunami wave heights vary along the coast, the ratios between the 50th and 90th percentiles do not change significantly except for particular areas. This is due to local tsunami amplification effects near-shore (e.g. shape and characteristic water depth of a bay). It is important to know both the average tsunami wave height and range of variations for tsunami disaster reduction.
The stochastic tsunami simulations are also useful for assessing the extent of onshore tsunami inundation by developing stochastic inundation maps. The probabilistic hazard maps can be further utilized in tsunami risk assessment and loss estimation (Goda & Song 2016). To demonstrate the procedure for developing stochastic inundation maps for the Nankai-Tonankai earthquake, results are shown for Kochi and Hamamatsu in Figures 12  and 13, respectively. The land elevation data are displayed in Figures 12a and 13a, respectively. The western section of the low-lying areas in Kochi is protected by hills/high grounds. On the other hand, the topography in Hamamatsu is flat along the coastal line. Although both locations have flat bathymetry offshore, the different land-side topographical features affect the predicted inundation maps, especially under extreme situations.
To develop stochastic inundation maps, a tsunami hazard parameter that is suitable for characterizing the city-level inundation extent needs to be chosen. In this demonstration, inundated areas with depths of more than 1 m are considered. It is noteworthy that the chosen parameter is relevant for tsunami inundation damage to wooden houses; different parameters may be adopted depending on the main features of the building portfolio of interest. Using the tsunami inundation results from stochastic tsunami simulations, city-level inundation values can be evaluated for individual source models and the cumulative distribution function (CDF) of the inundation parameter can be derived (Figs  12b & 13b). The developed CDF of the inundation area above 1 m depth is then used to define probabilistic hazard scenarios. For instance, risk managers may be interested in knowing typical (expected) scenarios and worst-case scenarios, and relevant probability levels may be chosen. In this study, two probability levels are selected for illustration: the 50th percentile and the 90th percentile. The inundation areas at the 50th and 90th percentiles are 7.5 and 29.9 km 2 , respectively, for Kochi (Fig.  12b), and are 8.1 and 64.5 km 2 , respectively, for Hamamatsu (Fig. 13b). It is noteworthy that the inundation areas in Hamamatsu can be significantly greater than those in Kochi (Figs 12b v. Fig. 13b) because the total low-lying areas susceptible to tsunami inundation in Hamamatsu are larger than those in Kochi (no major natural barriers; Fig. 12a v. Fig. 13a). Figure 12c, d show the 50th and 90th percentile tsunami inundation maps for Kochi, respectively. The counterparts for Hamamatsu are shown in Figure 13c, d. The comparison of the 50th and 90th percentile maps clearly indicate spatial expansion of the inundated areas and increased inundation depth along the shoreline. Consequently, the more extreme scenarios cause much significant damage and consequences to coastal cities and towns. The visual inspection of the changes of the inundation profiles from the 50th to the 90th percentile scenarios in Kochi and Hamamatsu clearly indicate that topographical features of the low-lying land have paramount influence on the tsunami inundation result. It is important to emphasize that such dramatic changes of the inundation results may be overlooked when a limited number of tsunami source models are used in tsunami hazard analysis. Using a wide range of stochastic source models, such cases can be captured and appropriate risk mitigation measures can be implemented.
Furthermore, the generated hazard maps can be related to the causative source models, helping define critical hazard scenarios. Figure 12e, f show slip distributions that correspond to the 50th and 90th percentile inundation maps for Kochi; Figure  13e, f are for Hamamatsu. The source models shown in Figures 12 and 13 indicate that the locations of the major asperities move closer to the site of interest in case more critical scenarios are considered (this is expected and similar observations were made for Figs 8 & 9). This information is useful for managing tsunami risk proactively, as local municipalities may be interested in preparing a set of tsunami hazard maps for evacuation purposes. Therefore, such stochastic tsunami scenarios and corresponding inundation maps for local communities are useful planning tools. On the other hand, from a wider perspective, emergency managers should prepare several tsunami scenarios and corresponding inundation maps for the region to assess spatial dependency of tsunami hazard at different locations for a range of critical scenarios.

Conclusions
Tsunami hazard assessment for future mega-thrust subduction earthquakes is challenging, because it involves various kinds of uncertainties and assumptions. The stochastic earthquake source models can be used in Monte Carlo tsunami simulation to quantify the uncertainty of tsunami hazard parameters at regional/local scales. This study developed stochastic source models for the M w 9-class Nankai -Tonankai earthquake based on the inversion-based source models for the 2011 Tohoku earthquake and the 2012 CDMC tsunami source model. The modelling procedure took into account uncertainties of the location of rupture nucleation point, the rupture propagation speed and the rise time to capture the kinematic fault rupture process. The applicability of the developed models was demonstrated by generating stochastic source models with different slip patterns, and by comparing the maximum tsunami wave-height contours and temporal waveheight profiles in the Nankai-Tonankai regions. Moreover, the models were employed to develop stochastic inundation maps in Kochi and Hamamatsu. The numerical examples highlighted that the developed tool is useful for exploring the effects of different scenarios on the triggered tsunamis and their impacts on coastal communities, and produces useful results for tsunami hazard mapping and preparedness. In particular, the stochastic tsunami simulations facilitate the uncertainty quantification and visualization of the tsunami hazard predictions.
The developed stochastic source models for the Nankai-Tonankai earthquake can be improved in several ways. Firstly, a critical review of the assumed model parameters should be conducted as the current models are primarily based on the source characteristics for the 2011 Tohoku-type earthquake, while other subduction earthquakes may also be applicable. For such purposes, a broad range of earthquake source models (e.g. those from the 2004 Indian Ocean event and the 2010 Chile event) should be explored to obtain the suitable sets of stochastic synthesis parameters, target slip parameters and kinematic rupture parameters.
Updating of the scaling relationships for the stochastic synthesis parameters is necessary (Mai & Beroza 2002;Murotani et al. 2013). For such a purpose, new scaling relationships for large tsunamigenic earthquakes by  can be adopted; the models have been developed based on an extensive set of 226 inverted rupture models from the SRCMOD database (Mai & Thingbaijam 2014). Subsequently, rigorous sensitivity analysis needs to be carried out by varying key model parameters and by evaluating the effects of these parameters on tsunami hazard predictions. In particular, the effects of changing source rupture parameters (e.g. rupture speed and rise time) should be investigated more thoroughly in future studies. Another important issue is to develop a formal method to incorporate the geological evidence of the anticipated earthquake fault rupture behaviour, such as segmentation, rupture limit along dip and asperity locations based on past events, in a subduction region of interest. This will enable the construction of more realistic stochastic source models. Such extended models are also useful for testing the effects of adopted assumptions on the tsunami hazard predictions.