Quantification of source uncertainties in Seismic Probabilistic Tsunami Hazard Analysis (SPTHA)

Abstract

We propose a procedure for uncertainty quantification in Probabilistic Tsunami Hazard Analysis (PTHA), with a special emphasis on the uncertainty related to statistical modelling of the earthquake source in Seismic PTHA (SPTHA), and on the separate treatment of subduction and crustal earthquakes (treated as background seismicity). An event tree approach and ensemble modelling are used in spite of more classical approaches, such as the hazard integral and the logic tree. This procedure consists of four steps: (1) exploration of aleatory uncertainty through an event tree, with alternative implementations for exploring epistemic uncertainty; (2) numerical computation of tsunami generation and propagation up to a given offshore isobath; (3) (optional) site-specific quantification of inundation; (4) simultaneous quantification of aleatory and epistemic uncertainty through ensemble modelling. The proposed procedure is general and independent of the kind of tsunami source considered; however, we implement step 1, the event tree, specifically for SPTHA, focusing on seismic source uncertainty. To exemplify the procedure, we develop a case study considering seismic sources in the Ionian Sea (central-eastern Mediterranean Sea), using the coasts of Southern Italy as a target zone. The results show that an efficient and complete quantification of all the uncertainties is feasible even when treating a large number of potential sources and a large set of alternative model formulations. We also find that (i) treating separately subduction and background (crustal) earthquakes allows for optimal use of available information and for avoiding significant biases; (ii) both subduction interface and crustal faults contribute to the SPTHA, with different proportions that depend on source-target position and tsunami intensity; (iii) the proposed framework allows sensitivity and deaggregation analyses, demonstrating the applicability of the method for operational assessments.

Probabilistic forecasting, Tsunamis, Earthquake interaction, forecasting, and prediction, Europe

1 INTRODUCTION

Computationally based Probabilistic Tsunami Hazard Analysis (PTHA) is performed by explicit numerical modelling of tsunamis generated by potential sources. This approach is often preferred to the empirical PTHA (Geist & Parsons 2006), built upon observations of past tsunamis, due to the scarcity of tsunami observations (Geist & Lynett 2014), which however need to be considered at least for ex post statistical testing, or as a complementary tool, for example, in a Bayesian framework (e.g. Parsons & Geist 2009).

In this work, we propose a methodology that is applicable for computationally based PTHA for any type of tsunami source. However, we emphasize only the implementation for tsunamis of tectonic origin. We refer to this component of PTHA as Seismic PTHA (SPTHA; Lorito et al.2015).

In the computationally based SPTHA methods, there are several open issues. We here try and approach two of the most relevant. First, all the potential seismic sources should be modelled, including the full aleatory variability of source parameters, which has an intrinsically high computational cost that needs to be limited with ad hoc simplifications (see discussion in Lorito et al.2015). Second, detailed information on all potential faults, including their geometry and earthquake occurrence rates, are highly inhomogeneous in space (Basili et al.2013), calling for an adaptive approach to fault modelling.

Considering only large earthquake sources in subduction zones is assumed in several cases a viable option (e.g. Annaka et al.2007; González et al.2009). However, discarding crustal seismicity could result in underestimating the hazard, at least in complicated source regions such as the Caribbean (Parsons & Geist 2009; Omira et al.2015), eastern Indonesia, the Philippines and New Guinea (Løvholt et al.2012b; Horspool et al.2014), the North East Atlantic, the Mediterranean and its connected seas (Baptista & Miranda 2009; Ozel et al.2011; Mitsoudis et al.2012). Conversely, constraining recurrence and geometry of all the possible crustal faults represents a challenge, in particular in terms of spatial completeness (Basili et al.2013), or associated computational burden (Geist & Lynett 2014; Lorito et al.2015). In some applications (e.g. Sørensen et al.2012), the completeness of seismic sources is dealt with the definition of seismic zones that identify seismogenic volumes in which seismicity is possible everywhere, at least to a possible given maximum, similarly to some Probabilistic Seismic Hazard Analysis (PSHA) approaches (e.g. Stucchi et al.2011; Stirling et al.2012; Giardini et al.2013; Hiemer et al.2014). However, this approach neglects most of the available information on the geometrical and rheological properties of the subduction zones and of the major faults, which influence at the first order the features of the tsunami, at least in the near-field (Geist 2009).

The quantification of epistemic uncertainty is typically addressed by means of logic trees (Geist & Lynett 2014), where a set of Mutually Exclusive and Collectively Exhaustive (MECE) alternative models are enumerated, and a subjective probability of being the best model is assigned to each of them (Scherbaum & Kuehn 2011). The development of such alternatives is very challenging from a theoretical and practical point of view (Bommer & Scherbaum 2008). Also, even if many practitioners interpret the outcome of logic trees through distributions and percentiles (Abrahamson & Bommer 2005; Field et al.2014), other practitioners claim that only the mean hazard is the true hazard, and distributions of outcomes do not have any probabilistic meaning (Musson 2005, 2012). Moreover, exploring all the alternatives in a logic tree may become very challenging from a computational point of view (e.g. Field et al.2005). This is particularly true for computationally based SPTHA. As a consequence, effective quantifications of epistemic uncertainty in SPTHA through logic trees are quite rare. For example, they have been totally neglected in the multi-hazard GAR15 assessment (UNISDR 2015), which nevertheless goes all the way to the risk assessment.

An alternative approach to the quantification of epistemic uncertainty is based on Bayesian methods (e.g. Parsons & Geist 2009; Grezio et al.2010, 2012; Knighton & Bastidas 2015), which jointly assess aleatory and epistemic uncertainties. In particular, prior distributions take into account epistemic uncertainty, either by imposing a subjective level of confidence on a given assessment (Parsons & Geist 2009; Grezio et al.2010, 2012), or through the analysis of epistemic uncertainty on model parameters (Knighton & Bastidas 2015). Then, prior distributions are combined with available and independent past data, to evaluate the final posterior distributions. Noteworthy, in the context of PSHA, Selva & Sandri (2013) reconciles logic tree and Bayesian approaches, by adopting logic tree outcomes to estimate prior distributions, to be then combined to past data.

Here, we propose a method for the quantification of both aleatory and epistemic uncertainties in SPTHA that is based on two pillars. The first pillar is the Event Tree (ET). In Lorito et al. (2015), we already proposed an ET for the treatment of aleatory uncertainty on a given subduction zone, to which we will refer here as ‘Interface Seismicity’ (IS). Here, we extend this ET in order to include the treatment of the crustal seismicity or, more generally, of the seismicity which does not occur on a subduction interface, here termed ‘Background Seismicity’ (BS). In a sense, we combine the two above-mentioned approaches to aleatory uncertainty in a more general scheme providing a specific treatment for each of the two kinds of seismicity. The main differences between IS and BS are that the expected maximum magnitudes for IS are typically much larger than that for BS; they usually follow different scaling laws; the amount and the quality of the geometrical information about IS are substantially different from the information on BS, which is expected to have a more sparse behaviour in terms of source geometry, and the knowledge about complex fault networks may be (perhaps intrinsically) incomplete. Such differences call for a different treatment of IS and BS. Note that the proposed (S)PTHA computational framework significantly deviates from the one proposed by Geist & Lynett (2014). The parameterization and the probability quantification of sources are here considered jointly, since these processes are strongly correlated, and a discrete ET is used in place of the hazard integral. Additionally, the tsunami propagation up to some depth offshore is explicitly separated by the inundation phase. This enables the computational reduction that allows scaling down regional to site-specific SPTHA (Lorito et al.2015).

The second pillar of the proposed methodology is an ensemble modelling approach to the treatment of the epistemic uncertainty (Marzocchi et al.2015). Ensemble modelling treats the alternative models as an unbiased sample, though not necessarily MECE, from a parent distribution describing the epistemic uncertainty. It is based on a clear taxonomy of uncertainty in hazard assessments (Marzocchi & Jordan 2014), and it solves some drawbacks of logic trees, like the interpretation of the statistics of outcomes and most of computational issues. It also provides a coherent theoretical framework compatible with Bayesian approaches, in which past data can be either used to update prior distributions (Selva & Sandri 2013) or to objectively estimate weights of alternative models (Marzocchi et al.2012). Here, we demonstrate that it enables us a full exploration of epistemic uncertainty in SPTHA.

In the next sections, we first introduce a general framework for SPTHA (Section 2), presenting the new method for treating aleatory variability (Section 2.1) and the treatment of epistemic uncertainty through ensemble modelling (Section 2.2). We then apply the proposed SPTHA framework to a case study for the Ionian Sea, Central Mediterranean (Section 3). In this application, we perform several sensitivity tests, also explicitly comparing the proposed method with the two end members of SPTHA described above (either considering subduction zones only, or using a seismicity homogenously distributed in volumes). Finally, we analyse the resulting Hazard Curves (HCs), simultaneously quantifying both aleatory and epistemic uncertainties through ensemble modelling. Our application has a purely explanatory goal and several over-simplifications have been made. Although the presented results cannot be regarded to as an actual hazard assessment and cannot yet be considered as an effective input for decision-making, our case study provides a guidance towards reconciling the throughout treatment of uncertainties with the practical feasibility of computationally based SPTHA.

2 THEORETICAL FRAMEWORK FOR SPTHA

The proposed SPTHA method is schematically summarized in Fig. 1. The assessment is based on four STEPS. We point out that these four steps are in principle applicable as well to PTHA for tsunamis due to any type of source, not only seismic.

Figure 1.

Flow diagram illustrating the computational scheme adopted for SPTHA. More details for each STEP have been described in the main text.

STEP 1 regards the development of the source model. The goal is defining a hierarchical discretization of the source parameters space within an ET (Newhall & Hoblitt 2002; Lorito et al.2015), to enable a full exploration of the aleatory variability. Source parameters are hierarchically ordered, so that each parameter realization may depend on the higher levels in the ET, but not on the lower levels. At each level, the expected natural variability is described through a discrete number of possible values and associated to a discrete Probability Density Function (dPDF), conditioned to the realizations at previous levels. The discretization is assumed fine enough and not affecting the final results (Baker & Cornell 2008). Ideally, a finer discretization would be required for portions of the parameter space to which the hazard is more sensitive (Lorito et al.2015). This is a potential advantage over more classical approaches based on the solution of the hazards integral (e.g. Geist & Parsons 2006). The result of STEP 1 is the definition of a complete set of source scenarios σ_s, and their annual rates λ(σ_s), through an ET which allows the exploration of the aleatory uncertainty. The epistemic uncertainty in STEP 1 is explored by defining alternative approaches for estimating the probability at each ET level, or even alternative ET formulations, resulting in a number of alternative assessments of λ(σ_s).

STEP 2 regards the propagation of tsunamis generated by each individual source identified in STEP 1. Since the number of such individual sources may be very large (on the order of 6 × 10⁶ for the Ionian sea, see Section 3.1.2), a set of pre-calculated tsunami scenarios, resulting from Gaussian-shaped initial sea level elevation centred at dense and equally spaced points all over the considered source regions, is stored in a database of tsunami waveforms. These scenarios are linearly combined for approximating the offshore tsunami impact, at a given isobath, due to any source considered in STEP1. The result of STEP 2 consists of the collection of the maximum wave height h_max at discrete points, in our application along the 50 m isobath, for each individual scenario defined at STEP 1. As for STEP 1, also in STEP 2 alternative implementations are possible, resulting in a number of alternative assessments of the h_max due to each individual source.

STEP 3 deals with tsunami amplification during the shoaling phase and inland impact assessment. Approaches in which high-resolution metric-scale inundation is simulated for each and every source scenario defined at STEP 1 are clearly unaffordable. Since two options are possible, we split this part in STEPS 3a and 3b. In either case, as for STEPs 1 and 2, alternative modelling approaches are possible, leading to a number of alternative assessments of the tsunami hazard intensity inland.

STEP 3a. For regional-scale hazard analysis, a crude approximation of maximum run-up can be obtained via amplification factors, depending on the features of the incident waves, of the coastal slope and on bottom friction. Examples are the Green's law (Synolakis 1987, 1991) and its modifications (see e.g. discussion in Geist 1999), and the technique used for the tsunami component of UNISDR (2015), based on Løvholt et al. (2012a). For long non-breaking waves, these amplification factors may be used to project the offshore h_max up to the shoreline (Carrier & Greenspan 1958; Synolakis 1987; 1991; Løvholt et al.2013), and this values can be considered either just an assessment at the shoreline or a proxy for the maximum run-up on the coast behind. Where relevant, for example, for coastal plains, shoreline h_max values can be input to empirical relationships considering also bottom frictional damping for obtaining inundation distance and consequently run-up.

STEP 3b. For site specific hazard, explicit inundation modelling is mandatory, in order to get the spatial distribution of the desired hazard intensity, such as maximum flow depth, momentum flux, or composite impact metrics obtained by combination of the fundamental variables (e.g. chap. 7 in TPSWG 2006). STEP 2 offshore results may provide the input for an ex post filtering procedure, which allows to subsequently perform high-resolution inundation modelling for a limited subset of scenarios (Lorito et al.2015). This technique makes the computational problem more affordable, while introducing a manageable additional source of uncertainty in the probabilistic maps. Other techniques aiming at reducing the computational burden associated to probabilistic inundation maps have been proposed as well (e.g. Thio et al.2010; Thio & Li 2015; Wei et al.2015). The result of STEP 3 consists of the collection of the inland tsunami hazard intensities for each individual scenario defined at STEP 1.

STEP 4. The combination of STEPS 1 and 2 provides the assessment of offshore HCs. STEP 3a (for maximum run-up) or 3b (for other metrics) assess inland HCs. In either case, the annual rate of exceedance of a threshold ψ for the hazard intensity Ψ at a location x, can be defined in terms of the annual rates λ(σ_s) of the seismic sources σ_s and the modelling of the tsunami propagation from the seismic source to x, that is in a general formulation:

\begin{equation}{\lambda ^{{\rm{Tot}}}}\left( {{\rm{\Psi }} \ge \psi ,{\rm{\ }}\boldsymbol{x}} \right) = \mathop \sum \limits_s \lambda \left( {{\sigma _s}} \right) \cdot {\rm{Pr}}({\rm{\Psi }} \ge \psi |{\sigma _s},\boldsymbol{x}). \end{equation}

(1)

Assuming all terms in eq. (1) are stationary over the exposure time ΔT and assuming a Poisson process, the HCs for the SPTHA at the target location x for the selected ΔT can be assessed as

\begin{equation}\Pr \left( {{\rm{\Psi }} \ge \psi ,\boldsymbol{x},{\rm{\Delta }}T} \right) = 1 - {\rm{exp}}\left( { - {{\rm{\lambda }}^{{\rm{Tot}}}}\left( {{\rm{\Psi }} \ge \psi ,\boldsymbol{x}} \right) \cdot {\rm{\Delta }}T} \right).\end{equation}

(2)

The two factors in eq. (1) are related, respectively, to the source (λ(σ_s)), and to the propagation from the source to the target location x (Pr(Ψ ≥ ψ|σ_s, x)). Evaluating eq. (2) for several values of the threshold ψ allows constructing the HCs at any location x. As discussed above, STEPS 1 to 3 can be repeated for each possible alternative implementation, regarding both source (e.g. different ways of assessing seismicity rates in STEP 1) and tsunami modelling (e.g. using different numerical codes in STEPS 2 and 3). Since the total set of scenarios σ_s to be modelled is set in STEP 1 within the ET(s), all the alternatives can be implemented by re-weighting the propagation results of STEPS 2 and 3, largely reducing the overall computational effort. STEP 4 is then devoted to rank all these alternatives and integrate them into a single model in the framework of ensemble modelling (Marzocchi et al.2015), where aleatory and epistemic uncertainties are consistently and simultaneously quantified. At this STEP, potential correlations among alternative models are treated, and the final results are prepared in terms of distributions that comprehensively represent all the uncertainties to be communicated to end-users. The ranks of models are expressed through weights assigned to each alternative formulation, ideally quantified for long-term analyses through experts’ elicitation experiments (e.g. SSHAC 1997; Saaty 2008).

The development of the ET in STEP 1 is discussed in details in Section 2.1. The development of the ensemble model for uncertainty assessment of STEP 4 is discussed in Section 2.2. The framework for the linear offshore tsunami propagation will be discussed in details in Molinari et al. (in preparation), particularly as it regards the uncertainty introduced by such an approach on the hazard intensity estimates. The basics of this framework are however briefly summarized in Section 3.1.2 for the application to the Ionian Sea. Uncertainty in numerical modelling of tsunami propagation will be instead here neglected. STEPS 3a and 3b have been discussed in details in Lorito et al. (2015). Here, in the application, we implement STEP 3a only adopting the Green's law for estimation of HCs at the shoreline (1 m isobath), which can be considered a very rough proxy for run-up, whereas inundation calculations are not performed.

2.1 STEP 1: development of the ET

The ET is focused on the source term λ(σ_s) in eq. (1), allowing a robust description of source variability (definition of the set of scenarios) and statistics (assessment of the mean annual frequency of each scenario). The ET for subduction zones proposed by Lorito et al. (2015) is here extended by considering both subduction Interface Seismicity (hereinafter, IS) and crustal Background Seismicity (hereinafter, BS).

The whole source area is divided into a number of statistically independent seismic regions R_i (regionalization), which may or may not contain a source of IS. In each region, an independent ET (Fig. 2) is developed. The first two LEVELS of the ET (Fig. 2a) are:

1. Magnitude M or corresponding seismic moment m₀;

2. Seismicity class C (here, indicating either IS or BS);

Figure 2.

Event Tree (ET) adopted for SPTHA. The part of the ET in common to all seismicity classes C_k is reported in panel (a). In this paper, C_k considers two classes: BS and IS. In panels (b) and (c), the BS-branch and the IS-branch are reported separately. To increase the readability, a generic number n of discrete values for the dPDF at each LEVEL is reported, even if it is not necessarily equal at all LEVELS. Also, at many branches the label ‘clone’ highlights that a subtree equal to the ones at the same LEVEL should be reported.

At LEVEL 1, the annual rates λ_i(M_j) are computed in each region R_i for a finite set of different magnitude levels M_j, each of them representing an actual interval of M. Different alternative procedures can be adopted for the assessment of λ_i(M_j). They are mainly based on fitting the total productivity and the frequency-size relationship (the a- and b-values of a GR distribution). Many different assumptions can be made, ranging from the choice of a distribution (e.g. truncated or tapered Pareto) to the setting of different procedures to quantify its parameters (e.g. maximum magnitude M_max), which may come from seismic catalogues and/or geometrical constraints. Since λ_i(M_j) is assessed jointly for all magnitude levels, it is theoretically possible to consider the statistical dependency of the distribution parameters (Bommer & Scherbaum 2008; Keller et al.2014). In addition, given the potential impact on SPTHA of the tails of the magnitude-frequency distributions, an important effort should be put for an operational assessment to constrain the rates of large events (e.g. Geist & Parsons 2014; Rong et al.2014); dealing with this topic is however beyond the scope of this study. The choice of starting the ET in each region from setting the frequency-size relationship is of paramount importance. Indeed, this enables an overall control of the energetic balance of each region, independently by all the other LEVELS and, in particular, by the separation in different seismicity classes.

At LEVEL 2, we separate the different classes of seismicityC_k. In our case, the two seismicity classes C_k are BS and IS. More classes could be added, if data support further distinctions and significant improvements are foreseen for tsunami modelling. The set of the seismicity classes represents a partition of the whole seismicity, so that |$\sum \nolimits_k \Pr ( {{C_k}{\rm{|}}{R_i},{M_j}} ) = 1$|⁠. This means that we define the probability that a given random earthquake in the region R_i of a given magnitude M_j is either in IS or BS, with probability |$\Pr ( {{C_k} = {\rm{IS|}}{R_i},{M_j}} )$| and |$\Pr ( {{C_k} = {\rm{BS|}}{R_i},{M_j}} )$| respectively. This LEVEL is not trivial only for the regions including IS. If the region does not include any subduction interface, the dPDF is simply |$\Pr ( {{C_k} = {\rm{IS|}}{R_i},{M_j}} ) = 0$| and |$\Pr ( {{C_k} = {\rm{BS|}}{R_i},\ {M_j}} ) = 1$|⁠, for all M_j. Otherwise, |$\Pr ( {{C_k} = {\rm{IS|}}{R_i},{M_j}} )$| is expected to be a function of M_j, with values likely closer to 1 for high magnitude events (e.g. M_j> 8), and close to a background level for relatively small magnitudes (e.g. M_j< 5). Even though this probability may be difficult to constrain, LEVEL 2 is fundamental for a separate treatment of IS and BS at the following LEVELS and thus for a more focused use of the available information. On one hand, its inclusion is itself a source of epistemic uncertainty. On the other hand, we argue that we may introduce a much larger source of epistemic uncertainty by neglecting it. These speculations are quantitatively tested in the application, in Section 3.

The overall goal of this first part of the ET is assessing the annual rates for all magnitudes M_j and all classes of seismicity C_k in each independent region. Indeed, through LEVELS 1 and 2, the overall rate of earthquakes can be written through the two separate contributions of IS and BS:

\begin{eqnarray}{\lambda _i}\left( {{M_j}} \right) &=& {\lambda _i}\left( {{M_j}} \right) \cdot [\Pr \left( {{C_k} = {\rm{IS|}}{R_i},{M_j}} \right)\nonumber\\ && +\, \Pr \left( {{C_k} = {\rm{BS|}}{R_i},{M_j}} \right)]\nonumber\\ &=& \lambda _i^{{\rm{IS}}}\left( {{M_j}} \right) + \lambda _i^{{\rm{BS}}}\left( {{M_j}} \right). \end{eqnarray}

(3)

After LEVEL 2, specific ET branches can be developed separately for each seismicity class C_k, aiming at describing all the aleatory variability of the sources in each class. Noteworthy, this formulation allows separating the contribution of BS and IS also in the overall hazard, as formulated in eqs (1) and (2).

The goal of each C_k-branch (subsequent to LEVELS 1 and 2, see Fig. 2) is to define a MECE set of possible scenarios representing all potential seismic sources in each region R_i, for each seismicity class C_k and each magnitude M_j. The MECE assumption implies that |$\sum \nolimits_l \Pr ( {\sigma _l^{( {{R_i},{M_j},{C_k}} )}{\rm{|}}{R_i},{M_j},{C_k}} ) = \ 1$| for all i, j and k, where with |$\sigma _l^{( {{R_i},{M_j},{C_k}} )}$| we indicate the l-th seismic source among all the possible ones in region R_i, of class C_k and with magnitude M_j. The overall mean annual rate of exceedance of a threshold ψ for the hazard intensity at a location x, can then be defined in terms of the separate contribution of all the possible seismic sources, so that

\begin{eqnarray} &&{\lambda ^{\rm{Tot}}\left( \rm{\Psi } \ge \psi ,\boldsymbol{x} \right) = \sum_k {\lambda ^{C_k}}(\Psi \ge \psi , \boldsymbol{x})}\nonumber\\ &&{\quad = \sum_{\boldsymbol{k}} \Bigg\{ \sum_{i,j,l} \Bigg[ \lambda _i^{{C_k}} \left( {M_j} \right)\Pr \left( {\sigma _l^{\left( {R_i},{M_j},{C_k} \right)}{\rm{|}}{R_i},{M_j},{C_k}} \right)}\nonumber\\ &&{\qquad{\rm{Pr}}({\rm{\Psi }} \ge \psi |\sigma _l^{\left( {R_i},{M_i},{C_k} \right)},\boldsymbol{x}) \Bigg] \Bigg\}} \end{eqnarray}

(4)

with C_k = BS, IS, i and j running over all regions and magnitudes respectively and l ranging over all MECE scenarios in each C_k-branch. Eq. (4) identifies the three main terms, corresponding to energy balance expressed through the earthquake rate |$\lambda _i^{{C_k}}( {{M_j}} )$|⁠, the probability of a specific earthquake scenario |$\Pr ( {\sigma _l^{( {{R_i},{M_j},{C_k}} )}{\rm{|}}{R_i},{M_j},{C_k}} )$| and that of hazard metric exceedance estimated with tsunami numerical propagation |${\rm{Pr}}({\rm{\Psi }} \ge \psi |\sigma _l^{( {{R_i},{M_j},{C_k}} )},\boldsymbol{x})$|⁠, respectively. Also, eq. (4) represents the development of eq. (1) through the presented ET, with the index s corresponding to all combinations of {i, j, k, l} (region R_i, magnitude M_j, seismicity class C_k and earthquake mechanism σ_l). The details on the formulations of BS- and IS-branches are reported in the following sections.

2.1.1 BS-branch

BS includes all the seismic sources that are outside the subduction interfaces of the active subduction zones forming IS (or, more in general, the interface of major fault systems). Therefore, it represents a broad class of sources occurring within the seismogenic layers. Within the BS-branch, we define the following five LEVELS (Fig. 2b):

1. Centre of the fault (x, y);

2. Depth of the centre of the fault z;

3. Fault geometry/earthquake mechanism, that is, strike, dip and an average rake (s, d, r);

4. Maximum rupture area A;

5. Seismic moment distribution on the finite fault |${}^{^{-\!-\!\rightharpoonup}}\!\!\!\!\!\!\!\!\!\!{\Delta m}$|⁠, defined in terms of local deviation from the average moment per unit area (m₀/A) and from the average rake r.

Note that the geometrical parameters of BS are treated through a purely probabilistic method, that is, no strong assumptions are made regarding for example the fault geometry. For simplicity, hereinafter we avoid introducing further indexes to indicate the discrete intervals for all above parameters.

At LEVELS BS-1 and BS-2, the position of the fault centre is established, both for its horizontal position (x, y) and its depth (z). The analysis of potential (x, y) may be performed on a regular grid identifying a discrete number of cells (e.g. Grezio et al.2012). The goal is the assessment of the spatial probability |$\pi _1^{{\rm{BS}}}\ = \Pr ( {x,y{\rm{|}}{R_i},{C_k} = {\rm{BS}}} )$|⁠, representing a dPDF over all cells in each region. Even if the probability value at one LEVEL is theoretically conditioned to the realization at all previous LEVELS, here we omit to report the magnitude M_j, since the spatial probability |$\pi _1^{{\rm{BS}}}$| is typically assessed independently from it. This kind of simplifications will be done at all LEVELS. Different methods have been proposed in literature, mainly based on the analysis of seismic catalogues (uniform, smoothed seismicity, smoothed seismicity plus faults, simple sampling from past occurrences, Grezio et al.2012; Sørensen et al.2012; Hiemer et al.2014). Also the assessment of the depth probability |$\pi _2^{{\rm{BS}}} = {\rm{Pr}}(z|{R_i},{M_j},{C_k} = {\rm{BS}},A,\ x,y)$| can be based on the analysis of the seismic catalogues. A uniform distribution is often adopted for offshore earthquakes because it is difficult to constrain this probability reasonably well and homogeneously enough in the considered domain (Grezio et al.2012; Sørensen et al.2012).

At LEVEL BS-3, a joint dPDF must be established for the focal mechanism, that is, strike, dip and average rake (s, d, r), conditioned to the parameters at the previous LEVELS. For BS, in each specific position any mechanism is theoretically possible. The joint dPDF |$\pi _3^{{\rm{BS}}}\ = \Pr ( {s,d,r{\rm{|}}{R_i},{C_k} = {\rm{BS}},x,y} )$| can be constrained by focal mechanism catalogues. In addition, information on sufficiently well-known local faults can be used to probabilistically constrain the strike, dip and rake of the potential earthquakes in that location.

At LEVEL BS-4, the dPDF for the rupture areaA is established. This area may be assumed rectangular for BS, characterized by a length L and a width W. Given that BS includes mainly crustal faults, the joint dPDF |$\pi _4^{{\rm{BS}}} = {\rm{Pr}}(A|{C_k} = {\rm{BS}},r)$| can be constrained by specific scaling laws for the different tectonic regimes (rake r, e.g. Wells & Coppersmith 1994).

At LEVEL BS-5, the heterogeneities on seismic moment distribution are modelled, assessing the probability |$\pi _5^{{\rm{BS}}} = \Pr (\,{{}^{^{-\!-\!\rightharpoonup}}\!\!\!\!\!\!\!\!\!\!{\Delta m}{\rm{|}}{M_j},{C_k} = {\rm{BS}},x,y,z,s,d,r,A} )$|⁠. At previous LEVELS, only the total seismic moment (LEVEL 1), the average focal mechanism (LEVEL BS-3) and the rupture area (LEVEL BS-4) are treated. From these values we can define an average moment per area unit m_0i/A and an average slip modulus u = m_0i/(A · μ), where μ is the rigidity in each unit area. Large variations of μ are not expected within the rupture area of the crustal faults, so it is reasonable to consider just one single value for the entire fault area. In any case, at LEVEL BS-5, it is possible to model the local deviations for such mean parameters on the finite fault. The heterogeneities in the moment and the mechanism distribution may be particularly important if the target site is in the near field, that is, within a distance of one or two fault lengths. This LEVEL can be simulated with different stochastic techniques (see e.g. Geist & Oglesby 2014). However, this aspect is expected to be more important for IS in subduction zones, characterized by stronger heterogeneities and larger earthquakes.

Once the discretization and the dPDFs are established at all LEVELS, the probability of each possible scenario can be computed as:

\begin{equation}\Pr \left( {\sigma _l^{\left( {{R_i},{M_j},{C_k} = {\rm{BS}}} \right)}{\rm{|}}{R_i},\ {M_j},{C_k} = {\rm{BS}}} \right) = \ \mathop \prod \limits_{n = 1}^5 \pi _n^{{\rm{BS}}}\ \end{equation}

(5)

This equation represents the earthquake mechanisms term in eq. (4) for C_k = BS.

2.1.2 IS-branch

IS includes all the potential interface earthquakes in active subduction zones (but potentially it can be extended to any well constrained fault). Therefore, the IS-branch is developed only in regions where such interfaces exist. IS represents a rather specific class of sources, and, generally, has more available information on geometrical parameters. This allows for a more specialized procedure to treat aleatory uncertainty. Within the IS-branch, we define only the following LEVELS (Fig. 2c):

1. maximum rupture area A on the interface plane and the potential centres of the fault (x', y');

2. seismic moment distribution on the finite fault, defined in terms of local deviation from the average moment per unit area (m₀/A) and from average rake r.

Note that the number of IS LEVELS is smaller than the BS LEVELS because of the different level of determinism that is possible in the analysis of the subduction interfaces. In particular, (i) the coordinates of the potential centres (x, y, z) are directly mapped on the known geometry, so a common reference system (x, y, z) = f(x', y') is defined on the interface; (ii) the definition of potential fault centres (x', y') is made jointly with the definition of potential rupture areas A, since they depend on each other; (iii) the ruptures follow the geometry of the interface (strike and dip) and the average convergence direction (rake), so they are set deterministically for IS (the ET deals with aleatory uncertainty only). Eventually, the uncertainty on such parameters may be modelled as epistemic uncertainty through alternative implementations, as discussed in Section 2.2.

At LEVEL IS-1, the assessment of potential rupture areaA for a given magnitude M_j and potential fault centres (x', y') is made jointly. Indeed, the possibility that a given point of the interface (x', y') is the centre of the rupture depends on the total rupture area for that given magnitude M_j (Lorito et al.2015). Conversely, the modelled rupture area A should be adapted to the complexity of the interface geometry, which varies locally and depends on the central point (x', y'). Therefore, a joint procedure is required to define the potential rupture areas A in each potential fault centre (x', y') and to assess their joint probability of occurrence |$\pi _1^{{\rm{IS}}} = {\rm{Pr}}( {x',y',A{\rm{|}}{R_i},{M_j},{C_k} = {\rm{IS}}} )$|⁠. Differently from BS, the IS rupture is not assumed rectangular especially for large magnitudes, when the rupture is supposed to eventually spread over all the available space. Moreover, the scaling laws are often violated by large tsunamigenic earthquakes, as it was for the M9 Tohoku 2011 earthquake (e.g. Romano et al.2014). At the very least, a large variability is expected around the mean values from the scaling laws for great subduction earthquakes.

At LEVEL IS-2, we may consider potential seismic moment distributions to evaluate the probability |$\pi _2^{{\rm{IS}}} = \Pr (\,{{}^{^{-\!-\!\rightharpoonup}}\!\!\!\!\!\!\!\!\!\!{\Delta m}{\rm{|}}{{{R}}_{i}},{{{C}}_{{k}}} = {\rm{IS}},{{\ x'}},{{y'}},{{A}}} )$|⁠. The same considerations that we reported for the LEVEL BS-5 are valid here, and heterogeneous moment distributions are expected to play a significant role for large subduction earthquakes. Additionally, the rigidity in the shallow regions of subduction zones may vary even up to a factor of 5, and then earthquakes with similar seismic moment may result in substantially different slip amount and tsunami excitation (e.g. Bilek & Lay 1999; Geist & Bilek 2001; Polet & Kanamori 2009). These systematic variations can be modelled in this framework, and we plan to address them in future studies.

Once discretization and dPDFs are established at all LEVELS, the probability of each possible scenario is assessed as:

\begin{equation}\Pr \left( {\sigma _l^{\left( {{R_i},{M_j},{C_k} = {\rm{IS}}} \right)}{\rm{|}}{R_i},\ {M_j},{C_k} = {\rm{IS}}} \right) = \ \pi _1^{{\rm{IS}}} \cdot \pi _2^{{\rm{IS}}}.\end{equation}

(6)

This equation represents the earthquake mechanisms term in eq. (4) for C_k = IS.

2.2 STEP 4: ensemble modelling

Epistemic uncertainty arising from the different alternatives briefly discussed in the previous sections for the various ET LEVELS is here treated in the framework of ensemble modelling, as introduced in Marzocchi et al. (2015) for PSHA. Aleatory uncertainty is quantified by the expected long-run frequencies of random events within the model of the system. Such frequencies are objective probabilities θ, and they can be potentially measured through a well-defined experimental concept. The experimental concept defines collections of data - observed and not observed yet - that are judged to be exchangeable. The long-run frequencies are determined by this data-generating process (Marzocchi & Jordan 2014). Hypotheses about aleatory variability (models for such frequencies) can be tested against observations by frequentist (error-statistical) methods. Epistemic uncertainties measure the lack of knowledge in the estimation of such frequencies. Models assessing θ are often based on expert opinions, thus epistemic uncertainty is described by subjective probabilities. Bayesian methods are appropriate for reducing epistemic uncertainties as new knowledge is gained through observation. The epistemic uncertainty arising from this framework describes ‘the center, the body, and the range of technical interpretations that the larger technical community would have if they were to conduct the study’ (SSHAC 1997, 2012).

To quantify such uncertainty, a finite set of different models {θ_i, ω_i} (i = 1, …, N) can be developed, where θ_i and ω_i are the outcome and the weight of the i-th model. The N different models describe one specific variable of interest θ, representing an estimation of a long-run frequency. In ensemble modelling, each model is considered as a sample of an unknown parent distribution f(θ) that describes the variable θ taking into account simultaneously for the aleatory variability and epistemic uncertainty:

\begin{equation}{\theta _i} \sim f\left( \theta \right) \equiv \left[ \theta \right].\end{equation}

(7)

In order to produce a consistent ensemble model, the sample {θ_i, ω_i} should represent a set of scientifically acceptable models (SSHAC 1997). Differently from logic trees, such models are not assumed MECE and can derive from one or more logic trees, or from a collection of models. The only requirement is that {θ_i, ω_i} represents an unbiased sample of the epistemic uncertainty (see discussion in Marzocchi et al.2015). The weight associated to each model should properly take into account the confidence (based on expert opinion and/or on quantitative evaluation of the forecasting performances) and the possible correlation between the models; examples of the assessments of these weights based on expert opinion are the SSHAC (1997) process and the Analytic Hierarchy Process (Saaty 2008).

Different variables of interest θ may be defined. In probabilistic hazard assessments, a typical experimental concept consists of collecting data on the exceedance of a selected hazard intensity level during N equivalent time intervals in one specific site or region. In this context, the exchangeability is referred to time intervals, meaning that they are considered equivalent and their position in time is exchangeable. In this case, specifically for (S)PTHA, it is convenient to set |${\theta ^{( \psi )}} = \Pr ( {\psi \ge \bar \psi ;\boldsymbol{x},{\rm{\Delta }}T} )$|⁠, as in eq. (2); then, the variable of interest is the hazard curve itself. To derive an appropriate set of models {|$\theta _i^{( \psi )}$|⁠, ω_i}, alternative formulations of all the terms in eq. (2) must be developed. Focusing on the epistemic uncertainty on the source term, we can consider for example alternative implementation of the ET discussed in Section 2.1. In practice, we can consider alternative assessments of the dPDF, alternative discretization procedures at the LEVELS, and alternative configurations of the subduction interface and the relative parameters. All these alternatives may be appropriately represented by means of an Alternative Tree (AT). Note that this is different from a Logic Tree (e.g. Geist & Parsons 2006; Bommer & Scherbaum 2008), since no probabilistic interpretation of branches is considered. Of course, incompatible branches cannot be combined, but this choice does not pose practical or theoretical problems, because a model outside the AT is just another sample for the ensemble with its associated weight. Since the weights are not probabilities they do not have to sum to 1.

For a small sample size, the distribution [θ^(ψ)] can be set through a parametric distribution (e.g. a Beta distribution), replacing such few probability outcomes with a continuous distribution that describes the aleatory variability and the epistemic uncertainty. For a large enough sample size, a non-parametric distribution fitting the samples may be also adopted (if judged scientifically sound, see discussion in Marzocchi et al.2015). In any case, if the number of alternative models is too large to be computationally feasible, a subset of alternative models can be sampled from the original set of alternatives considering the relative weights. This subset may represent a sample for the ensemble. In all cases, the aleatory uncertainty may be assessed through the central tendency values of [θ^(ψ)] (e.g. mean or median), while the epistemic uncertainty may be represented through a confidence interval (e.g. the interval 16th–84th percentiles). The distributions [θ^(ψ)] enable a meaningful test of the hazard results against real independent observations (Marzocchi & Jordan 2014; Marzocchi et al.2015) and a worthy communication of the degree of epistemic uncertainty to the decision makers (Paté-Cornell 1996; SSHAC 1997, 2012).

3 APPLICATION: THE CASE STUDY OF THE IONIAN SEA (CENTRAL MEDITERRANEAN SEA)

This case study emphasizes the seismic source variability analysis in SPTHA, through implementation of STEPS 1 and 4 of Fig. 1 for seismic sources. With the aim of demonstrating the proposed methodology, many choices are here made just for illustrative purposes.

In Section 3.1, we introduce the specific SPTHA settings used for the case study, going through STEPS 1 to 4. In Section 3.2, we analyse the sensitivity of the HCs to alternative source models and to the separation between IS and BS. In Section 3.3, we finally present the results of the SPTHA case study in terms of an ensemble model that simultaneously accounts for both aleatory and epistemic uncertainty.

3.1 Setting the SPTHA for the Ionia Sea

The application area is the Ionian Sea, in the central-eastern part of the Mediterranean Sea, a region tectonically dominated by the interaction of the Africa, Eurasia and Aegean plates. The African oceanic crust is the World's oldest (280 My) and is subducted beneath the European Plate in the Calabrian Arc and beneath the Aegean Plate in the Hellenic Arc (see Fig. 3a). Although the two slabs belong to the same plate, the Calabrian and Hellenic subduction zones are different in terms of both geometry and behaviour. The Hellenic slab dips to the northeast in the Ionian Sea and to the north in the stretch from Crete to Rhodes, at a shallow angle (ca. 20°–30°, e.g. Gesret et al.2011), whereas the Calabrian slab dips to the northwest and is much steeper (ca. 70°–80°, e.g. Chiarabba et al.2008). Convergence rates are on the order of 35 mm yr⁻¹ in the Hellenic Arc (e.g. Reilinger et al.2006) and 1–5 mm yr⁻¹ in the Calabrian Arc (e.g. Devoti et al.2008). Both subduction zones feature wide and thick accretionary wedges. Away from the subduction of the African oceanic crust, the collision of continental crust is active along the coast of northern Greece and Albania, in Sicily, and in the Taranto Gulf. A promontory of continental crust of the African Plate, often referred to as the Adria Microplate, indents the European Plate towards the Adriatic Sea. Intracontinental deformation occurs in the Sicily Channel and carries on towards the Gulf of Sirte (e.g. Corti et al.2006).

Figure 3.

(a) Application area for the Ionian Sea case study, highlighting the source area and the target area. (b) Zoom on the selected target area, highlighting the position of the receivers (small yellow dots) and localities for which hazard curves have been plotted and commented in the main text (red dots). (c) Zoom on the source area, with tectonic regionalization of the Ionian area.

Our SPTHA target area is a segment of the Southern Italian coasts, from the south-easternmost point of Sicily to the east side of Apulia region (see Fig. 3b). In particular, 398 points along the coasts are selected, spaced ∼2 km on the average (yellow dots in Fig. 3b). Detailed analyses are reported on 12 selected coastal points for selected locations, namely Siracusa, Catania, Taormina, and Messina in Sicily; Bova Marina, Locri, Soverato, Crotone, and Sibari in Calabria; Taranto, Leuca, and Brindisi in Apulia (red dots in Fig. 3b).

3.1.1 STEP 1: the event tree for exploring uncertainties on sources

For the sake of conciseness, we report in the main text the most important steps of the application, while additional technical details regarding the SPTHA computations are reported in Appendix A in the Supporting Information.

In Figs 3(a) and (c), the tectonic regionalization for the application source area is shown. The regions are defined according to the tectonic setting. In this figure, we highlight the two subduction interfaces in the application area, namely the Calabrian and the Hellenic arcs. To guarantee homogeneity in the assessments at LEVELS 1 and 2, all regions including parts of the same subduction interface are merged to obtain macro-regions. The macro-region corresponding to the Calabrian arc is formed by regions 12 and 26; the one corresponding to the Hellenic arc is formed by regions 7, 29, and 35. More details on the regionalization are reported in Appendix A.2 (Supporting Information). The assessment at LEVELS 1 and 2 is made using the macro-regions, then at the following LEVELS for the single regions.

At LEVEL 1, the goal is the assessment of the mean annual rates of the earthquakes at the different magnitudes levels in each region/macro-region separately. We define a set of n₁ magnitudes M_j (j = 1, 2, …, n₁) covering the complete range of expected magnitudes within all the source area. Each magnitude level M_j represents the most probable magnitude (the mode) in an interval [M_j, M_{j + 1}].

Several possible methods may be adopted to assess the rates λ(M_j): classical or Bayesian joint assessments of the parameters of the frequency-size distribution (e.g. Keller et al.2014); separated assessments of mean annual rates of earthquakes and of the parameters of a Pareto distribution (e.g. Kagan 2003). For simplicity, we adopt the latter, assessing λ(M_j) for each region through

\begin{eqnarray}\lambda \left( {{M_j}} \right) &=& \lambda \left( {M\ > {M_t}} \right) \cdot \Big\{ \Phi \left( {{M_{j + 1}};{M_t},{M_{\max }},\beta } \right)\nonumber\\ && -\, \Phi \left( {{M_j};{M_t},{M_{\max }},\beta } \right) \Big\},\end{eqnarray}

(8)

where M_t represents a lower threshold for the magnitudes (e.g. Kagan 2002a,b), and it is set to 4.5; λ(m  > M_t) is the mean annual rate for the magnitudes greater than M_t; Φ(M; M_t, M_max, β) is a cumulative distribution function for the magnitudes with parameters M_t, M_max, and β, defined for M > M_t.

We adopt three alternative assessments of λ(M > M_t), two functional forms for Φ with eight alternative choices for setting the distribution's parameters. We assign different integer weights to each alternative choice indicating their relative credibility. In total, we consider 48 alternative assessments (see details in Appendix A.3.1 in the Supporting Information).

At LEVEL 2, the goal is the separation of IS and BS. We adopt the following formulation:

\begin{eqnarray}\Pr \left( {{C_k} = {\rm{IS|}}{R_i},{M_j}} \right) = {a_i}\left( {{M_f}} \right) + \left( {1 - {a_i}\left( {{M_f}} \right)} \right) \cdot f\left( {{M_i} - {M_f}} \right),\!\!\!\!\!\nonumber\\ \end{eqnarray}

(9)

where |${a_i}( {{M_f}} ) = \frac{{{n_{{\rm{IS}}}}( {m > {M_f}} )}}{{{n_{{\rm{IS}}}}( {m > {M_f}} ) + {n_{{\rm{BS}}}}( {m > {M_f}} )}}$|⁠, is the fraction of the total number of events being IS at a given magnitude level M_f, and f(M_i − M_f) is a monotonically increasing function with f(0) = 0, and f(x ≫ 1) = 1.

The parameter a_i(M_f) is assessed counting the IS events in the EMEC seismic catalogue (Grünthal et al.2010, Grünthal & Wahlström 2012) with magnitude M > M_f (see Appendix A.1, Supporting Information). In particular, the events in the catalogue are assigned to IS or BS depending on their hypocentre location with respect to the known interface geometry. Events with hypocentre inside a vertical buffer around the subduction interface are assigned to IS, and all the others are assigned to BS. We consider 5 alternative choices for this buffer: 0, 5 km, 10 km, 15 km and ∞. The first and the last choices correspond to assume that the whole seismicity is BS (a_i = 0 and f(x) = 1, for all x) or IS (a_i = 1), respectively. Also, we adopt two alternative choices for both M_f (5.0 and 6.0) and f(x) (linear or exponential between M_f and 8) providing 12 alternatives (see Appendix A.3.2, Supporting Information). Each alternative allows defining eq. (9) for all macro-regions containing a subduction interface. The earthquake rates of BS are finally divided into each region, according to the fraction of BS events with M > 5 in the catalogue. For all the regions not containing any subduction interface, LEVEL 2 is trivial, and |$\Pr ( {{C_k} = {\rm{IS|}}{R_i},{M_j}} )$| is set to 0.

At LEVEL BS-1, the probability |$\pi _1^{{\rm{BS}}} = {\rm{Pr}}(x,y|{R_i},{C_k} = {\rm{BS}})$| is studied over a regular grid with cells of equal area, approximately 25 km × 25 km. The probability |$\pi _1^{{\rm{BS}}}$| is estimated either following a uniform distribution, or using a model of smoothed seismicity (set as in Hiemer et al.2014). Both estimations are made considering only the events with M > 4.5 and assigned to BS at LEVEL 2. This means that, in the macro-areas with a subduction interface, the two mentioned methods are combined with the four different buffers adopted at LEVEL 2 (excluding buffer = ∞, since no BS is foreseen in this case), providing 8 alternatives (see details in Appendix A.3.3, Supporting Information).

At LEVEL BS-2, the probability |$\pi _2^{{\rm{BS}}} = {\rm{Pr}}(z|{R_i},{M_j},{C_k} = {\rm{BS}},x,y)$| is studied over a number of discrete depths depending on the magnitude level. In particular, a mean seismogenic depth W_s is assumed equal to the crust width of 27 km, as derived from the CRUST 1.0 model (Laske et al.2013). The number of depth levels is based on the average width of the crustal faults for the given magnitude level M_i, evaluated from Wells & Coppersmith (1994). Then, the depth is divided in n_i equal intervals from 1 km depth (top fault) to a depth of W_s (bottom fault), with probability equal to 1/n_i in each interval. Here a uniform probability distribution is adopted due to the difficulty of constraining reasonably well and homogeneously enough this probability for the offshore earthquakes over the considered domain. More details can be found in Appendix A.3.4 (Supporting Information).

At LEVEL BS-3, the potential combinations of the strike, dip and rake angles are considered. The parameter space is transformed following Selva & Marzocchi (2004), in order to better preserve distances. This transformation provides the strike S between 0 and 180, the dip D between 0 and 180, and the rake R between −180 and 180 (more details in Appendix A.3.5, Supporting Information). Each interval is discretized with 4 strike, 9 dip and 4 rake intervals, for a total of 144 variations for the dPDF. A two-layer Bayesian model is then adopted to assess |$\pi _3^{{\rm{BS}}}$|⁠, assuming a Dirichlet distribution for the prior and a multinomial distribution for the likelihood. This functional choice is a rather common assumption for this kind of assessments (e.g. Gelman et al.1995). At the first layer, a dPDF is obtained at the region level. To obtain this distribution, we set a semi-informative prior based on a priori judgment about the possible combinations of angles (mainly, dip and rake). This prior is then updated considering the two alternative catalogues: AllCMT and EMMA (see Appendix A.1 in the Supporting Information; Pondrelli et al.2011; Ekström et al.2012, and references therein). Both focal planes presented in the catalogues are considered with weight ½. This assigned weight represents the uncertainty on the data (see discussion in Selva & Sandri 2013). At the second layer, the obtained posterior dPDF, again a Dirichlet distribution, is adopted as regional prior for the assessment in each cell (x, y). In case of the presence of known faults in one cell, the regional prior is combined with this specific information available for the cell. Alternative probability distributions are finally sampled from the obtained distribution in each cell, accounting for epistemic uncertainty in their definition. More details can be found in Appendix A.3.5 (Supporting Information).

At LEVEL BS-4, we make a simplification by neglecting the potential aleatory variability on the fault rupture area A, and adopting the best estimate of the scaling laws. We also assume the same scaling laws irrespective of the specific faulting style (strike- or dip-slip). Therefore, we assume that the area is equal to the mean value of the scaling law provided by Wells & Coppersmith (1994), with the best guess aspect ratio.

At LEVEL BS-5, the heterogeneities on the seismic moment distribution are not implemented, nor are those regarding the rigidity μ. This means that a uniform moment release is assumed all over the rupture area A with probability |$\pi _5^{{\rm{BS}}} = 1$|⁠.

At LEVEL IS-1, that is, for subduction earthquakes, the assessment of the rupture area and its centre is made jointly with the following procedure: (a) a grid (x', y') is set over the interface; (b) each point (x', y') corresponds univocally to one volumetric position (x, y, z), that is, for each point (x,y) is fixed at the local depth z of the slab interface; (c) for each magnitude M_j and centre (x', y'), we consider a first guess rectangular area from the scaling laws (Strasser et al.2010); of course, other scaling laws for subduction earthquakes could be considered; (d) several possible geometries are attempted to fill the seismogenic depth and the one that better approximates the initial area is selected. In order to accommodate the area for the given magnitude and the given point (x^′,y^′), we require that the ratio between the selected triangular elements is greater than 25 per cent in the first guess rectangular area. If this condition is not possible the event is discarded. For any given M_j, all acceptable areas for the first guess points (x', y') are assumed equiprobable. The probability is |$\pi _1^{{\rm{IS}}} = \Pr ( {x,y,z,A{\rm{|}}{R_i},{M_j},{C_k} = {\rm{IS}}} ) = \frac{1}{{{n_j}}}$|⁠, where n_j is the number of the accepted areas for magnitude M_j. Two possible alternative models are considered to set the described procedure. One considers the rupture confined within the classical seismogenic depths (here called ‘nucleation zone’), and the other allows rupture propagation to the conditionally stable trench zone (Lay et al.2012; here called ‘propagation zone’, see Appendix A.3.7 in the Supporting Information for additional details). In our model, the larger the magnitude, the more the slip propagation is allowed at shallow depths. This roughly corresponds to the assumption that the rupture energy controls the propagation within the decoupled trench zone.

At LEVEL IS-2, we here only consider uniform slip equal to the mean slip value from scaling laws, with probability |$\pi _2^{{\rm{IS}}} = 1$|⁠, for both alternative models of LEVEL IS-1. More details are reported in Appendix A.3.7 (Supporting Information). Additionally, we point out that this model does not allow either for shallow slip amplification (we only use uniform slip) or for ‘tsunami earthquakes’ with slip confined only in the shallow zone (Polet & Kanamori 2009). Specific models should be developed to account for this kind of events, which is not an easy task, but likely a mandatory one in an operational assessment, as tsunami impact is sensitive to the details of the slip distribution particularly in the near-field of the source.

3.1.2 STEP 2 and 3a: the linear propagation framework and h_max extrapolation

For this application, we model more than 6 × 10⁶ scenarios, arising from STEP 1. Tsunami propagation from source to target is then obtained by linear combination of pre-computed tsunami scenarios. This is a common strategy to significantly reduce the number of numerical simulations needed to represent a very large amount of input parameter combinations. We use a set of about 53 000 Gaussian-shaped elementary tsunami sources covering the whole Mediterranean Sea, representing the initial sea-level elevation (e.g. Miranda et al.2014; Tsushima et al.2014). For each of them, we perform a numerical simulation. We use a nonlinear shallow water multi-GPU code (HySEA, Gonzalez Vida et al.2015) over a 30 arc-sec bathymetric grid (SRTM30+, http://topex.ucsd.edu/WWW_html/srtm30_plus.html) covering the whole Mediterranean Sea. Our database contains, for each elementary source, the tsunami time histories (8 hr sampled every 30 s) stored along the 50 m isobath (about 13 000 points at an average distance of 2 km). The target points considered for our case study are shown in Fig. 3(b) (yellow and red dots).

In order to simulate a given scenario in the ET, we model a fault dislocation with Okada (1985). The Kajura (1963) filter is then used for the sea-bottom/water-surface transfer of the dislocation to obtain the static initial condition for sea level elevation. Then, we perform a rapid estimation of the coefficients that approximate at best the tsunami initial condition by linear superposition of the Gaussian functions. The same coefficients are used to linearly combine the corresponding pre-computed time histories at the 50 m isobath. Eventually, maximum elevations are extracted and extrapolated at the coast by use of the Green's law. A detailed description of this database and its performances is out of the scope of this paper, and it will be reported in Molinari et al. (in preparation).

In this paper, the potential uncertainties on the propagation part are completely neglected, so that a Dirac delta distribution centred at the estimated h_max is assumed for the propagation term of SPTHA (⁠|$\Pr ( {\Psi \ge \psi {\rm{|}}{\sigma _s},\boldsymbol{x}} )$| in eq. (1) is the a step function); Random coupling with tidal phases is totally neglected (e.g. Omira et al.2015). Furthermore, no epistemic alternatives are implemented, to consider for example uncertainties in the bathymetric models, or implementations of dispersion and wave breaking in tsunami numerical models. This is done for purely explicative reasons, since this paper is more focused in exploring uncertainty in the source parameters space.

3.1.3 STEP 4: the ensemble model

The main results of SPTHA are the HCs at the target points, evaluated assuming a Poisson distribution (eqs 2 and 4) with an exposure time ΔT of 50 yr. By combining all the alternative formulations at all the LEVELS of the ET, we obtain a total of 10 752 alternative formulations. Given that no alternatives are considered for the propagation (STEPS 2–3), we have a total of 10 752 alternative formulations for the SPTHA.

This number is very high, and then only a subset, a sample of these alternatives, is considered to obtain the final ensemble model. In particular, we first compute the rank of the alternatives by multiplying the weights relative to each branch of the AT. Then, N_s = 100 alternative configurations are sampled out from the 10 752 alternatives. This number is chosen deliberately small, in order to constrain the computational effort and, at the same time, to enable a sufficient quantification of the order of magnitude of the confidence interval through low and high percentiles (16th and 84th percentiles, in this application). The sampling is performed according to the total weight of each model, so that high-weight models are more likely to be sampled than low-weighted models. Applying eq. (2) to the alternative models, we obtain N_s HCs that represent the input sample for the ensemble model of eq. (7). Since intermediate models between the formulated alternatives are expected to exist, a monomial ensemble distribution is preferred. In particular, we adopt a Beta distribution with the parameters set through the Maximum Likelihood Estimation. This choice is quite common in the hazard literature (e.g. Selva & Sandri 2013; Marzocchi et al.2015), and it is considered sufficient for the explicative nature of this application.

The obtained ensemble models finally quantify simultaneously and consistently the uncertainties (both aleatory and epistemic) on SPTHA.

3.2 Sensitivity tests

Two different kinds of sensitivity tests are here presented. We first analyse the sensitivity of SPTHA results to some of the implemented alternative formulations. Then, we analyse the sensitivity of SPTHA by exploring how much the results are affected by the new combined treatment of IS and BS that we propose, compared to more classical approaches where only IS or BS (including the subduction zones, like in PSHA) were considered.

The sensitivity to alternative implementations is performed considering the N_s models sampled as described in Section 3.1.3. We consider the total mean annual rate curves λ^Tot of eq. (4) at 12 target points along the coastline (Fig. 3b). The results are reported in Fig. 4. We test the sensitivity of SPTHA to (i) two different approaches for the estimation of the b-value at LEVEL 1 of the ET (see Appendix A.3.1, Supporting Information); (ii) two different approaches allowing or not coseismic slip in the ‘propagation zone’ for subductions (see Appendix A.3.7, Supporting Information). We consider this to be likely two of the most important factors in controlling the SPTHA, and then we expect a measurable impact on our results; hence, they are used here primarily for checking the consistency of the proposed approach. A more systematic analysis for all the involved parameters (e.g. Knighton & Bastidas 2015) should be however performed for an operational assessment.

Figure 4.

Sensitivity of the annual rates of exceedance to model implementation at LEVEL 1 (b-value) and LEVEL IS-1 (extent of the shallow rupture in subduction zones) of the ET for the selected points (Fig. 3) along the Italian coastline.

In Fig. 4, the colours of the sampled models are assigned according to the implemented alternative models for the b-value estimation at LEVEL 1 of the ET, and for slip propagation at LEVEL IS-1. As probably expected, a clear separation is observed due to the choice on the b-value; that is, the curves for b = 1 (red and blue curves) feature generally higher values than for the case in which b is estimated from seismicity data. The opposite is observed at some sites for intermediate annual rates (and ARPs), where the black curve (‘propagation’ allowed at shallow depths in the subduction zones, though with b-value estimated from data) crosses through and overcomes the curves for b = 1, for annual rates between ∼10⁻⁴–10⁻⁶. This happens with different details for most of the sites ranging from Bova Marina to Leuca (see also Fig. 3), indicating that the sensitivity to the alternative models for the subduction zone is enhanced in the near-field of the Calabrian Arc.

The second group of sensitivity tests regards the importance of the proposed approach for combining IS and BS seismicity in the model for SPTHA. Compared to the more classical procedure where only subduction events are considered (e.g. González et al.2009), or to cases in which events are considered equiprobable in a volume (Sørensen et al.2012), one of the most innovative part of the proposed methodology is the separation between IS and BS to allow a more focused treatment of the available information. A coherent test has to be developed to evaluate the impact of this procedure on the hazard results. With this goal, we consider the best guess models (highest weights) for all LEVELS but for LEVEL 2, obtaining the following configurations:

1. BS-only: a buffer = 0 is assumed, so that all events are modelled as BS;

2. IS-only: a buffer = ∞ is assumed, so that all events are modelled as IS (in the regions where a subduction interface is not present, no seismicity is modelled);

3. IS+BS: the 12 alternatives of LEVEL 2 are considered to combine IS and BS.

In this way, we end up with 14 models for comparing the result of BS-only (model 1), IS-only (model 2) and IS+BS (models 3–14) within a common and homogeneous framework.

In Fig. 5, we report the results of this comparison for the total mean annual rates of exceedance λ^Tot at the same 12 target points of Fig. 4. The SPTHA of IS+BS combinations results very similar to each-other, demonstrating that the epistemic uncertainty modelled through the 12 alternatives considered in LEVEL 2 has a relatively small impact on the SPTHA results.

Figure 5.

Sensitivity of total mean annual rates of exceedance to the different approaches used in SPTHA. In our application, considering only subduction zones (IS) would lead to hazard underestimation, while treating all sources as background (BS), that is, letting them occur everywhere in the seismogenic volume as in PSHA, leads to an overestimation. For comparison, all the models shown in Fig. 4 are also reported in the background (light grey).

The results of BS-only seem to overestimate the mean annual rates for all the intensity levels, especially for the highest hazard intensities. We report also the results from the N_s models of Fig. 4 (light grey), even if not all those models are completely coherent with BS-only and IS-only cases, which consider only the most likely alternative for all LEVELS but LEVEL 2. The curves of BS-only overcome the N_s models at all the sites but for the lowest intensity levels. Thus, the excessively large aleatory variability modelled in the BS-only has the effect of maximizing the potential tsunami impact at any target site. This effect may be due to the existence of potential earthquakes with high magnitudes for any possible faulting mechanism in the entire regions. For example, the BS-only considers an earthquake with magnitude 9 with a fault centred in the corner of the region at depth W/2, for each possible faulting mechanism, even if there is no fault that may cause such an event. On the opposite, in all the other models magnitude 9 events are forced to be on the subduction. Consequently, in the BS-only model, all events maximizing the tsunami impact at any target site exist, even those that are not physically possible. As a result, a general overestimation of the hazard at all ARPs occurs in BS-only cases.

The SPTHA of IS+BS is instead in general higher than that of IS-only, at least for intermediate and long ARPs. The effect is less pronounced for sites that are closer to the Calabrian Arc subduction zone (Bova Marina to Leuca), where the curves partially overlap. This shows that, at least in complex tectonic environments, crustal seismicity need to be taken into account. It is possible that for high intensities (i.e. for long ARPs) these differences would be to some extent compensated by including heterogeneous slip distributions and shallow slip amplification in subduction zones, which we did not consider in our case study.

In summary, the comparison of IS-only, BS-only and IS+BS shows that: (i) the epistemic uncertainty introduced for separating the IS and BS contributions is significantly smaller than the uncertainty introduced when only BS is considered; (ii) the combination of IS and BS is always important (particularly if the target site is not directly in front of a subduction zone); (iii) given the directivity of the tsunami propagation and its sensitivity to the source geometry, the excessive extension of the aleatory variability may have a very significant impact on the hazard results, as in the case of the BS-only model, potentially beyond the epistemic limits.

3.3 SPTHA results: quantification of the impact on hazard of aleatory and epistemic uncertainty on sources

In Fig. 6, we report the main results of the case study. The ensemble model is shown through its statistics, namely the mean and median HCs, and the 16th–84th percentiles confidence intervals. For comparison, also the N_s input models in the previous section are reported in light grey.

Figure 6.

Mean, median, 16th and 84th percentile hazard curves (in red) from ensemble modelling, showing the exceedance probability in 50 yr as a function of h_max. The grey curves in the background are the hazard curves obtained from the individual models sampled from the ensemble.

The visual analysis of these curves indicates that in most of the cases the ensemble of alternatives is reasonably represented by the statistical description of the ensemble model. However, in some cases the highest hazard curves fall outside the confidence interval identified by the 16th and the 84th percentiles. These curves are likely related to models with a relatively low subjective credibility (weight). Hence, their impact on the ensemble distribution is relatively small.

The spatial distribution of the HCs can be explored also considering either Hazard Maps, which report the intensity corresponding to the selected probability thresholds, or Probability Maps, which report the exceedance probability for the selected intensity thresholds. These maps are obtained by cutting the HCs through either horizontal thresholds (for Hazard Maps) or vertical thresholds (for Probability Maps). The results are reported at different statistics for the epistemic variability: the mean and median values (to represent the central value of the aleatory uncertainty) and the 16–84 percentiles confidence intervals (to represent the epistemic uncertainty).

The resulting Hazard Maps for the threshold corresponding to 1 and 10 per cent in 50 yr (corresponding to 475 yr and 4975 yr of ARP) and the Probability Maps for the thresholds equal to 1 and 5 m of h_max are reported in Figs 7 and 8, respectively. By comparing for example the hazard maps, a correlation between the hazard intensity at the different ARP and the relative position of the target sites respect to the two subduction interfaces can be inferred to be a general stable trend along the target coasts. Indeed, the biggest intensities are observed along the most Southern coast of the Calabrian peninsula and the Eastern coast of Sicily, and on the Southernmost tip of Apulia. This is likely the result of a complex combination of the relative position of the targets and the sources, and of the tsunami propagation. Possibly, some wave focussing occurs at several sites (e.g. Bova Marina, Crotone, Leuca, for which the HCs of Fig. 6 are also more sustained), which was already observed in previous studies (e.g. Tiberti et al.2008; Basili et al.2013).

Figure 7.

Hazard maps for 10 per cent in 50 yr (left column) and 1 per cent in 50 yr (right column) at the different statistics of Fig. 6 of the epistemic uncertainty.

Figure 8.

Probability maps for 1 m (left column) and 5 m (right column) at the different statistics of Fig. 6 of the epistemic uncertainty.

These features, and in particular the relative importance of the different source regions on the SPTHA results can be analysed through the hazard disaggregation (Bazzurro & Cornell 1999) as:

\begin{equation}\Pr \left( {{R_i}{\rm{|\Psi }} \ge {\rm{\psi }}} \right) = \ \frac{{\lambda _i^{{\rm{Tot}}}\left( {{\rm{\Psi }} \ge {\rm{\psi }},{\rm{\ }}\boldsymbol{x}} \right)}}{{\mathop \sum \nolimits_i \lambda _i^{{\rm{Tot}}}\left( {{\rm{\Psi }} \ge {\rm{\psi }},{\rm{\ }}\boldsymbol{x}} \right)}}\end{equation}

(11)

In Fig. 9, we report the results for Ψ = 1 m, where we split the contribution of the sources in three groups: contributions from the tectonic regions including the Hellenic Arc (regions 7, 29 and 35); those including the Calabrian Arc (regions 12 and 26); and those from all the other source regions. We applied also here the concept of ensemble modelling to the probability values obtained through eq. (11), thus quantifying also epistemic uncertainty through its statistics. In Appendix A (Supporting Information), we report also the corresponding figures (Supporting Information Figs A10 and A11) for two further thresholds of Ψ = 3 m and Ψ = 5 m.

Figure 9.

Disaggregation per source region for intensity Ψ > 1 m; the corresponding figures for Ψ > 3 and 5 m are reported in the Supporting Information (Figs A10 and A11). The uncertainty emerging from the ensemble model is represented, for each region, through the mean value (cross), the median value (circle), and 16–84 confidence intervals (error bars), and with different colours and symbols for the Hellenic Arc, the Calabrian Arc, and for the other regions, in order to highlight their relative contributions at each site.

The results show that the Hellenic and Calabrian Arcs regions (the ones containing also the subduction interfaces) tend to dominate the hazard in most of the selected target points, in some cases in combination with crustal sources, particularly for relatively short ARPs; the situation becomes more complex by considering intermediate or long ARPs. However, in all target points, and in particular in the most Eastern sites (from Sibari to Brindisi), it is possible to observe the non-negligible contribution of several purely crustal regions. A large variability on the results of disaggregation is expected if different thresholds Ψ are selected. For example, the relative importance of purely crustal regions increases with higher thresholds, as shown in Appendix A (Supporting Information). It can be noted indeed that the general influence of the regions in the near-field (with only BS, reported in black) tends to increase, as the intensity level increases. At the same time, as the number of sources that may cause larger tsunamis diminishes, the impact of epistemic uncertainty on their rates tends to increase, resulting in larger error bars in Supporting Information Figs A10 and A11. Note also that in several cases, the mean values fall outside the error bars, indicating distributions highly asymmetric with fat tails towards large values.

Another kind of disaggregation is performed by grouping the sources per-cell, instead of per-region as above. An equation similar to eq. (11) is then used, and we obtain the probability that if a SPTHA threshold is exceeded at a given site, this is due to a specific source cell, allowing for a visual analysis of the regional distribution of the potential tsunamigenic sources relevant for a specific site and tsunami intensity, and for a classification of their relative importance. Note that, in each cell, both the contributions of BS and IS are considered, by adding the contribution in each cell of IS scenarios with centres inside each cell. The results are shown in Fig. 10 for Ψ = 1 m. In Appendix A (Supporting Information), we report also the corresponding figures (Supporting Information Figs A12 and A13) for two further thresholds of Ψ = 3 m and Ψ = 5 m. For simplicity, in all these maps only the means of the ensemble models are plotted. Again, for each site, the analysis of different intensity thresholds shows that in some cases the relative importance of near-field and more distant sources, and of crustal and subduction sources can be different, and that bathymetric features controlling tsunami propagation contribute significantly in shaping the probability patterns. It can be also noted that the cells causing such large intensities tend to diminish in number, and to concentrate in the near-field with respect to the targets. In addition, the propagation features are even more marked for the highest intensities (in Appendix A, Supporting Information), since the number of influencing cells diminishes, and the propagation patterns of the remaining sources emerges more clearly.

Figure 10.

Disaggregation per cell for intensity Ψ > 1 m; this map is a spatial distribution of the causative source probability for a tsunami at the site which exceeds the selected threshold. The corresponding figures for Ψ > 3 and 5 m are reported in the Supporting Information (Figs A12, A13). Red crosses indicate the location of the target in each panel.

The impact of the epistemic uncertainty is important also in the definition of either intensity or probability values of interest at single locations. In Fig. 11, we report the explicit quantification of the computed uncertainty by defining the intensity corresponding to an exceedance probability value of 10 per cent in 50 yr (475 ARP) and by defining the exceedance probability corresponding to an intensity value of 5 m for the selected points along the Italian coastline. The results show that the selection of any single value may be misleading, since it is strongly affected by the epistemic uncertainty in the hazard assessment. On the contrary, these distributions are more informative, and also enable for quantitatively testing the hazard results against the observed frequency of tsunamis at specific sites, whenever sufficiently complete catalogues exist (Marzocchi & Jordan 2014). For example, we can count the number of time windows in which at least one exceedance of 5 m has been observed at a given site, and compare the observed frequency of exceedance with the distribution reported in Fig. 11 for the same site. Performing a quantitative statistical test, we can evaluate whether both aleatory variability and epistemic uncertainty have been well captured into the actual implementation of SPTHA.

Figure 11.

Uncertainty on the definition of reference intensity and probability values obtained from horizontal (left panel) and vertical (right panel) cuts of the ensemble hazard curves; on the left, the ensemble distribution of the intensity values for an exceedance probability of 10 per cent in 50 yr; on the right, the ensemble distribution of the exceedance probability for an intensity of 5 m. Bars indicate discretized probability density functions, while dashed lines the cumulative distributions functions. To increase readability, the mode of histograms is normalized to 1.

4 FINAL REMARKS

In this paper we present a general procedure for the SPTHA, organized in four steps:

1. development of the source model by an ET with a hierarchic discretization of the (source) parameters space, developing also alternative implementations of the ET and of its probabilistic formulation to explore the epistemic uncertainty;

2. computation of the tsunami wave propagation for all potential sources by a linear model adopting a kernel of Gaussians, evaluating h_max at the 50 m isobaths;

3. (optional) quantification of the tsunami impact at the shore-line or inland, either by run-up through amplification factors at coastal sites, or by explicit inundation evaluation for a subset of sources selected ex post;

4. integration of the alternative realizations of STEPs 1–3 through ensemble modelling, in order to provide unbiased assessment of both aleatory and epistemic uncertainty.

The proposed ET allows hierarchic discretization of the source space. Each LEVEL considers a discrete number of values that jointly define the tsunamigenic seismic sources. The potentially large aleatory variability in the parameter space is efficiently accounted for through a series of conditional dPDF. The specific structure of the ET allows for quantifying all inherent uncertainties and obtaining accurate and as precise as possible results. More specifically, the proposed ET allows to:

1. fully control the energy release in each region;

2. separate the subduction interface (or in principle any well constrained fault system) from the background seismicity (here crustal sources) by considering two seismicity classes (IS and BS);

3. specifically use a large set of different source of information in each seismic region, spanning from the seismicity and the focal mechanism catalogues, the analytical and the empirical laws and models, to the dataset of known faults;

4. implement a potentially large number of alternatives at all steps of the hazard analysis for the quantification of the epistemic uncertainties.

It is also important to stress that when known faults are present within the cells that are treated as background, this information is used and the probability for the earthquake mechanism occurrence in the cell peaks towards values imposed by the presence of the fault.

Ensemble modelling is adopted to quantify all the uncertainties that arise from different and alternative formulations of the SPTHA, with the only requirement that the set of alternative models represents an unbiased sample of the epistemic uncertainty. The result of the assessment is the ensemble distribution on the HCs. The ensemble distribution represents simultaneously all the uncertainties in the hazard (both aleatory and epistemic), and may eventually be summarized through central tendency estimators (e.g. the mean, to represent the best estimated aleatory variability) and confidence intervals (e.g. 16th–84th percentiles, to represent the overall epistemic uncertainty). In addition, the ensemble allows one to propagating such uncertainty to all the secondary estimates, such as hazard or probability maps for selected intensity or probability of interest, or disaggregation analyses, and to quantitatively test the results against past data.

Noteworthy, this procedure is in principle applicable more in general for PTHA, not only for seismic sources. Of course, the treatment of non-seismic sources in a probabilistic framework, for example the landslides would require a specific implementation at several steps. The proposed method also provides a framework for linking region-wide to site-specific SPTHA.

The presented SPTHA methodology is here applied in the Ionian Sea, a relatively small area in the central-eastern Mediterranean Sea, characterized by two subduction zones and crustal faults in a diverse tectonic setting. We used a great deal of seismic and tectonic data and various earthquake rupture models. Through this application, it is shown that our approach provides an efficient way to fully explore both aleatory and epistemic uncertainty: we explore through the ET more than 6 × 10⁶ alternative tsunami scenarios, and we combine them through more than 1 × 10⁴ alternative probabilistic implementations. To derive this set of alternative implementations, we consider alternative formulations of the ET with alternative assessments of the dPDFs, and alternative configurations of the subduction interfaces and their related parameters.

The results of the application consist of (i) two sensitivity analyses (on some alternative SPTHA formulations, and on the separation between BS and IS); (ii) SPTHA results along the southern Italian coastline for an exposure time of 50 yr, along with disaggregation analyses, and relative uncertainty. These results provide interesting and general insights as summarized below:

1. Constraining the aleatory uncertainty in the IS is of primary importance. The separation of IS and BS allows for a full use of all the available information on main structures (subduction interface geometry and kinematics, in our application), without neglecting the potential contribution of less known faults. Distinguishing between IS and BS allows for a full control of the aleatory variability on sources. Not doing so may lead to significant bias on the hazard. On the one hand, considering only IS may underestimate the hazard, especially for target sites not directly located in front of the subduction interface. On the other hand, considering only BS may lead to a significant overestimation of the hazard.

2. An important mixture exists of both IS and BS at all target sites. In particular, in our case study, a significant contribution of BS to the total hazard is observed for both small and large intensities, especially in the sites not directly sitting in front of subduction interfaces. Region disaggregation shows that the main contribution to the hazard at all target sites comes from the Hellenic Arc, from the Calabrian Arc, and from several other purely crustal regions in the mid-near field. In the most Eastern sites the impact of such purely crustal regions becomes comparable (or even larger, e.g. in Brindisi) with the impact of subduction regions. The relative impact of purely crustal regions seems to increase, as the hazard intensity of interest increases.

3. The full propagation of all the uncertainties allows an explicit estimation of the uncertainty on all secondary results of hazard curves, such as hazard/probability maps, reference intensities, etc. In the same way, all sources of uncertainty may be tracked back based on their impact on the results. Here we provide several examples of sensitivity tests, but based on the proposed models a large number of potential tests are thinkable for future analyses.

4. Cell disaggregation helps in decoding the complexity of this assessment, in which millions of tsunami scenarios and thousands of alternative formulations are considered. In particular, cell disaggregation allows defining in a robust manner the source areas of interest for any given target and any given reference intensity, accounting for both efficiency in the tsunami generation and propagation, and likelihood of events in each area.

We thank William Power, an anonymous Reviewer and the Editor Saskia Goes for their thoughtful insights. We are grateful the EDANYA Research Group at University of Malaga for providing the HySEA code for tsunami simulations. We acknowledge financial support by: Italian Flagship Project RITMARE, EC FP7 ASTARTE (Grant agreement 603839) and STREST (Grant agreement 603389) projects, Italian FIRB-‘Futuro in Ricerca’ project ‘ByMuR’ (Ref. RBFR0880SR), INGV-DPC Agreement, Annex B2.

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this paper:

APPENDIX A: Details in SPTHA assessments

A.1 Basic input data.

A.2 Tectonic framework.

A.3 Probabilistic assessments for the Event Tree in each region.

A.4 Disaggregation results at higher intensity levels.

Table A1. Magnitude of completeness in the EMEC regions extracted for the tectonic regions of Fig. 1 (main text).

Table A2. Regionalization nomenclature for the area in Fig. 3 (main text).

Table A3. Number of depth intervals used for all cells for BS, with the exception of the case BS-only.

Table A4. Number of depth intervals and values used for cells belonging to regions in subduction zones for BS-only assessment.

Figure A1. Cumulative frequency–magnitude curves of each tectonic region using the EMEC catalogue. The completeness magnitude is 4.5 for all regions, while the relative completeness years are indicated for each region in the legend.

Figure A2. Representation (mesh of triangular elements) of the two subduction zones of the Calabrian and Hellenic Arcs. Red dots are the geometrical centres of the earthquake ruptures used to spatially explore the slab interface. The blue and red polygons represent the boundaries of the nucleation and propagation zones, respectively. The colour-shaded polygons are the various crustal regions (BS) overlapping the subduction interfaces and which are considered in various combinations with the relevant subduction for the calculation of earthquake occurrence rates.

Figure A3. The grid adopted at LEVELS B1–B3: it is composed by non-conformal equal-area cells of 25 × 25 km² (see the text, for more details).

Figure A4. Alternative Tree for LEVEL 1. To simplify the graphics, we reported all the branches only in few cases, while the sub-tree structure is substituted with ‘clone’ in all the other cases.

Figure A5. Resulting cumulative frequency size distributions, for the 48 alternative models in two randomly selected macro-regions (the macro-region of the Calabrian arc, containing regions #12 and #26, and region #14). The distributions are compared with the estimation data (light blue, completeness level of 4.5), the data from the same catalogue with a completeness level of 5.5 (dark blue), and with the mean model (red).

Figure A6. Alternative Tree for LEVEL 2. To simplify the graphics, we reported all the branches only in few cases, while the subtree structure is substituted with ‘clone’ in all the other cases.

Figure A7. Resulting frequency of IS events with respect to BS evens, for the 12 alternative models. The functions are compared with the data for the three depth buffers (blue), when data are available, and with the mean model (red). Note that in many cases (for the largest magnitudes) less than 5 events are present, so the plotted points are not significant.

Figure A8. Spatial distribution for the whole region, multiplied by the mean model for mean annual rates of exceedance of M = 4.5, for two alternative models: a uniform distribution in each region (left panel) and smoothed seismicity (right panel), both computed from a depth buffer of 5 km at LEVEL 2. Note that, if less than 50 events are found in the area, a uniform distribution is adopted also in the smoothed seismicity branch.

Figure A9. Mean values of the posterior distribution for the region #7. The results are reported for all the combinations of depth buffer for LEVEL 2 (excluding ∞) and the alternatives at LEVEL BS-3 (4 × 2 = 8 rows, from up to down, we consider buffers of 0, 5 km, 10 km and 15 km with catalogue AllCMT first, then for Emma). In the different columns, we group the results for the different rake intervals (from left to right, centred in −90, 0, 90, 180).

Figure A10. Region disaggregation results for an intensity threshold of 3 m. In the regions in which IS is present, both IS and BS contributions are considered.

Figure A11. Region disaggregation results for an intensity threshold of 5 m. In the regions in which IS is present, both IS and BS contributions are considered.

Figure A12. Cell disaggregation results for an intensity threshold of 3 m. In each cell, both IS and BS events with fault centres in each cell are considered. The red crosses report the locations of the target considered in each panel.

Figure A13. Cell disaggregation results for an intensity threshold of 5 m. In each cell, both IS and BS events with fault centres in each cell are considered. The red crosses report the locations of the target considered in each panel.

(Supplementary Data)

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.