The formulation presented herein is aimed at developing the earthquake–tsunami hazard curves for a specific location. Let IM represent the intensity measures of interest, such as spectral acceleration (S_a), inundation height (h), flow velocity (v), flux momentum, and tsunami force. Assuming a Poissonian arrival time process, the probability to observe an earthquake–tsunami sequence having intensity measure values IM equal to or greater than the specific values im in t years is (1)P(IM≥im|t)=1−exp[−λ(IM≥im)⋅t] where λ(IM ≥ im) is the mean annual rate at which the intensity measures IM will exceed specific values im at a given location. The rate λ(IM ≥ im) can be expressed as a filtered Poisson process [e.g., Parsons and Geist (2008)]: (2)λIM≥im=λM≥Mmin⋅∫ PIM≥im|θ⋅ Sθ|M⋅fM⋅dM λ(M ≥ M_min) is the mean annual rate of occurrence of the seismic events with magnitudes greater than the minimum magnitude considered in the magnitude–frequency distribution. P(IM ≥ im|θ) is the probability that the joint intensity measures IM will exceed prescribed values im at a given coastal location for a given set of source parameters θ. S(θ|M) represents the scaling relationships (or prediction models) of the uncertain earthquake source parameters conditioned on the magnitude. f (M) is the magnitude–frequency distribution.

Five phases are defined: Phase 1 – fault model and earthquake occurrence, Phase 2 – source parameter characterization and stochastic slip synthesis, Phase 3 – earthquake simulation, Phase 4 – tsunami simulation, and Phase 5 – development of earthquake–tsunami hazard curves. Detailed descriptions for each of these phases are presented in the following. Figure 2 shows the computational framework of the methodology.

The first step is the identification of all seismic sources capable of producing damaging ground motions and tsunami inundation at a site. In this study, a curved surface is considered. Specifically, a 2011 Tohoku-type fault is analyzed with a source zone of 650 km along the strike and 250 km along the dip (Figure 3A); such geometry is capable of accommodating a M9 earthquake that is consistent with the maximum magnitude adopted for the magnitude–frequency distribution. The fault plane geometry is the extended version of the source model by Satake et al. (2013). Note that extremely large earthquakes that span across multiple seismotectonic segments are not considered (e.g., simultaneous rupture of the off-the-Tohoku subduction region and the off-the-Hokkaido subduction zone). To implement the stochastic synthesis method of earthquake slip distribution, the fault plane is discretized into many sub-faults; a 10-km mesh with variable dip based on Satake et al. (2013) is generated. Such a discretization allows simulating accurately the slip distribution corresponding to a seismic event with M7.5 (i.e., the smallest central magnitude value considered for the magnitude–frequency distribution, as shown later), involving at least 5-by-5 sub-faults.

To describe the earthquake sizes in the target region, i.e., the term f (M) in Eq. 2, a truncated Gutenberg–Richter relationship (Gutenberg and Richter, 1956) is adopted, and its CCDF is given by (3)G(M)=1−10−bM−Mmin1−10−bMmax−MminMmin<M<Mmax where M_min and M_max are the minimum and maximum moment magnitudes, respectively. For the simulation, it is convenient to convert the continuous distribution of magnitudes into a discrete set of values (M_min, …, M_i, …, M_max), assuming that they are the only possible magnitudes; such probabilities are computed as follows: (4)PMi=GMi+0.5⋅ΔM−GMi−0.5⋅ΔM where ΔM is the discretization interval. The discrete term presented in Eq. 4 is used in Eq. 2 instead of f (M).

For the analyses, M_min and M_max are set to 7.375 and 9.125, and a discretization interval of 0.25 is adopted. This means that seven central magnitude values, i.e., 7.5, 7.75, 8.0, 8.25, 8.5, 8.75, and 9.0, are considered to calculate the corresponding conditional probabilities as in Eq. 4. The minimum magnitude value is chosen, since small-to-moderate earthquakes rarely generate significant tsunamis, and their contributions to the tsunami hazard are negligible (Annaka et al., 2007). For the Tohoku case study, a b-value equal to 0.9 is adopted (Headquarters for Earthquake Research Promotion, 2013).

Once the magnitude interval is selected and the major source area containing all possible rupture scenarios is defined, the mean annual rate of occurrence of earthquakes with magnitudes greater than or equal to 7.375 falling in that area, i.e., the term λ(M ≥ M_min) in Eq. 2, can be calculated. In order to perform such a calculation, the NEIC earthquake catalog (http://earthquake.usgs.gov/earthquakes/search/) is used. Figure 3B shows the events reported in the database that fall in the considered rupture area, recorded in the period between 1976 and 2012, having a depth varying between 0 km and 60 km, and considering a magnitude range between 5 and 9. According to the data analysis (Figure 3C), the estimated rate λ(M ≥ 7.375) is equal to 0.183. Figure 3D shows the occurrence probabilities for the discrete set of magnitude values (i.e., 7.5, 7.75, 8.0, 8.25, 8.5, 8.75, and 9.0). Note that the probability mass function shown in Figure 3D is normalized (conditional) with respect to the occurrence rate for the minimum magnitude event.

To take into account uncertainties related to the rupture process, multiple random slip fields are simulated (Figure 2A). The simulation procedure is based on a spectral synthesis method (Goda et al., 2014; Fukutani et al., 2016), where the earthquake slip distribution is characterized by wavenumber spectra (Mai and Beroza, 2002; Lavallée et al., 2006). Scaling relationships that evaluate the source parameters as a function of moment magnitude are needed for stochastic tsunami simulation (e.g., rupture size and spectral characteristics of the rupture). In this study, new global scaling relationships for tsunamigenic earthquakes are employed. These relationships are obtained on the basis of 226 inverted source models in the SRCMOD database (Mai and Thingbaijam, 2014). The details of the adopted scaling laws can be found in Goda et al. (2016).

The following relationships are employed to obtain the rupture width (W) and length (L): (5)log10W=−0.4877+0.3125M+0.1464×εW (6)log10L=−1.5021+0.4669M+0.1717×εL where the numbers multiplying by the ε terms are the SDs of the regression errors. The two geometrical dimensions are used to create the rupture area, which is randomly located inside the predefined subduction fault plane. Subsequently, a slip distribution realization with desired properties is obtained using a stochastic synthesis method (Goda et al., 2014). First, a random field, having quasi-normal distribution with a desired spatial correlation structure, is generated 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 considered (Mai and Beroza, 2002): (7)P(k)∝CLzCLx(1+k2)HN+1 where k is the wavenumber (i.e., reciprocal of the wavelength). 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 and are determined from the following scaling relationships: (8)log10CLz=−1.0644+0.3093M+0.1592×εCLz (9)log10CLx=−1.9844+0.4520M+0.2204×εCLx

On the other hand, the Hurst number NH determines the spectral decay in the large wavenumber range and can be modeled as a bimodal random variable that takes a value of 0.99 with probability of 0.43 or a value sampled from the normal distribution with mean equal to 0.714 and SD equal to 0.172 with probability of 0.57 (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 right-tail characteristics, in agreement with those observed in the finite-fault models, using the Box–Cox parameter λ (Box and Cox, 1964). It can be modeled as a normal random variable with mean equal to 0.312 and SD equal to 0.278.

Finally, the generated slip distribution is further adjusted in order to have a mean slip (D_a) and maximum slip (D_m), according to the values calculated from the scaling relationships for the given magnitude: (10)log10Da=−5.7933+0.7420M+0.2502×εDa (11)log10Dm=−4.5761+0.6681M+0.2249×εDm

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 distributed according to a multivariate normal distribution (Goda et al., 2016). The linear correlation matrix of the regression errors ε is given by (12)                  εW            εL           εCLz           εCLx           εDa           εDm ρε=εWεLεCLzεCLxεDaεDm[ 1.0000.1390.8260.035−0.680−0.545 1.0000.2490.734−0.595−0.516  1.0000.288−0.620−0.564   1.000−0.374−0.337    1.0000.835     1.000 ]

Therefore, values of W, L, CL_z, CL_x, D_a, and D_m can be simulated jointly in the stochastic source simulation. The central estimates and the confidence interval (16th and 84th percentiles) of the scaling relationships are shown in Figure 4. The same figure also shows simulated data (green dots) and associated statistics (colored circles), which are obtained from the stochastic source modeling. Magnitude values for simulated data are not perfectly aligned at the seven discrete values; in fact, the simulation algorithm allows a tolerance band of ±0.05 around each magnitude value.

The preceding procedure of earthquake source characterization is innovative with respect to the literature. In particular, a common physical process for concurrent earthquake and tsunami threats is considered, i.e., the common fault rupture scenario that is modeled through the generic stochastic slip scenario on the subduction fault plane. The adoption of a common physical process facilitates the probabilistic investigation of the dependency between shaking and tsunami hazard parameters.

Ground motion prediction equations are extensively used as an effective way to predict seismic intensity measures for a given earthquake scenario (Wald et al., 2006). To account for seismic intensities at multiple locations that occur simultaneously for a given event, GMPEs together with spatial correlations in the regression residuals can be treated as statistical prediction models (Goda and Atkinson, 2010; Goda, 2011). This feature is particularly important in extending the seismic hazard assessment into a risk assessment of a portfolio of buildings/infrastructures. In this study, only the intra-event SD is propagated through the simulation procedure; such a choice is consistent with the simulation scenario of a single fault plane.

Two GMPEs that are applicable to subduction zones are used for the seismic simulations (Figure 2B). The first GMPE (Abrahamson et al., 2016) was developed with a global dataset of earthquakes in subduction zones and has been modified by adding the 2010 Maule Chile and 2011 Tohoku Japan earthquakes to the initial database. The basic functional form of the model for interface subduction events is (13)ln(Sa)=θ1+θ4⋅ΔC1+θ2+θ3⋅(M−7.8)⋅lnR+C4⋅expθ9⋅(M−6)+θ6⋅R+fMAG(M)+fFABA(R)+fSITEPGA1000,VS30+σ⋅ε where ln is the natural logarithm, R is the closest distance to the rupture area, V _S30 is the shear wave velocity in the uppermost 30 m of soil column, PGA₁₀₀₀ is the median peak ground acceleration (PGA) value corresponding to V _S30 = 1000 m/s, σ is the total SD, and ε is the Gaussian error term, represented by 0 mean and unit SD. The SD is period-dependent; it is obtained by the combination of intra-event (ϕ) and inter-event (τ) SDs. The magnitude function is (14)fMAGM=θ4⋅M−7.8+ΔC1+θ13⋅10−M2forM≤7.8+ΔC1θ5⋅M−7.8+ΔC1+θ13⋅10−M2forM>7.8+ΔC1 where ΔC₁ is the term representing the epistemic uncertainty in the break of the magnitude scaling and allows adjusting the GMPE for large interface events that were not originally considered in the earthquake database. f _FABA(R) represents the forearc/backarc scaling term; it is equal to 0 for forearc or unknown site, and this is applicable to the case of this study. Finally, the model for site response scaling is given by (15)fSITE=θ12⋅lnminVS30,1000Vlin−b⋅lnPGA1000+c+b⋅lnPGA1000 + c⋅minVS30,1000VlinnVS30 < VlinfSITE=θ12⋅lnminVS30,1000Vlin+b⋅n⋅lnminVS30,1000VlinVS30≥Vlin

All the model coefficients for Eqs 13–15 can be found in Abrahamson et al. (2016).

The second GMPE by Morikawa and Fujiwara (2013) is suitable for M9 earthquakes in Japan. Two formulations are proposed: one is expressed with a quadratic magnitude term, while the other considers a linear magnitude term. In this study, the quadratic formulation is used since the correction factors (presented in the following) that are included in the quadratic formulation reduce the regression SD. The functional form for interface events is (16)logSa=a1⋅min(M,8.2)−M2+b1⋅R+c1−logR+d1⋅10e1⋅min(M,8.2)+σ⋅ε where log is the base-10 logarithm, a₁, b₁, c₁, d₁, e₁, and σ are period-dependent regression coefficients and can be found in Morikawa and Fujiwara (2013). In Morikawa and Fujiwara (2013), no distinction is made between intra-event and inter-event SDs. The two SDs ϕ and τ can be determined by splitting the total variance into the intra-event and inter-event components based on the ratios of intra-event and inter-event variances to the total variance presented in Zhao et al. (2006). The estimate of the prediction equation is further modified using three additional correction terms: amplification due to the deep sedimentary layers (G_d), amplification due to shallow soft soils (G_s), and anomalous seismic intensity distribution due to the position of the site of interest with respect to the volcanic front (AI). In this study, only the second correction term is taken into account and is given by (17)logGs=ps⋅logminVSmax,VS30350 where VSmax is a period-dependent regression parameter that can be found in Morikawa and Fujiwara (2013).

Figures 5A–D compare the two GMPEs for three seismic intensity parameters [i.e., PGA, S_a(T = 0.3 s), and S_a(T = 3 s)], for M equal to 7.5 and 9.0 and for V _S30 equal to 300 m/s. Figure 6 shows the acceleration response spectra obtained using the two GMPEs, considering two values of closest distance (i.e., 50 and 200 km), and the same values of magnitude and V _S30 described before. Significant differences between the two GMPEs can be observed, especially for the large value of magnitude. Moreover, the Morikawa–Fuijwara GMPE tends to attenuate PGA and short-period spectral acceleration faster than the Abrahamson et al. GMPE with the distance, while the opposite trend occurs for the long-period spectral accelerations.

For seismic simulations, three main inputs are required: event magnitude, distance from the rupture, and shear wave velocity for the considered site. Regarding the distance from the rupture, the GMPEs presented above are both based on the closest distance between the location of interest and the rupture area (Figure 7A). To optimize the computation of the shortest distance, the distances between the coastline location and each discretized element of the 2011 Tohoku-type fault are precomputed and stored (Figure 7B). As an example, Figure 7C shows the distances computed for 500 stochastic scenarios, which are considered in Section “Results.” It is worth noting that the minimum distance is circa 50 km, corresponding to the depth of the fault plane under the considered location. As observed in Goda and Atkinson (2014), the source-to-site distance is affected by the location and size of the fault plane, which in turn is determined by the magnitude of the event. In Figure 7C, it can be observed that the greater magnitude value results in smaller variability of the closest distance (i.e., distribution function has a steeper slope). This is because the rupture plane can move more freely within the overall fault plane (Figure 2A) when the earthquake magnitude is small. For the shear wave velocity, the USGS global V _S30 map server is used (Wald and Allen, 2007).

Finally, to generate shake maps of intensity measures IM, the multivariate lognormal distribution can be adopted. The median values of IM at sites of interest are calculated from the GMPE, whereas their variances are based on the intra-event components. The prediction errors ε in the GMPE are spatially correlated; the correlation coefficient matrix has diagonal elements equal to 1 and off-diagonal elements equal to the correlation coefficient ρ. The correlation coefficient can be calculated using the following equation (Goda and Atkinson, 2010): (18)ρi,jΔ=maxγ⋅exp−α⋅Δβ−γ+1,0 where Δ is the distance between the points i and j, while α, β, and γ are period-dependent model parameters that can be found in Goda and Atkinson (2010).

For each stochastic event, a tsunami simulation is carried out in order to compute the maximum inundation intensity measure (Figure 2C). To optimize the computational time, the subduction plane is discretized into sub-faults of 10 km × 10 km (Figure 3A), and for each sub-fault, the seafloor displacement corresponding to 1 m of slip is calculated using analytical equations by Okada (1985) and Tanioka and Satake (1996). Subsequently, for each simulated earthquake slip (i.e., event), the overall seafloor displacement field is estimated by scaling and summing the seafloor deformation fields of all individual sub-faults that make up the event.

Tsunami modeling is then carried out using a well-tested numerical code of Goto et al. (1997) that is capable of generating offshore 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. The numerical tsunami calculation is performed for duration sufficient to model the most critical phases of tsunami waves (i.e., 2 h). The integration time step is determined by satisfying the CFL condition; it depends on the bathymetry/elevation data, and their grid sizes and is typically between 0.1 and 0.5 s. For the simulation, it is possible to obtain the maximum tsunami intensity measures of interest (i.e., tsunami height, tsunami velocity, etc.) for one or more specific locations along the coast. The results can also be used to evaluate aggregate tsunami hazard parameters, such as inundation areas above a certain depth.

A complete dataset of bathymetry/elevation, coastal/riverside structures (e.g., breakwater and levees), and surface roughness is obtained from the Miyagi prefectural government. The data are provided in the form of nested grids (1350 m–450 m–150 m– 50 m), covering the geographical regions of Tohoku. The ocean-floor topography data are based on the 1:50,000 bathymetric charts and JTOPO30 database developed by Japan Hydrographic Association and based on the nautical charts developed by Japan Coastal Guard. The tidal fluctuation is not taken into account in this study. The elevation data of the coastal/riverside structures are primarily provided by municipalities. 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, Homma’s overflowing formulae are employed. In the tsunami simulation, the bottom friction is evaluated using the Manning’s formula. The Manning’s coefficients are assigned to computational cells based on national land use data in Japan: 0.02 m^−1/3s for agricultural land, 0.025 m^−1/3s for ocean/water, 0.03 m^−1/3s for forest vegetation, 0.04 m^−1/3s for low-density residential areas, 0.06 m^−1/3s for moderate-density residential areas, and 0.08 m^−1/3s for high-density residential areas.

For each value of magnitude, the simulations are used to evaluate the term P(IM ≥ im|M) for the location of interest. Such probability is represented by the CCDF of the IM (Figure 2D). Specifically, the CCDF of the IM (i.e., spectral acceleration or tsunami inundation) is obtained as the Kaplan–Meier estimator (Kaplan and Meier, 1958), for which the variance can be calculated through the Greenwood’s formula (Greenwood, 1926), and therefore, a confidence interval around the central estimate can be obtained. In this study, the 95% confidence interval is considered.

The curves obtained in the previous step for each magnitude are then multiplied by the probabilities corresponding to the related magnitude and eventually are summed up (Figure 2E). Also in this case, three curves are obtained, one corresponding to the central value and two for the confidence interval. The final hazard curves, representing the mean annual rate of occurrence of specific values of earthquake–tsunami intensity measures, are obtained by multiplying the previous three conditional curves (for each hazard) by the occurrence rate of events with magnitudes greater than the minimum magnitude considered in the magnitude–frequency distribution.

Short or incomplete records lead to biased estimation of the hazard parameters, especially when conventional statistical methods are used (Lamarre et al., 1992). To investigate the effect of the number of simulations on the final hazard estimation, a bootstrap procedure is carried out by randomly sampling m values from the original sample containing n elements (with m ≤ n). This provides a pool of different samples of independent and identically distributed random variables, whose distribution function is the same as that of the original sample. For each generated sample, statistics of the parameter of interest (e.g., mean, median, and different percentiles) are then computed. The ensemble of such estimates can be used to quantify the uncertainty in the parameter value.

Figures 8 and 9 show five percentiles (i.e., 5th, 25th, 50th, 75th, and 95th) of the spectral acceleration and wave height, respectively; such intensity measures are calculated for the site in Sendai by considering different magnitude values (i.e., 7.5, 8.0, 8.5, and 9.0) as a function of the number of simulations. Moreover, Figure 8 shows bootstrap results for the seismic case considering two different GMPEs. The analysis is carried out for the maximum sample size of n = 500 simulations. The bootstrap procedure is then applied considering the number of simulations m varying between 1 and 500. For each trial number of simulations m, 1000 Monte Carlo samples are generated, and then the percentile curves are obtained as the mean value of such simulations. Based on Figure 8, it can be concluded that for the considered seismic case, 200 simulations is sufficient to observe stable percentiles for all magnitude values and for both GMPEs considered. Also in this case, for increasing values of magnitude, there is a decreasing trend of variability of the simulated values due to the reduction in variability of the closest distance (Figure 7C).

Tsunami simulation results show that the 50th percentile curves are stable after 100 simulations for all the considered magnitude values. To obtain stable estimates of the high percentiles, a larger number of simulations are needed (the red dotted line in Figure 9). 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 what was observed for the rupture distance (Figure 7C), i.e., when the magnitude is relatively small, the variability of the inundation intensity measures are increased because the rupture area can move more freely within the fault plane. In turn, when the magnitude is large, the fluctuation of the rupture area is more constrained.

Considering the results shown in Figures 8 and 9, 300 stochastic simulations are carried out for the final coupled multi-hazard simulation process in the following. The calculation can be completed in less than 1 week using a conventional workstation with parallel processing.

For each value of 7 magnitudes (i.e., 7.5, 7.75, 8.0, 8.25, 8.5, 8.75, and 9.0), 300 sets of the source parameters θ are generated using the scaling relationships by Goda et al. (2016). Figure 4 shows the simulated source parameters (the green dots). Simulated data are in agreement with the source parameter distributions (i.e., green dots are well clustered within the confidence interval of the scaling relationships). Then, 300 simulations are carried out for the earthquake hazard and tsunami hazard analysis, starting from the same stochastic source models.

The CCDFs (Figure 2D) in terms of PGA, S_a(T = 0.3 s), and S_a(T = 3 s) are shown in Figures 10A,C,E, respectively, for all the magnitude values analyzed. The analogous CCDFs for the tsunami wave height are presented in Figure 11A. Figures 10B,D,F and 11B show the CCDFs, weighted by the probability values obtained from the discretized Gutenberg–Richter relationship (Figure 3D).

As shown in Figure 2E, for each IM [i.e., PGA, S_a(T), h, etc.], the summation of the curves presented in Figures 10B,D,F and 11B, multiplied by λ(M ≥ 7.375) = 0.183, leads to the final hazard curves. Figures 12A–C shows the final hazard curves, and the 95% confidence interval, for PGA, S_a(T = 0.3 s), and S_a(T = 3 s) obtained using the two GMPEs, i.e., Abrahamson et al. (2016), represented with blue lines, and Morikawa and Fujiwara (2013), represented with red lines. In the same figures, the mean seismic hazard curves are represented with black lines. Similarly, Figure 12D shows the final tsunami hazard curve and its 95% confidence interval that is very tight around the central estimate curve. It is noteworthy that the steep slope of the final tsunami hazard curve for wave heights greater than 10 m is because 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 proposed procedure also facilitates the construction of uniform hazard spectra (UHS) for seismic hazard. By repeating the seismic simulations for several spectral accelerations (Figure 13A) and considering specific values of the mean annual rate (e.g., 2, 5, and 10% in 50 years), it is possible to obtain the UHS, as shown in Figure 13B. The hazard curves and spectra are jagged because of the limited number of simulations (i.e., 300). Since the seismic simulations are less time consuming with respect to the tsunami simulations, it is possible to increase the number of simulations; in particular, for each scenario, several thousands of simulations can be conducted with a low additional computational effort. Figure 13C shows the seismic hazard curves obtained by performing 300 × 5000 simulations (i.e., 5000 simulations for each stochastic simulation). By increasing the number of simulations, the UHS become smooth (Figure 13D), and the confidence intervals around the central estimate hazard curves become narrow.

As a byproduct of the procedure, the earthquake–tsunami hazard deaggregation is obtained. Deaggregation shows the relative contributions of dominant seismic scenarios to the specified hazard levels and can be represented in terms of distance and magnitude. The deaggregation of the two hazards for the same mean annual rate of occurrence is demonstrated. Figure 14 shows the deaggregation results for Sendai by considering PGA (Figures 14A,C) and tsunami inundation height (Figures 14B,D) corresponding to two values of mean annual rate of occurrence (i.e., 63 and 10% in 50 years). The 10% in 50-year hazard level (corresponding to an event with 475-year return period) is commonly used to describe the life safety limit state, whereas the 63% in 50-year hazard level (corresponding to an event with 50-year return period) corresponds to the damage control limit state (CEN, 2004). It is worth noting that only large magnitude events contribute to higher values of IMs. Moreover, as observed before, the larger the magnitude is, the less the distance is influential. Finally, it is interesting to observe that the combinations of magnitude and distance that affect the seismic hazard and tsunami hazard differ significantly.