Simulation-Based Probabilistic Tsunami Hazard Analysis: Empirical and Robust Hazard Predictions

Probabilistic tsunami hazard analysis (PTHA) is the prerequisite for rigorous risk assessment and thus for decision-making regarding risk mitigation strategies. This paper proposes a new simulation-based methodology for tsunami hazard assessment for a specific site of an engineering project along the coast, or, more broadly, for a wider tsunami-prone region. The methodology incorporates numerous uncertain parameters that are related to geophysical processes by adopting new scaling relationships for tsunamigenic seismic regions. Through the proposed methodology it is possible to obtain either a tsunami hazard curve for a single location, that is the representation of a tsunami intensity measure (such as inundation depth) versus its mean annual rate of occurrence, or tsunami hazard maps, representing the expected tsunami intensity measures within a geographical area, for a specific probability of occurrence in a given time window. In addition to the conventional tsunami hazard curve that is based on an empirical statistical representation of the simulation-based PTHA results, this study presents a robust tsunami hazard curve, which is based on a Bayesian fitting methodology. The robust approach allows a significant reduction of the number of simulations and, therefore, a reduction of the computational effort. Both methods produce a central estimate of the hazard as well as a confidence interval, facilitating the rigorous quantification of the hazard uncertainties.


Introduction
The 2004 Sumatra and the 2011 Tohoku tsunami events revealed the extreme threat posed by catastrophic tsunamis for coastal regions around subduction areas globally. These well-documented events have increased the scientific interest towards the tsunami research, resulting in rapid development of sophisticated computational tools and new methods as well as physical observation networks (e.g., tsunami early warning systems). Accurate predictions of tsunami hazards are essential for enhancing the preparedness and resilience against future tsunami disasters. In particular, rigorous quantification of tsunami hazard uncertainties at various spatial scales is critically important.
According to González et al. (2009), tsunami hazard assessment methodologies can be classified into three broad categories: (a) probabilistic tsunami hazard analysis (PTHA), (b) worst-case scenario approach, and (c) sensitivity analysis. In this work the focus is upon PTHA. PTHA has many common features with probabilistic seismic hazard analysis (PSHA, Cornell 1968), on the other hand there are important differences: (a) tsunami intensity measures are obtained through numerical inundation simulations and not by means of empirical relationships, such as ground motion prediction equations (GMPEs); (b) being based on numerical inundation simulations, PTHA automatically considers the spatial correlation among tsunami hazard estimates at different locations; and (c) in conventional PSHA, scaling relationships relating magnitude to seismic source characteristics are not explicitly considered because GMPEs used in PSHA do not account for such features, whilst details of earthquake rupture process (e.g., slip distribution and rise time) have major influence on tsunami simulations. Limitations in the adoption of GMPEs has been already emphasized also in recent studies proposing physically based PSHA (Convertito et al. 2006;Hutchings et al. 2007) in which the empirical relations are substituted with several wave propagation simulations based on empirical Green's functions.
The existing PTHA can be classified into three categories. In the first category, PTHA is conducted using tsunami catalogs (Burroughs and Tebbens 2005;Kulikov et al. 2005;Tinti et al. 2005;Orfanogiannaki and Papadopoulos 2007), whereas in the second category, different ''scenario-based'' PTHA methods are suggested (Geist and Dmowska 1999;Downes and Stirling 2001;Farreras et al. 2007;Liu et al. 2007;Power et al. 2007;Yanagisawa et al. 2007;Burbidge et al. 2008;González et al. 2009;Løvholt et al. 2012). In the third category, a combination of the previous two is considered (Geist 2005;Geist and Parsons 2006;Annaka et al. 2007;Thio et al. 2007;Burbidge et al. 2008;Parsons and Geist 2009;Grezio et al. 2010Grezio et al. , 2012Horspool et al. 2014;Fukutani et al. 2016). Recently, De Risi and  have developed a new simulation-based multihazards approach for ground motions and tsunamis. The method combines PSHA and PTHA, and computes earthquake and tsunami hazard curves for specific locations, starting from the same source characteristics. The major novelties of the developed method include: (1) slip distribution of earthquake rupture is taken into account in the assessment, unlike conventional uniform slip distribution over a fault plane, (2) uncertainties of earthquake source parameters, i.e., fault width, fault length, mean slip, etc., are modeled using probabilistic prediction models of these parameters, accounting for their variability and dependency, (3) a wide range of magnitude scenarios is considered by characterizing regional seismicity of the target region in terms of occurrence rate of major earthquakes and their relative frequency, and (4) inland inundation of incoming and receding tsunami waves are simulated, rather than stopping at offshore locations, producing more accurate and realistic estimates of tsunami hazard parameters. The adoption of realistic heterogeneous slip distributions in a stochastic simulation framework is a key point for reliable results, as also emphasized in Li et al. (2016). It is important to highlight that the computational requirements for the simulation-based method are relatively high, in comparison with conventional methods. De Risi and  suggested that for the region and site considered in their analyses, about 300 tsunami simulations per magnitude scenario were necessary for reliable estimates of the tsunami hazard. Therefore, it is desirable to be able to reduce the number of simulation runs.
In this paper, the methodology proposed by De Risi and  is firstly adopted and extended to carry out tsunami hazard assessments at regional scale, rather than a single location. This extension is particularly useful for city-level tsunami hazard mapping. Secondly, a new Bayesian statistical method is proposed to improve the statistical robustness of the tsunami hazard prediction. It is noted that De Risi and Goda (2016) employed a classical statistical method to develop a tsunami hazard curve by treating tsunami simulation results as empirical data. Although this is straightforward, it requires a large number of simulations to obtain stable tsunami hazard estimates. The Bayesian method produces so-called robust hazard curves by fitting suitable analytical probabilistic models (e.g., log-normal distribution and Pareto distribution) using less data with respect to the empirical approach. The computational effort of the Bayesian method can be high when fitting involves numerical evaluations of multi-dimensional integrals related to a large number of parameters, which are treated as uncertain variables in the Bayesian method. In the case of tsunami hazard curves, it can be shown that the number of parameters to represent a tsunami hazard curve is relatively small (less than or equal to three), therefore, the Bayesian procedure becomes feasible with a reasonable computational time, which is significantly less than the tsunami simulation runs for numerous source scenarios. Thereby, the adoption of the Bayesian approach allows to reduce significantly the number of simulations required to have reliable tsunami inundation results (up to one-third of the original number). For both empirical and robust methods, three hazard curves are obtained, one corresponding to the central value and two for the confidence interval. In both cases, the confidence interval reflects the limited number of simulations. To demonstrate the developed methodology, the procedure is applied to the Tohoku region of Japan, where the subduction fault plane is well defined and information on regional seismicity is available. Specifically, tsunami hazard for the plain-type coastal region of Miyagi Prefecture, is investigated.

Earthquake Occurrence Model
A classical occurrence model is the memoryless Poisson process (McGuire 2004), generally adopted for long-term hazard assessment, whereas a renewal model (Matthews et al. 2002) may be applied for time-dependent hazard assessment based on recent seismic activity. For the sake of simplicity, in the following, a Poissonian process is assumed. It is noteworthy that the simulationbased framework proposed in this study can incorporate a more general case, such as a renewal occurrence model. Let IM represent the tsunami intensity measure of interest, such as inundation height (h) or flow velocity (v), the probability to observe the first occurrence of a tsunami having intensity value equal to or greater than the specific value im in t years is: where k(IM C im) is the mean annual rate at which the tsunami intensity measure IM will exceed the specific value im, at a given location. In analogy to the PTHA methodology by Parsons and Geist (2009), the rate k(IM C im) can be described as a filtered Poisson process: In Eq. (2), N Sources is the number of subduction seismic sources that are capable of generating a tsunami considered in the analysis. k i (M C M min ) is the mean annual rate of occurrence of the seismic events having magnitudes greater than the minimum magnitude M min for the ith source, whereas f i (M) is the magnitude-frequency distribution characterizing the th siource. The term S(h|M) represents the functional distribution of the uncertain source parameters conditioned on the earthquake magnitude. P i (IM C im|h) is the probability that the tsunami intensity measure IM produced by the ith source will exceed a prescribed value im at a given coastal location for a given set of tsunami source parameters h.
To obtain the terms on the left-hand side of Eq. (2), four phases are defined: (1) definition of input data (i.e., geometrical characteristics of each seismic source and magnitude-frequency distribution), (2) definition of earthquake source scaling relationships and stochastic source model generation, (3) tsunami inundation modeling, (4) statistical analysis of simulated tsunami results and final convolution. The description for each of these phases is presented in the following.

Sources Characterization and Magnitude-Frequency Distribution
The first step is the identification of all seismic sources that are capable of producing damaging tsunami inundation at a site. The sources that are of interest for PTHA are usually located in subduction zones at convergent plate boundaries, and they are known from seismological studies. For instance, detailed geometrical information on subduction zones is available online (Hayes et al. 2012). With respect to the location of interest for which PTHA is performed (e.g., the star in Fig. 1a), tsunamigenic sources can be divided into near-field and far-field (Fig. 1a). Tsunamis triggered by near-field seismic sources can be regarded as the main contributors of the tsunami impact, and they should be studied in detail. In comparison to local tsunamis, a simpler parameterization is usually sufficient for far-field tsunamis because seismic moment, source mechanism, and radiation pattern are more influential than slip distribution within a rupture plane (Geist andParsons 2006, 2016).
The procedure proposed herein simplifies the geometry of the tsunamigenic sources by defining one or more curved surfaces having a rectangular shape (Fig. 1b). Slip distributions of past events from literature can be used to determine the source geometry (Fig. 1b). The seismic source must be able to accommodate the maximum magnitude (M max ) and should be consistent with the magnitude-frequency distribution. Extreme values of magnitude (e.g., larger than 9) should be considered carefully, since such large earthquakes are rare and may span across multiple rupture segments. These events may well be modeled by the characteristic magnitude model (Youngs and Coppersmith 1985).
Since the methodology considered in this study is based on the stochastic synthesis of simulated slip distributions representative of realistic seismic events, a discretization of the fault plane into many sub-faults (Fig. 1c), generally having variable dip, is required. The dimension of the sub-faults must allow  a Subduction zones, near-field and far-field sources. b Simplified source geometry and c source meshing. d Spatial distribution of earthquakes around the world according to the NEIC catalog. Schematic representations of e Gutenberg-Richter relationship and f discrete probability mass based on the fitted Gutenberg-Richter relationship an accurate modeling of the slip distribution corresponding to the minimum magnitude (M min ) that is considered in the magnitude-frequency distribution. This minimum value of magnitude should also take into account that small-to-moderate earthquakes rarely generate significant tsunamis and their contributions to the tsunami hazard are negligible ).To describe the earthquake size in a target region, i.e., the term f(M) in Eq.
(2), a truncated Gutenberg-Richter relationship (Gutenberg and Richter 1956) can be adopted: where the b value is calibrated on the basis of historical events available from earthquake catalogs (e.g., the NEIC earthquake catalog, http://earthquake. usgs.gov/earthquakes/search/, Fig. 1d). Subsequently, the mean annual rate of occurrence of earthquakes with magnitudes greater than or equal to M min falling in that area can be calculated from the fitted Gutenberg-Richter relationship (Fig. 1e). For simulation, it is convenient to convert the continuous distribution of magnitudes into a discrete set of values (M min , …, M i , …, M max ), by adopting a specific discretization interval DM. The discrete probability can be calculated as follows: The probability mass function (pmf, Fig. 1f) for the discrete values of magnitude presented in Eq. (4) is normalized (conditional) with respect to the occurrence rate for the minimum magnitude event. Since the discrete magnitude is considered, the integral in Eq. (2) is replaced by a summation.
It is noteworthy that the adoption of the GR model with a Poisson occurrence process and estimating model parameters based on short earthquake catalogs may not produce the reliable estimate of the longterm recurrence rate for large earthquakes ([M8.5) given the lack of major historical events in modern instrumental catalogs (e.g., 869 Jogan earthquake for the Tohoku case; see Sawai et al. 2012). The extrapolation of the fitted magnitude-recurrence model should be considered carefully (Pisarenki and Rodkin 2010).

Scaling Relationships of Earthquake Source Parameters and Stochastic Source Models
The proposed simulation-based method generates a certain number of stochastic source models to take into account uncertainty related to the rupture process. The simulation is based on the probabilistic models of earthquake source parameters ) and the spectral synthesis method (Goda et al. 2014;Fukutani et al. 2016), characterizing the earthquake slip distribution by wavenumber spectra (Mai and Beroza 2002). Specifically, for each discrete value of magnitude M i (target magnitude hereafter), many samples of the source parameters are necessary to define a slip field on the fault plane comprehensively. Herein, scaling relationships that evaluate the source parameters (e.g., rupture size and spectral characteristics of the slip) as a function of moment magnitude are used for stochastic source generation. Such scaling relationships are obtained on the basis of 226 inverted source models in the SRCMOD database (Mai and Thingbaijam 2014).
Specifically, the two geometrical parameters, i.e., rupture width W and length L, are used to create the rupture area, which is randomly located inside the pre-defined subduction fault plane. Subsequently, a slip distribution is represented as a constrained random field based on desirable seismological features, which are characterized by anisotropic wavenumber spectra (Mai and Beroza 2002), and a realization of such earthquake slip distribution is obtained using a stochastic synthesis method (Goda et al. 2014). The simulation of the random slip distribution is carried out using a Fourier integral method (Pardo-Iguzquiza and Chica-Olmo 1993). The amplitude spectrum of the target slip distribution is specified by a theoretical power spectrum, while the phase spectrum is represented by a random phase matrix. For the amplitude spectrum, the von Kármán model is adopted (Mai and Beroza 2002). According to the von Kármán model, the correlation lengths (CL z along the dip and CL x along the strike) are important source parameters that define the spatial heterogeneity of small wavenumber components in the spectrum. On the other hand, the Hurst number NH determines the spectral decay in the large wavenumber range. All three parameters that describe the slip heterogeneity can be simulated using the scaling relationships by Goda et al. (2016). The obtained complex Fourier coefficients are transformed into the spatial domain via 2-D inverse Fast Fourier Transform. The synthesized slip distribution is then scaled non-linearly to achieve suitable righttail characteristics, in agreement with those observed in the finite-fault models, using the Box-Cox parameter k (Box and Cox 1964) that is modeled as a normal random variable. Finally, the generated slip distribution is further adjusted to have a mean slip (D a ) and maximum slip (D m ), which also scale with respect to the magnitude; the slip parameters can also be simulated using the scaling relationships by Goda et al. (2016).
It is important to note that the error terms of the source parameters W, L, CL z , CL x , D a , and D m mentioned above are considered to follow a multivariate normal distribution , therefore, values of these source parameters can be simulated jointly in the stochastic source generation. Joint random sampling of the source parameters ensures overall consistency in the simulated parameters, leading to a potential reduction of the hazard uncertainties. Nevertheless, due to uncertainty in the source parameters, sampled values of W, L, and D a may result in a seismic moment (i.e., lWLD a ) different from the target magnitude M i . Therefore, to avoid sampled values of W, L, and D a that do not match the target magnitude, consistency of the simulated magnitude with the target is ensured in determining an acceptable source model; in case the calculated moment magnitude does not fall within a certain range, the simulated combination of W, L, and D a is discarded and the sampling is repeated. A tolerance band of ±dM around each magnitude value can be used to define the acceptance criterion in this regard; in this study, dM is set to 0.05.

Tsunami Modeling
As mentioned above, with respect to PSHA, intensity measures are not usually assessed using empirical prediction equations, such as GMPEs. In PTHA, for each stochastic event, the maximum inundation intensity measure for a specific location needs to be computed through numerical simulations of the physical phenomena from the source to the location of interest.
The first step is the calculation of the initial water surface elevation for a given earthquake slip scenario. Specifically, assuming incompressibility of water, initial water surface elevation from the mean sea level can be considered to be equivalent to the seabed displacement field induced by the slip on the fault plane. Such a displacement field can be evaluated using analytical formulae for elastic dislocation (Okada 1985;Tanioka and Satake 1996). To optimize the computation of seafloor dislocation, the seafloor displacement field induced by a unity slip for each sub-fault can be computed in advance, by creating a database of seabed displacement fields. To obtain the total effects of the ith slip distribution, each displacement field is scaled based on the slip in the sub-fault and summed.
Tsunami modeling is then carried out using a suitable numerical code that is capable of generating offshore tsunami propagation and inundation profiles by evaluating non-linear shallow water equations with run-up. See Dutykh et al. (2011) for a summary of available computer codes in literature specific or adaptable to tsunami analysis (e.g., ComCot, FUN-WAVE, MOHID, TIDAL, TUNAMI, and SWAN, among the others). To run tsunami simulations, a comprehensive database of bathymetry/elevation, coastal/riverside structures (e.g., breakwater and levees), and surface roughness is required.

Empirical Tsunami Hazard Curve Based on Tsunami Simulation Results
A conceptual representation of the calculation of the empirical tsunami hazard curve is presented in Fig. 2a. For each seismic source and for each magnitude, simulated intensity measures are used to evaluate the term P(IM C im|M) for the locations of interest (Fig. 2b). Such probability is represented by the complementary cumulative distribution function (CCDF) of the resulting IM. Specifically, IM is represented with the Kaplan-Meier estimator (Kaplan and Meier 1958), being the central estimate: where N sim is the number of simulations. In addition, a confidence interval around the central estimate can be obtained by calculating the variance of the data through Greenwood's formula (Greenwood 1926): The hazard curves obtained in the previous step for each magnitude (i.e., conditional hazard curves) are then multiplied by the probability corresponding to the related magnitude, and eventually are summed up (Fig. 2c). Also in this phase, three curves are obtained, one corresponding to the central value and two for the confidence interval.

Robust Hazard Curves
A disadvantage of the empirical method for developing a conditional tsunami hazard curve for a given magnitude is that a greater number of simulations are required to obtain stable high/low percentile values of tsunami wave height (e.g., 10th and 90th percentiles). In this section, a Bayesian model fitting method is presented. The method replaces the empirical conditional hazard curve with an analytical model using a smaller number of simulation results. Adopting analytical probability distributions has two more advantages. Firstly, parametric models reduce numerical costs with respect to non-parametric models, and improve robustness by avoiding overfitting problems (Zentner 2017). Secondly, for some analytical distributions, it is possible to perform hazard-vulnerability integration in a closed form (Cornell et al. 2002). Such a computationally efficient method is particularly useful, when tsunami inundation simulations are run over land areas represented by high-resolution digital elevation data. Figure 2d shows a procedure adopted for the derivation of the robust hazard curves. In this case, the conditional hazard curves are fitted with an analytical model that is the most suitable for each magnitude (Fig. 2e). The final convolution (Fig. 2f) remains identical to the case of the empirical hazard curve. Specifically, a three-step approach can be followed. Firstly, candidate probabilistic models that are suitable to describe tsunami hazard parameters are selected among those available in literature. Then, an evidencebased Bayesian model selection is carried out to find the most suitable probability distribution. Finally, a robust model is calculated by integrating the analytical distribution of the tsunami hazard and the joint posterior distribution of the hazard curve statistics obtained from the evaluation of the model evidence. Taking advantage of the model selection results, a Bayesian model averaging can also be carried out, which can be compared with the robust model.

Candidate Probabilistic Models
Many studies suggested that the log-normal distribution is the most suitable in fitting tsunami wave heights observed along a given coastal line (Van Dorn 1968;Kajiura 1983;Go 1997;Choi et al. 2002;Kim et al. 2014). On the other hand, other studies indicated that tsunami wave heights can be approximated by different probability distributions (Go et al. 1985;Mazova et al. 1989;Kim et al. 2014). Therefore, further investigations are warranted. In this study, seven distributions, which are widely applied in modeling extreme events, are considered: the Exponential distribution (EXP, Eq. (7)), the Log-Normal distribution (LN, Eq. (8)), the Log-Cauchy distribution (LC, Eq. (9)), the Generalized Pareto distribution (GP, Eq. (10)) with threshold equal to zero, the Generalized Extreme Value distribution (GEV, Eq. (11)), the three-parameter Log-Normal distribution (LN3, Eq. (12)), and the Generalized Logistic distribution (GLO, Eq. (13)).

Mmin Mmax
For each magnitude an empirical CCDF is built EXP is a one-parameter distribution; LN, LC, and GP are two-parameter distributions; and GEV, LN3, and GLO are three-parameter distributions. In this section, the symbol h represents the parameters of the probability distributions and should not be confused with the earthquake source parameters.

Bayesian Model Selection
Conventionally, the best-fit probability distribution for the peak tsunami height can be determined through various types of goodness-of-fit test, such as v 2 test, Cramer-von Mises test, Kolmogorov-Smirnov test, Anderson-Darling test, the probability plot coefficient (PPCC) method, and the L-moment ratio diagram (Kim et al. 2014). In addition, Akaike Information Criterion (AIC, Akaike 1974) and Bayesian Information Criterion (BIC, Schwarz 1978) are also applicable in the model selection. In this study, a Bayesian model selection is considered with respect to other methodologies, Bayes' Theorem, at the model class level, automatically enforces model parsimony without the need for terms penalizing a case with the larger number of uncertain parameters.
Denoting available empirical data to be fitted with the ith model T i among N models by D, the probability of the model T i can be obtained as follows: where p(T i |D) is the posterior distribution of the model T i and can be used to select the most probable one among the considered N models. p(T i ) is the prior distribution of the model T i and can be taken equal to 1/N. Finally, p(D|T i ) is called evidence for the model T i provided by the data D. According to Muto and Beck (2008), the log-evidence can be expressed as the difference of two terms: where p(D|h i ,T i ) is the likelihood function depending on the adopted model T i and its parameters h i . Specifically, assuming that data are independent, the likelihood is the product of the probability of the data D given the parameter h i . p(h i |T i ) and p(h i |D,T i ) are the prior and the posterior of the model parameters, respectively. The prior distribution of the model parameters represents the information available on h i prior to the estimation, while the posterior of the model parameters is obtained according to the Bayesian paradigm (Box and Tiao 1992): It is worth noting that in Eq. (15), the first term of the right-hand side is the posterior mean of the loglikelihood function which gives a measure of average goodness-of-fit of the model T i to the data, and the second term is the Kullback-Libler information (Kullback and Leibler 1951) (or relative entropy), which is a measure of the information gain about T i from the data. Therefore, Eq. (15) explicitly accounts for a trade-off between the data-fit of the model and its model complexity, i.e., how much information it takes from the data (Cheung and Beck 2010).

Bayesian Robust Hazard Curve
The evidence-based assessment facilitates the calculation of robust hazard curves (Jalayer et al. 2015). Once the most suitable model T i is identified, the posterior distribution of its parameters p(h i |D,T i ) can be used to build the robust predictive probability density function of future response X through total probability theorem: The robustness of the model comes from the fact that it takes into account the uncertainty in parameters of the model T i , reflecting the limited number of data D (i.e., simulations). It is worth noting that the robust hazard curve is the expected value of a prescribed probability model taking into account the posterior distribution of the model's parameters. Therefore, the variance of the model can be calculated as follows: The calculated variance can be used to define a confidence interval around the robust model.

Model Averaging
Equation (14) can be used not only for model selection, but also for response prediction based on all models taken into account. Let X denote the quantity to be predicted. According to Total Probability Theorem, the probability of X given the data D can be calculated as follows: instead of using the single best model for prediction as in the robust hazard curve. The operation presented in Eq. (19), which is a weighted average of the models adopting evidence as weight, is also called posterior model averaging in the Bayesian literature (Raftery et al. 1997) or hyper-robust predictive model. In other words, models are weighted according to their degree of belief and based on the quality of the information they deliver.

Case Study
As a case study, the plain-type coastal area of Miyagi Prefecture, in the Tohoku region, Japan, is investigated (Fig. 3). This region has been severely affected by the 2011 Tohoku tsunami. Specifically two sets of locations are focused on in this investigation: 44 points located on the coastal line (elevation about zero) in front of tsunami protection barriers, and 44 points, having the same latitude of the first ones, located inland behind tsunami protection barriers. The mean average height of these protection barriers is about 5.30 m, with a minimum and maximum height of 1 and 10 m, respectively.
To reduce the computational efforts and to focus on the methodological aspect, only a specific seismogenic context, i.e., near-field source in the Tohoku region of Japan (Fig. 1a), is taken into account. Nevertheless, the procedure can be extended to consider all possible sources of interest for the Tohoku region and can be applied to other subduction zones. Furthermore, only geophysical uncertainty is considered herein. Specifically, a 2011 Tohoku-type fault is analyzed with a source zone of 650 km along the strike and 250 km along the dip (Fig. 3); this is an extended fault plane of the source model for the 2011 Tohoku earthquake (Satake et al. 2013). The fault model can accommodate a M9 earthquake, consistent with the maximum magnitude adopted for the magnitude-frequency distribution. Values of magnitude larger than 9 are neglected since simultaneous rupture of the off-the-Tohoku subduction segment and the off-the-Hokkaido subduction segment are not considered. For the stochastic synthesis of simulated seismic events, a 10-km mesh with variable dip is generated. Such discretization allows accurate modeling of the slip distribution that corresponds to a M7.5 seismic event, involving at least 5 by 5 subfaults.
With respect to the magnitude frequency distribution for the Tohoku region, the Japanese Headquarters for Earthquake Research Promotion proposes a hazard model based on the Poisson process, combined with the Gutenberg-Richter magnitude-frequency model. Specifically, a b value equal to 0.9 is adopted (Headquarters for Earthquake Research Promotion 2013) in Eq. (3). To obtain the pmf presented in Eq. (4), seven discrete values of magnitude are considered. Adopting a discretization interval of 0.25 and considering 7.5 and 9.0 as the smallest and largest central discrete values of moment magnitude, the minimum and maximum values of magnitude to consider in the truncated Gutenberg-Richter relationship are 7.375 and 9.125, respectively. Therefore, the seven central magnitude values are 7. 5, 7.75, 8.0, 8.25, 8.5, 8.75, and 9.0. According to the NEIC earthquake catalog (Fig. 1b), the mean annual rate of occurrence of earthquakes with magnitude greater than or equal to 7.375 in the source area is calculated. In this study, the events reported in the database that fall in the considered major rupture area, recorded in the period 1976-2012, having a depth varying between 0 and 60 km, and considering a magnitude range between 5 and 9, are considered. According to the Gutenberg-Richter fitting (Fig. 4a), the rate k(M C 7.375) is estimated to be 0.183. Figure 4b shows the pmf for the discrete values of magnitude. Note that the Gutenberg-Richter model presented herein shown in Fig. 4a is similar to the magnitude-recurrence model adopted by the Headquarters for Earthquake Research Promotion (2013). Figure 5 shows the adopted scaling relationships for four source parameters; the central estimates and the 16th/84th percentiles are presented with continuous and dashed lines, respectively. Simulated data from the stochastic source modeling (grey dots) and their statistics (colored circles) are also shown. Magnitude values for simulated data are not perfectly aligned at the seven discrete values; a tolerance band of ±0.05 around each magnitude value is allowed (see the detail shown in the red dashed circle in Fig. 5). Examples of the slip distributions generated for the M7.5 and M9.0 scenarios are shown in Fig. 6a, c.

Stochastic Source Models and Tsunami Simulations
Once the slip distribution is calculated, tsunami simulations can be carried out. Tsunami modeling is performed using a well-tested numerical code (Goto et al. 1997 Latitude ( tsunami propagation and inundation profiles by evaluating non-linear shallow water equations with run-up using a leapfrog staggered-grid finite difference scheme. The run-up calculation is based on a moving boundary approach, where a dry/wet condition of a computational cell is determined based on total water depth relative to its elevation. To catch the most critical phases of the tsunami waves, the numerical tsunami simulation is performed for a time window of 2 h following the seismic event. The Courant-Friedrichs-Lewy condition determines the integration time step DT; it is a function of bathymetry/elevation data and their resolution. In this study, DT is equal to 0.5 s. Results from simulations are the maximum tsunami intensity measures of interest (i.e., wave height, flow velocity, etc.) for one or more specific locations along the coast; aggregate tsunami hazard parameters, such as inundation areas above a certain depth, can also be evaluated.
The Miyagi prefectural government provided a complete dataset containing information about (a) bathymetry and elevation, (b) tsunami defense coastal/riverside structure (e.g., breakwater and levees), and (c) surface roughness. Data are organized as nested grids (1350-450-150-50-m) for the entire geographical regions of Tohoku; therefore, the digital elevation model used for the inundation provided by Miyagi Prefecture has a resolution of 50-m.
More specifically, the ocean-floor topography resolution is 1:50,000; it is based on the bathymetric charts and JTOPO30 database developed by the Japan Hydrographic Association on the basis of nautical charts developed by the Japan Coastal Guard. Municipalities in Miyagi Prefecture provided detailed data about the elevation of the coastal/riverside tsunami protection structures. Such coastal/riverside structures are represented as vertical walls along one or two sides of each interested computational cell. Homma's overflowing formulae are adopted for evaluating the volume of water that overpasses the walls. No tidal fluctuation is considered in this study. Ocean bottom friction is estimated through Manning's formula. Japanese land use data are adopted to assign Manning's coefficients to computational cells; the following assumptions are made: 0.02 m -1/3 s for agricultural land, 0.025 m -1/3 s for ocean/water, 0.03 m -1/3 s for forest vegetation, 0.04 m -1/3 s for low-density residential areas, 0.06 m -1/3 s for moderate-density residential areas, and 0.08 m -1/3 s for high-density residential areas. Figure 6b, d show the results in terms of maximum tsunami wave height for the case study area, corresponding to the two slip simulations for M7.5 (Fig. 6a) and M9 (Fig. 6c), respectively. Results show that tsunami intensity increase significantly with magnitude, both in terms of inundation height and inundated areas.

Effects of Number of Simulations
It is well known that short or incomplete data can lead to biased estimation of the hazard parameters when conventional statistical methods are adopted (Lamarre et al. 1992). In this section, a bootstrap procedure is carried out to study the effect of the number of simulations on the final hazard estimation. Bootstrap provides, through a Monte Carlo simulation, the sub-sampling of a pool of m independent and identically distributed random variables from an initial sample of n elements (with m B n), having a distribution function identical to the empirical distribution function of the original sample. For each generated sample containing m elements, mean, median, and different percentiles can be computed. Such estimates can be used to quantify the uncertainty associated with parameter values.
All points along the coast (Fig. 6) are investigated; for illustration, results for points 8, 20, and 35, located on the coast in front of the barriers, are shown in Fig. 7. For each location, five percentiles (i.e., 5th, 25th, 50th, 75th, and 95th) of the wave height are presented as a function of the number of simulations  Fig. 7). In particular, 300 simulations are necessary for M7.5, 250 simulations for M8.0, and 200 simulations for M8.5 and M9.0. Such a decreasing trend with the magnitude is consistent with the physical process: when the magnitude is relatively small, the rupture area can move more freely over the fault plane (see Fig. 6), increasing the variability of the inundation intensity measures. In turn, when the magnitude is large, the fluctuation of the rupture area is more constrained (i.e., the major slip area tends to occupy the entire subduction plane). This concept is still valid for points located behind tsunami protection barriers; in the latter case the minimum values of magnitude are shifted to values for which the barrier starts to be inefficient (e.g. [8).

Empirical Tsunami Hazard Curves and Hazard Maps
For each value of seven magnitudes (i.e., 7.5, 7.75, 8.0, 8.25, 8.5, 8.75, and 9.0), according to the bootstrap analysis, 300 sets of the tsunami source parameters h are generated using the scaling relationships by Goda et al. (2016) and 300 tsunami simulations are performed. In total, 2100 stochastic simulations are performed in the specific case study.
The simulated results are then treated according to the empirical method presented in Sect. 2.5. The empirical CCDFs in terms of tsunami wave height for point 20 in front of and behind the tsunami protection barriers are presented in Fig. 8a, d, respectively, for seven magnitude values analyzed. Figure 8b, e show the conditional hazard curves, weighted by the probability values obtained from the discretized Gutenberg-Richter relationship (Fig. 4b). The summation of the curves presented in Fig. 8b, d, multiplied by k(M C 7.375) = 0.183 (Fig. 4a), leads to the final empirical hazard curves (Fig. 8c, f).
It is possible to observe that for the point in front of the protection barriers, the slope of the final hazard curve for wave height greater than 5 m is very steep. This is the direct consequence of less variability of tsunami inundation for large values of magnitude. Moreover, the tsunami height cannot be so high in the Sendai plain areas unlike ria-type coastal areas (e.g., Onagawa and Kesennuma), where the wave amplification due to topographical effects is significant.
The comparison of the hazard curves for the two points with the same latitude in front of and behind the protection barriers, having different elevations, shows another important result: when the location is further inland or at higher elevation and when the location is protected by tsunami defense, contributions of medium values of magnitude become less important (values between 7.5 and 8.25 in the examined case, Fig. 8a, d). Therefore, the maximum mean annual rate of occurrence decreases, and the extension of the flat part of the hazard curve increases progressively, going from the coast to inland (Fig. 8c, f) since only higher values of magnitude produce significant effects at inland locations. A consequence of this result is that for inland locations a smaller number of simulations can be adopted by ignoring smaller magnitude cases. Figure 9a, c show the hazard curves in terms of tsunami wave height for all points along the coast, for points in front of and behind the tsunami protection barriers, respectively. The results for points 8, 20, and 35 (same locations considered in Fig. 7) are highlighted. Figure 9b, d show the hazard curves for the same points in terms of tsunami depth that is the tsunami wave height corrected for the local topography/elevation; such correction partially removes the flat part of the curves for the inland locations. No major differences can be noted between Fig. 9a, b, since the considered locations have elevation approximatively equal to zero. Some of the hazard curves shown in Fig. 9a present a flat part that is lower than the maximum mean annual rate (i.e., 0.183); such a result is observed especially for coastal locations close to harbor protection structures (i.e., breakwater, reef, etc.). Structures protecting harbors from storms have demonstrated to be effective in protecting also from tsunamis generated by small magnitude seismic events.
The comparisons of the hazard curves for points 8, 20, and 35 on the coastline and inland (Fig. 10) lead to an interesting observation. Tsunami hazard curves for points 8 and 20, located behind the protection barriers, have little contributions from lower values of magnitude (as shown in Fig. 8), therefore, hazard curves in terms of both tsunami wave height and inundation depth are lower than hazard curves for the corresponding points located in front of the protection barriers. Conversely, the hazard curve in terms of tsunami wave height for point 38, located behind the tsunami barriers (dotted brown line in Fig. 10a), is greater than the hazard curve in front of the tsunami barriers (continuous brown line in Fig. 10a). This result is due to the insufficient effectiveness of protection barriers; in fact, in this case, the low values of magnitude still contribute to the final hazard, but only with high values of wave height, specifically higher than local elevation. Generally, the probability of simulated tsunami wave height lower than the local topography is equal to zero, by definition of tsunami wave height. This effect is obviously lost for curve in terms of tsunami inundation depth (Fig. 10b). This observation highlights the need of adopting depth-based tsunami fragility curves for the risk assessment, especially for the comparison of different mitigation strategies.
An alternative hazard representation is shown in Fig. 11; specifically, tsunami wave height (Fig. 11b, e) and tsunami inundation depth (Fig. 11c, f) Fig. 11d-f show results behind the protection barriers. It is possible to observe that for some points the depth values are different from zero only for high values of return period. This is related to the presence of protection structures and/or to the topographical effects. Finally, the hazard computation procedure can be extended to obtain uniform tsunami hazard maps. The term 'uniform' refers to the same annual probability of exceedance of the tsunami intensity values for all points in a given geographical region. Figure 12 shows the uniform tsunami hazard maps for five values of return period: 30, 50, 475, 975, and 2475 years. These maps are obtained by repeating the simulations shown for a single point, for a set of grid points covering the area of interest. As tsunami intensity measure, maps show both tsunami wave height (Fig. 12a-e) and inundation depth (Fig. 12f-j). It is possible to observe that values of tsunami

Robust Tsunami Hazard Curves
In Sect. 4.2, it has been indicated that 300 simulations are sufficient to obtain stable results for several locations along the analyzed coast in terms of high and low percentile values of tsunami wave height (e.g., 10th and 90th percentiles). In the following, only 100 simulations are used for the construction of the robust hazard prediction through the Bayesian procedure, and the results are then compared with the empirical results based on 300 simulations. This demonstrates the possibility to reduce the number of simulations by adopting analytical models. For illustration, results for the point 20 located in front of the protection barriers are obtained and shown below.

Bayesian Model Selection
The first step for the model selection is the computation of the posterior joint distribution p(h i |D,T i ) of the parameters h i for each considered model. This is calculated according to Eq. (16). The adopted data are the tsunami wave height results from 100 stochastic tsunami simulations, chosen randomly from the 300 simulations presented before. For the prior distribution p(h i |T i ), in absence of other information, independence of the parameters is assumed, therefore, the joint prior distribution is decomposed into marginal distributions of independent variables. Figure 13 shows, for example, the posterior distributions of the three parameters that are necessary to define the GEV model (Eq. (11)) for the seven magnitude values. Once the posterior distributions for all parameters of all selected models are obtained, it is possible to calculate the evidence according to Eq. (15) and then the probability of the models according to Eq. (14). Figure 14 shows the probabilities of the seven considered probability models. It can be observed that the GEV model is preferable in all cases. It is noteworthy that the LN model and the three-parameter GLO distributions are equally competitive. Table 1 shows the maximum likelihood parameters obtained from the posteriors corresponding to the best models that are based on the Bayesian model selection.

Robust Hazard Curves and Model Averaging
Once the best model for each magnitude is identified and the posterior distributions of the parameters are determined, it is possible to build the robust hazard curve according to Eqs. (17) and (18). Figure 15a shows the final robust hazard curve (central estimate and confidence interval). Results indicate that the hazard confidence interval (dashed red lines) around the central estimate (solid red line) contains the hazard results obtained from the 300 simulations (the jagged blue line). It is worth noting that there are some discrepancies between the robust hazard curve and the curve representing the empirical data; specifically, the robust hazard curve slightly overestimates the mean annual rate of occurrence at low annual rate levels. This is mainly due to the probability contribution of the right tails of the analytical distributions associated with the mediumhigh values of magnitude. Finally, model averaging, according to Eq. (19), is carried out, and the results are shown in Fig. 15b, presenting the hyper-robust curves obtained by averaging the different probability models according to their evidence values. Also in this case, the confidence interval associated with the final hazard curve contains the empirical results based on the 300 simulations. However, no significant fitting improvement can be observed with respect to the simple robust model shown in Fig. 15a because for the examined case GEV presents the highest value of evidence and, therefore, the contribution of the other models is not significant.

Conclusions
Classical and robust probabilistic tsunami hazard assessment procedures, based on the new PTHA procedure proposed by De Risi and Goda (2016), have been investigated and compared. Similarities and differences with respect to PSHA were emphasized. New global scaling relationships of earthquake source parameters for tsunamigenic events were used to generate a wide range of earthquake scenarios corresponding to pre-defined discrete magnitude values. For each scenario, an inundation simulation is carried out. Eventually, the inundation results are  h -tsunami wave height (m)

Figure 15
Final a robust and b hyper-robust tsunami hazard curves combined to obtain tsunami hazard curves and tsunami hazard maps. A study on the sufficient number of simulations required to obtain stable high and low percentiles of the inundation results was presented.
Two types of hazard curves were presented. The empirical hazard curve was built according to traditional statistical procedures. The robust hazard curve was obtained using an advanced Bayesian methodology that allowed reducing the number of simulations by adopting analytical probabilistic models. Both hazard curves can produce a confidence interval, allowing the potential propagation of the hazard uncertainties in the risk assessment.
The procedure was applied to the plain-type coast of the Tohoku region (Japan), and 44 points, which are located on the coastline between the cities of Shinchi and Sendai, were considered for assessing the tsunami hazard. 300 simulations were performed to obtain the empirical hazard curve. The number of simulations can be reduced to 100 when the robust Bayesian fitting is adopted. From the Bayesian model fitting, it can be concluded that for the considered locations, the Generalized Extreme Value distribution is a suitable model for tsunami inundation height and depth for all investigated magnitude values. A potential further development of the procedure proposed in this study is the improvement of the Bayesian integration by the implementation of advanced integration strategies, such as Markov Chain Monte Carlo simulations.
The hazard curves calculated for a lattice of points facilitated the generation of uniform tsunami hazard maps corresponding to different probabilities of occurrence in a given time window. Such maps are useful for urban planning and for tsunami evacuation. The observation of the hazards for all points along the coast facilitated the understanding of the efficacy of protection barriers against medium-size tsunamis.