Effect of Shallow Slip Amplification Uncertainty on Probabilistic Tsunami Hazard Analysis in Subduction Zones: Use of Long-Term Balanced Stochastic Slip Models

The complexity of coseismic slip distributions influences the tsunami hazard posed by local and, to a certain extent, distant tsunami sources. Large slip concentrated in shallow patches was observed in recent tsunamigenic earthquakes, possibly due to dynamic amplification near the free surface, variable frictional conditions or other factors. We propose a method for incorporating enhanced shallow slip for subduction earthquakes while preventing systematic slip excess at shallow depths over one or more seismic cycles. The method uses the classic k−2 stochastic slip distributions, augmented by shallow slip amplification. It is necessary for deep events with lower slip to occur more often than shallow ones with amplified slip to balance the long-term cumulative slip. We evaluate the impact of this approach on tsunami hazard in the central and eastern Mediterranean Sea adopting a realistic 3D geometry for three subduction zones, by using it to model ~ 150,000 earthquakes with Mw\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{w}$$\end{document} from 6.0 to 9.0. We combine earthquake rates, depth-dependent slip distributions, tsunami modeling, and epistemic uncertainty through an ensemble modeling technique. We found that the mean hazard curves obtained with our method show enhanced probabilities for larger inundation heights as compared to the curves derived from depth-independent slip distributions. Our approach is completely general and can be applied to any subduction zone in the world.


Introduction
A relatively high rate of great seismic events (M w C 8.0) characterized the last two decades. Most of these events occurred along subduction zones and triggered some of the strongest ever-recorded tsunamis (e.g., 2004 M w 9:2 Sumatra-Andaman, 2010 M w 8:8 Maule, and 2011 M w 9.1 Tohoku). Some of these great earthquakes revealed unprecedented rupture features, for example the Tohoku earthquake that produced an unexpectedly large amount of slip (* 50 m) just at the trench, resulting in a huge tsunami (e.g., Romano et al. 2014;Lorito et al. 2016;Lay 2018). Before this earthquake, it was commonly stated that the accretionary sedimentary wedges could not accumulate sufficient strain to produce a large co-seismic slip (e.g., Hyndman et al. 1997;Moore and Saffer 2001). Even smaller events, such as the 2010 M w 7:8 Mentawai earthquake (classified as a tsunami earthquake, e.g., Yue et al. 2014), produced larger than expected tsunami waves due to relatively large slip at shallow depths. It was also observed that shallow subduction earthquakes tend to have a longer normalized source duration than deeper ones, which was explained by inverse dependence of the rigidity and/or stress drop with depth (Bilek and Lay 1999;Geist and Bilek 2001). Moreover, the depth-dependent frequency radiation recorded during great earthquakes (Wang and Mori 2011;Lay et al. 2012), featuring generally higher-frequency seismic radiation zones at depth, has been interpreted in the framework of geometrical and structural segmentation of the slab and variation of thermal properties with depth (Satriano et al. 2014). Some numerically simulated dynamic effects may indeed favor the up-dip propagation of subduction rupture. These include the bi-material effect (Rubin and Ampuero 2007;Ma and Beroza 2008;Scala et al. 2017), the interaction with radiation reflected and/or converted from the free surface (Nielsen 1998;Lotto et al. 2017;Scala et al. 2019), and the variability of frictional properties (Kozdon and Dunham 2013;Murphy et al. 2018).
All the above observations and their interpretation are relevant for tsunami hazard, as it is well-established that tsunamis and tsunami hazard are sensitive to slip complexity, including shallow slip features, not only in the near-field of the source, but also at regional distances (e.g., Geist 2002;McCloskey et al. 2007;González et al. 2009;Li et al. 2016).
Several strategies have been proposed to produce heterogeneous slip distributions for tsunami hazard calculations. The stochastic slip distributions are either computed from pre-defined statistical distributions (LeVeque et al. 2016;Sepúlveda et al. 2017) or constrained by models of real earthquakes (Mori et al. 2017). These strategies allow for the definition of large slip distribution ensembles to explore the slip distribution uncertainty. This variability in the source models can be propagated in the computation of tsunami (inundation) scenarios and hazard assessment.
Nevertheless, to our knowledge, a shallow slip amplification has not yet been included in tsunami hazard models (Grezio et al. 2017). Some techniques which use sets of data-driven slip distributions (Davies and Griffin 2018;Goda et al. 2014) or pseudo-dynamic kinematic seismic source descriptions (Song and Somerville 2010;Song et al. 2013) may intrinsically address the issue since they are based on seismic or tsunami observations. Moreover, dynamically modified k -2 kinematic source models have been proposed to account for shallow slip amplification and provided enhanced tsunamigenic potential (Murphy et al. 2016). While accounting for enhanced shallow slip, these approaches do not deal with the frequency of occurrence of shallow slip amplification. The enhanced shallow slip should be imposed in a manner such that the time-integrated rate of slip across the whole fault plane reflects the convergence rate and plate coupling over long time periods, that is, among other things, avoiding an unrealistic and unjustified slip accumulation at shallow depths (e.g., Wang and Dixon 2004).
In this work, we propose a methodology to produce sets of kinematic k -2 slip distributions considering slip amplification at shallow depths on the plate interface, while re-balancing the spatial distribution of the consequent long-term slip. For individual earthquakes, we consider systematically enhanced shallow slip. That is, the slip amplitude increases as the rigidity decreases, reaching its maximum at relatively shallow depths, apart from a less strongly coupled zone at the very shallowest part. However, since the long-term spatial slip rate needs to be compatible with the convergence rate and coupling along the interface (e.g., Nalbant et al. 2013), the enhanced slip for a single shallow earthquake should be balanced by an overall lower rate of occurrence for such events. In this way, the earthquakes simulated by this model produce reasonably uniform (hence realistic) long-term slip accumulation, enabling their incorporation in long-term hazard models.
We address the importance of such a model through a long-term seismic-probabilistic tsunami hazard assessment (S-PTHA; hazard for tsunamis generated by earthquakes only). The sensitivity of S-PTHA is demonstrated by using either our method featuring enhanced shallow slip, or the usual depthindependent heterogeneous slip distributions. To perform this sensitivity analysis we use the Mediterranean basin as a case study, this basin contains several subduction structures and has hosted several destructive tsunamis in the past (e.g., the 365 Crete tsunami). We employ a simplified version of the TSUMAPS-NEAM (http://www.tsumaps-neam. eu/) seismicity model. The model is simplified since we consider only subduction events while neglecting all crustal seismicity. In this area, the hazard due to crustal seismicity represents an important component of the total S-PTHA (Sørensen et al. 2012;Selva et al. 2016). One should eventually consider also hazard from non-seismic sources (e.g., Tonini et al. 2011;Grezio et al. 2012Grezio et al. , 2015Urlaub et al. 2018;Paris et al. 2019). As a further simplification in comparison to the TSUMAPS-NEAM model, we consider a narrower range of alternative models to describe the epistemic uncertainty, both as far as the seismicity rates and the inundation models are concerned.
In the following sections, we first introduce our assumptions on the depth-dependence of rigidity and seismic coupling. Then, we synthetically describe the subduction zones in the Mediterranean, which are needed to illustrate the procedure for constructing the seismic slip distributions. Finally, we conclude by explaining the different steps involved in the longterm S-PTHA case study in the Mediterranean Sea.

Seismic Moment and Rigidity/Coupling Variation with Depth
The seismic moment of an earthquake can be defined as: where d is the slip amplitude and l the local rigidity over fault area A, whereas l Á d h i is the spatial average of slip times rigidity.
When the continuous quantity defined by Eq. (1) is modeled on a discretized mesh, the seismic moment can be approximated by the following summation of discrete quantities: where N is the number of mesh cells, whereas l n and A n represent the local rigidity and the surface size of the nth cell, respectively. For a given value of the seismic moment M 0 , the rigidity variations at different locations on the fault constrain the slip value at the same location; for example, the slip is larger in (shallow) low-rigidity areas. We here assume that the rigidity only varies with depth, while we neglect lateral (i.e., alongstrike) variations.
Assuming constant stress drop, to fit the observation that the normalized source duration generally decreases with depth (Bilek and Lay 1999), the rigidity has been proposed to exponentially vary with depth as well, as shown in Fig. 1a, according to: where l is rigidity in GPa, z is depth in km, and from the regression, the parameters a and b take the values of 0:5631 and 0:0437, respectively. A rigidity-depth dependence is also featured by the preliminary reference earth model (PREM, Dziewonski and Anderson 1981;Fig. 1a). However, as pointed out by Geist and Bilek (2001), to match observed tsunami amplitudes generated by several tsunami earthquakes (M w 7:7; 1992 Nicaragua earthquake, Satake 1995; M w 7:5; 1996 Peru earthquake, Heinrich et al. 1998;Tanioka and Satake 1996), the average rigidity of the shallow portion of the subduction zone needs to be somewhat higher than in the Bilek and Lay (1999) end-member case, yet still lower than implied by the PREM. In our analysis, we adopt an intermediate rigidity profile (Fig. 1a), computed as the average between the two aforementioned profiles.
To evaluate the long-term slip accumulation (in multiple seismic cycles), slip on the subduction faults should also account for the aseismic contribution produced by creeping and other non-tsunamigenic, non-radiative episodic slip events (e.g., Scholz 1998), including the contribution of after slip and viscoelastic relaxation phenomena (e.g., Sun and Wang 2015). The balance between the seismic and aseismic slip can be described as controlled by the seismic coupling factor, that is the seismic fraction of total plate convergence in a given time period, from now on called ''coupling'' for simplicity. The total longterm slip D T released at the position x can be modelled as: where d T x ð Þ and d T x ð Þ are the seismic and aseismic contributions to the slip at x, respectively, and k x ð Þ 1 K Á K x ð Þ is the coupling factor, modelled as the product betweenK, the average absolute coupling value of the specific subduction zone, and K x ð Þ, the relative coupling variation at position x. Whatever arbitrary choice we make about the coupling variation, it affects both the slip distribution of single events and the expected cumulative long-term slip. In this work, as a first approximation, a variation of the relative coupling with depth has been imposed: K z ð Þ is assumed to be a constant function along most of the subduction interface, while it monotonically decreases towards the shallowest and deepest boundaries of the seismogenic zone. The decrease at depth is modelled starting from 40 km depth (Tichelaar and Ruff 1993), possibly corresponding to the 350°isotherm, which may be among the factors controlling seismicity cut-off at depth (e.g., Hyndman 2013). Note, however, that the exact down-dip seismogenic limit is relatively unimportant as far as tsunamigenesis is concerned. At shallow depths, the decrease of the coupling occurs where the average dip angle decreases below 10°, generating an almost flat sedimentary wedge (Polonia et al. 2011;Gallais et al. 2012;Maesano et al. 2017), which makes it difficult to define an exact trench location in the specific tectonic context used here as a case study, as discussed in the next section. This depth-dependent coupling function produces a maximum coupling up to the shallowest seismogenic portion of the subduction, possibly allowing for shallow slip amplification. Figure 1b illustrates the function K z ð Þ used in this work, with which we mimic a gradual coupling decrease toward the boundaries of the seismogenic zone.

Subduction in the Mediterranean
The case study focuses on three subduction zones within the Mediterranean basin, namely the Calabrian Arc, the Hellenic Arc, and the Cyprus Arc (Fig. 2).
The Calabrian Arc seismogenic and tsunamigenic potential is much debated; some authors suggested that the Calabrian Arc hosted the tsunamigenic earthquakes of 1693, M7.3, or of 1905, M7.0. For the source of the tsunami associated to the 1693 earthquake, however, there is no consensus on whether it was caused by the subduction interface, by a crustal fault, or by a landslide (Piatanesi and Tinti 1998;Tinti et al. 2001;Gutscher et al. 2006;Gerardi et al. 2008;Argnani et al. 2012). The Hellenic Arc subduction is considered by some authors as the locus of at least two M8 ? earthquakes in 365 and 1303, which generated destructive tsunamis (Guidoboni et al. 1994;Papazachos and Papazachou 1997;Tinti et al. 2005;Lorito et al. 2008;Papadimitriou and Karakostas 2008;Ganas and Parsons 2009;De Martini et al. 2010;Maramai et al. 2014). Some others claim that these events were generated on shallower crustal faults embedded in the accretionary wedge (Shaw et al. 2008;Stiros 2010;England et al. 2015).
Available estimates of the convergence rates in these subduction zones are very variable, and we recall just some of them here. The Calabrian Arc convergence rate is estimated at 5 mm/year by Devoti et al. (2008) and recently at 1.5-1.6 mm/year if creeping, or 2:7 À 3:0 if temporarily locked by Carafa et al. (2018). The Hellenic Arc convergence rate is estimated at 23 mm/year in the Ionian Islands (Hollenstein et al. 2008), 35 mm/year in the western Figure 1 a Rigidity profiles as a function of the depth computed from the end-member case of Bilek and Lay (1999) (blue curve) and interpolating the PREM (Dziewonski and Anderson 1981) values (green line). The average between the two profiles is shown in red. b A priori hypothesis of relative coupling as a function of depth. We impose a homogeneous coupling decreasing toward the upper and deeper limits of the seismogenic depth interval. Symbols in the diagram titles for rigidity and relative coupling are the same as in Eq. (3) and in the description of the Eq. (4), respectively part (Reilinger et al. 2006;Noquet 2012), and 10 mm/ year in the eastern part (Reilinger et al. 2006 (Howell et al. 2017) in the western part, and at 7-8 mm/year (Wdowinski et al. 2006), 5-8 mm/year (Noquet 2012) in the eastern part. The values of the coupling coefficients in these three subduction zones are also highly debated. In the Calabrian Arc two competing interpretations-ranging from partially-locked to unlocked or inactive-were recently proposed (Carafa et al. 2018;Nijholt et al. 2018). In the Hellenic Arc interpretations range from full locking (Ganas and Parsons 2009) to low coupling (Shaw and Jackson 2010;Vernant et al. 2014), but also the presence of important along-strike coupling variations was proposed (Laigle et al. 2004).
Rates used afterward for S-PTHA are inherited from the TSUMAPS-NEAM Project. We refer to the Project documentation for further details. We just list for completeness here that, depending on the alternative TSUMAPS-NEAM modelling strategy, either the SHEEC-EMEC seismic catalog (Stucchi et al.  Davies et al. (2017). The geometry of the three slabs was initially derived from the European Database of Seismogenic Faults (EDSF; Basili et al. 2013) and then modified according to newer data where available. In particular, the Calabrian Arc is replaced by the more recent model by Maesano et al. (2017) which is derived from the interpretation of a dense network of seismic reflection profiles integrated with the analysis of the seismicity distribution with depth. The Hellenic Arc is the same as that in the EDSF, but we verified its consistency with recent works by Sodoudi et al. (2015) and Sachpazi et al. (2016). The Cyprus Arc was slightly modified in consideration of the results of recent works (Bakirci et al. 2012;Salaün et al. 2012;Howell et al. 2017;Sellier et al. 2013a, b) that are based on seismic reflection profiles and tomographic and seismological data and constrain the geometry of the western part of the slab.
According to the classification by Clift and Vannucchi (2004), these three subduction systems are of the accretionary type (e.g., Barbados, Nankai, or Makran), as opposed to erosional type (e.g., Mexico, Tonga, or Kermadec). Due to their relatively slow convergence rate, old age, and presence of a thick sedimentary cover onto the lower plate, they show huge accretionary wedges and no clear evidence of a trench, as it is filled with the accreted deformed sediments. In this configuration, there is a large portion of the shallow part of the lower plate in contact with highly-heterogeneous rocks along a lowangle interface, where one can expect large earthquake ruptures to propagate (e.g., Lallemand et al. 1994;Gutscher and Westbrook 2009). Such configuration is sometimes also considered to favor the occurrence of very shallow slip (e.g., Bilek and Lay 2018, and references therein). Due to the large extension of the accretionary prisms and the gentle dip of the interface (see Fig. 2b) our depth-dependent coupling profile determines that high-amplitude coseismic slip is released just below the sea bottom.
Starting from these subduction geometries, we built 3D triangular meshes with * 15-km element size for all of them, using the Cubit mesh generator (http://cubit.sandia.gov). These geometries and their modeling should be compared to the Slab2 model that was very recently published by Hayes et al. (2018). The three meshes are illustrated in Fig. ESM1 in the Electronic Supplementary Material highlighting the depth variation over the modeled interfaces.

Rupture Areas and Stochastic Slip Distributions
We model tsunamigenic earthquakes in a range from M w ¼ 6:0 to M w ¼ 9:0 subdivided into 18 bins. Some previous sensitivity analysis showed that the lower limit of this range might still have a significant effect on the tsunami hazard depending on some particular local conditions (Selva et al. 2016).
We used the TSUMAPS-NEAM parameterization (http://www.tsumaps-neam.eu/documentation, see also Lorito et al. 2015;Selva et al. 2016) and a similar approach for hazard computation, though we considered fewer epistemic alternatives. Here, for example, we consider only one seismogenic depth interval for each subduction zone, and only one earthquake scaling relation (Strasser et al. 2010). The surface extension A is a priori constrained to the expected value of the respective magnitude-size scaling empirical relation, neglecting any predictive uncertainty regarding the fault length and width. Note that the surface extension is assumed independent of the average depth of the scenario, and therefore compatible with the hypothesis of uniform stress drop (Bilek and Lay 1999).
For each modeled magnitude and each possible position in the seismogenic domain of the mesh, a ruptured surface is built by starting from a geometrical center and then by iteratively adding more and more neighbor mesh cells. This procedure is arrested when the selected area exceeds the expected value from the selected scaling relation. Duplicated surfaces may be generated and they are considered only once. The surfaces whose centroid is farther than 0:1 Á ffiffiffi A p from the initial cell centroid are discarded. This selection inevitably leads to less numerous surface sets as the earthquake magnitude increases. No constraint is imposed on the length L (distance along strike) and width W (distance along dip), and thus the set of ruptures explores a wide range of aspect ratios L=W. It is generally observed that for most of the selected areas L [ W due to the larger extension along strike of the subduction zones.
We verified that, for each magnitude bin, the selected set of surfaces covers the irregular seismogenic surface rather homogeneously, apart from some tapering towards the edges of the seismogenic domain. Figures ESM2(a -b-c) in the Electronic Supplementary Material show the number of events generating slip within each cell of the Calabrian Arc mesh at three magnitude bins.
Within each of the identified rupture areas, five slip distributions are computed to explore the earthquake aleatory variability using a stochastic composite source model (Zheng et al. 1994;Ruiz et al. 2011). This model is based on the random spatial distribution of overlapping circular dislocations of different sizes over the pre-defined slipping surface. These individual dislocations will henceforth be referred to as ''sub-asperities''. The number of sub-asperities for a given size is defined to ensure that the slip spectral amplitude decays as k 2 (where k represents the radial wavenumber). The number of asperities of a given size is given by a power law relationship such that the cumulative distribution of sub-asperities against radius is: where p is a fractal dimension constrained by the imposed seismic moment and stress drop (Zheng et al. 1994). In Eq. (5) the fractal dimension is 2 ensuring the k -2 decay of the slip spectral amplitude. In our model, the sub-asperities have radii ranging between R min $ 5Dx and R max $ 0:35W, with Dx being the average of the mesh-cell linear sizes and W the width as inferred from the selected scaling relation. Each sub-asperity contains an individual slip distribution based on the Eshelby's (1957) circular crack slip function (Ruiz et al. 2011). The distribution of circular sub-asperities over non-planar faults is ensured by the implementation of a multi-lateration scheme that allows for the distance across non-planar surfaces to be accurately calculated (Herrero and Murphy 2018). Once all the sub-asperities have been placed on the fault surface, they are summed together producing a slip distribution that has the expected spectral amplitude k -2 decay.
The location of each sub-asperity is randomly chosen according to a Probability Density Function (PDF). In this approach, this PDF is in turn imposed as a combination of two PDFs. The first PDF is depth-independent and is either a Gaussian or a sum of several Gaussian functions. Both the number of Gaussian functions (from 1 to 4) and their centers are randomly drawn from a uniform distribution. The Gaussian function(s) provide a slight focusing of the slip, allowing exploration of the variability of the size and slip amplitude of the main sub-asperity. The second PDF is based on the distribution of rigidity and coupling with depth. The role of this second PDF is central since it is used to include the shallow slip amplification.
So, to obtain the five slip distributions previously mentioned, the first PDF is calculated five times, for all the rupture areas defined at all available positions for all earthquake magnitudes on each considered subduction zone. The detailed description of the features of the second PDF will be provided in the next sub-section.

Slip Weight Function
The final depth-dependent PDF is built by combining the Gaussian PDF with a Slip Weight Function (SWF), which is a function of rigidity and coupling. The average rigidity profile (Fig. 1a) allows us to define a rigidity value as l n ¼ lð z n Þ; where the subscript n refers to the n-th cell and z n represents the average depth of the n-th cell. Similarly, the coupling associated with each cell can be defined as K n ¼ Kð z n Þ: Figure 3 shows the assumed distributions of rigidity l n (panel a) and coupling K n (panel b) for the Calabrian Arc. For a single earthquake, it is reasonable to expect the slip to be larger where the rigidity is smaller, and the coupling is larger. Therefore, we defined: with SWF n representing the cell-discretized Slip Weight Function and C f is a normalization factor defined such that P N n¼1 SWF n ¼ 1; where N is the total number of cells on the seismogenic portion of the subduction interface. Once a specific rupture Vol. 177, (2020) Effect of Shallow Slip Amplification Uncertainty 1503 surface is extracted, the restricted SWF n is normalized and hence it is the second depth-dependent PDF. The SWF n for the Calabrian Arc is shown in Fig. 3c. In Fig. 4, a scheme for the k -2 slip distribution computation is presented for a M w = 8.6 event on the Calabrian Arc. Figure 4a is an example of random multiple Gaussian PDF extraction, whose features were described in Sect. 2.2. The left-hand side of Fig. 4 summarizes the steps leading to the definition of one of the slip distributions for the case with depth-dependent rigidity and coupling. Hereafter, we refer to the set of slip distributions generated in this way as the ''depth-dependent set''. Figure 4b shows the SWF n defined within the ruptured area. Figure 4c is the normalized product between the random Gaussian PDF (panel a) and the SWF n (panel b). This PDF is used to modulate the distribution of the sub-asperity centers that represents the phase of the k -2 distribution.
For comparison, for each slip distribution in this set, we also compute a corresponding depth-independent k -2 slip map by considering uniform rigidity (l = 33 GPa) and coupling on the fault. In this case, the slip distribution depends only on the Gaussian PDF. Hereafter, we refer to this set of slip distributions generated in this simpler way as the ''depthindependent set''. The right-hand side of Fig. 4 shows that in this case the sub-asperity location is modulated only by the random Gaussian PDF.
Figure 4e, f show the slip distributions computed starting from the two different schemes. For the depth-dependent set, the effect of the variable SWF n included in the k -2 PDF is to enlarge the shallow high-amplitude patch of slip along the strike direction. Moreover, the smaller value of the shallow rigidity with respect to the reference one (i.e., l = 33 GPa) contributes to the amplification of the maximum values within the patch. Since homogeneous coupling is imposed for the depth-independent set, it is worth noting that the slip decrease toward the shallower boundary is only due to the tapering effect of the Eshelby's (1957) slip distributions (Fig. 4f).
For lower magnitudes, due to the smaller rupture area in comparison to the mesh size, it is difficult to define the k -2 sub-asperities distribution properly. In the configuration presented in this work, it is not possible for M w \ 8.5. Hence, for smaller magnitudes, no stochastic selection of slip distribution parameter is performed. Sketch of the steps for the definition of the slip distributions. Lefthand column. a A random Gaussian function is combined with b the SWF n to define c a depth-dependent PDF controlling the location of the sub-asperities over the mesh. Right-hand column. d The PDF coincides with the random Gaussian function due to homogeneous rigidity and coupling. e Sample slip distribution belonging to the depth-dependent set. f Sample slip distribution belonging to the depth-independent set. For the same stochastic slip distribution, the depth-dependent SWF n leads to a wider alongstrike extension of the shallow slip asperity. The absolute maximum slip value is larger in panel e due to the smaller rigidity at shallow depths Hence, for the depth-independent set, a uniform slip is imposed for the smaller magnitudes as: where l is the uniform rigidity, and A is the rupture area. Conversely, for the depth-dependent set, for M w \ 8.5, we compute a normalized seismic momentM o ¼ PÑ n¼1 l n SWF n A n ; whereÑ is the number of rupturing cells. Considering the real seismic moment of the event M 0 , the slip d n within each cell n is estimated as: show examples of a M w ¼ 7:5 earthquake with homogeneous slip distribution for the depth-independent set and a SWF n -derived distribution for depth-dependent set, respectively.

Balancing Slip Probability
We defined two sets of slip distributions (either depth-dependent or depth-independent). To check the cumulative slip over the long term, that is to verify whether there is progressively larger unrealistic slip accumulation at shallow depths over multiple events, we computed the mean slip per earthquaked n , that is: where d ni is the slip (in meters) generated in the nth cell by the i-th distribution for a given magnitude M w . The probability P M w ð Þ is computed from the cumulative tapered Pareto distribution (Kagan et al. 2010, Eq. 2). P Sl i jM w ð Þ represents the conditional probability for the slip distribution Sl i , given the magnitude M w . The meand n is separately computed for the two sets, considering all the sampled magnitude bins and, within each bin, all the N M w slip distributions.
Cumulating this mean slip over a large enough number of earthquakes, we obtain, up to a multiplicative constant, a proxy of the slip rate. For example, multiplyingd n by the product k Á year, with k representing the mean annual rate of the considered events (larger than M w ¼ 6:0 in our case) and year being an arbitrarily large number of years (e.g., the number of years after which we could expect at least one event at the maximum magnitude), a spatial pattern of the released co-seismic slip over the long term can be computed. Therefore, Eq. (7) must provide ad n having a pattern compatible with the a priori hypothesis made on the coupling.
As a first attempt, we assume that, for a given magnitude, all the earthquakes have a uniform probability of occurring anywhere on the fault, that is P Sl i jM w ð Þ¼ 1 N Mw 8i. With this ansatz in Eq. (7), the resulting long-term slip from the depth-independent set turns out to be approximately uniform andd n does not show any particular zone of slip accumulation on the subduction interface. A tapering towards the edges of the seismogenic zone emerges. This tapering is due both to the smaller number of events rupturing close to the boundaries as compared to those rupturing in the middle of the fault, and to the intrinsic tapering of the slip distribution. The Electronic Supplementary Material (Fig. ESM4), showŝ d n , normalized over the multiple n locations considered, from the depth-independent set of slip distributions, for the Calabrian Arc.
When the same ansatz is used for the depthdependent distributions set, the systematic shallow slip amplification generates a spatial concentration of accumulated slip around the area where the SWF n is maximum. This concentration of slip is highlighted in Fig. 5a, where we show the normalizedd n for the Calabrian Arc. Subdividing the seismogenic area into a series of along-strike sections (e.g., the black rectangle in Fig. 5a) we computed n as the mean of the normalizedd n within each section. Figure 5b shows the variability ofd n as a function of the average rigidity of the strike section, l. For relatively small rigidity values ( l\30 GPa), corresponding to the shallower depths, there is a drop ind n l ð Þ. This is caused by the same reasons already discussed for the depth-independent set, but also by the near-trench decreasing coupling in the definition of the SWF n (see Eq. (6)). Hence, the shallow coseismic slip is depleted to a certain extent at the locations where the slip is being partially accommodated aseismically in the less coupled zone. For rigidity larger than 30 GPa, a systematic decrease of thed n l ð Þ is observed that is approximately linear in the semilogarithmic plot of Fig. 5b. The maximumd n l ð Þ in Fig. 5b corresponds to the along-strike section highlighted by a black rectangle in Fig. 5a, b. However, in the zone where a relative coupling K z ð Þ ¼ 1 is imposed, the total long-term slip should We compute these parameters considering only the points at those depths where K z ð Þ ¼ 1 (and b [ 0). The best-fit solution is shown in Fig. 5b. From this regression, we can determine a ''correction factor'' to our initial ansatz, which is applied to make the mean in the Eq. (12) approximately uniform with depth. From simple geometrical considerations, we can define a Figure 5 a Normalized stack of the slip computed from Eq. (7), used as an estimate of the total long-term slip. The domain is subdivided into alongstrike sections. An example of such section is the black rectangle in panel (a). From each of these along-strike sections the mean of thed n is obtained. bd n l ð Þ, that is along-strike mean ofd n , plotted versus the average rigidity within each section. Each point refers then to a different along-strike section and the black rectangle in panel (b) highlights the value corresponding to the section enclosed by the black rectangle in panel (a). The black dashed line is the best-fit solution from which the parameter b is extracted (see text for details). c Here thed n is computed when the P Sl i jM w ð Þof Eq. (9) is used: no clear trend against the depth emerges as an effect of the imposed long-term slip balancing. dd n z ð Þ is plotted for non-balanced (blue dots) and balanced (red dots) long-term stack of the slip. To compare the results with the a priori coupling hypothesis, the function K z l. This latter equation can be re-arranged as: Incorporating this correction into the definition of P Sl i jM w ð Þ, Eq. (8) can be used to re-normalize the conditional probability of each slip distribution as a function of the average rigidity of the scenario itself: where now l i is the average rigidity of the ith scenario defined for the magnitude M w , and N M w is still the number of slip distributions defined for that magnitude. It is straightforward to verify that the discrete distribution of Eq. (9) is normalized and can be regarded as a PDF.
We further observe that this procedure generates a new distribution d n l ð Þ that under-corrects the decreasing trend with a non-zero parameter b, that is, this trend is still present after the correction. This occurs because the regression (black dashed line in Fig. 5b) is based on the local rigidity value of each cell, whereas the balancing can be performed only on a non-local property, that is the probability of occurrence of a particular slip distribution (computed from the average rigidity of the slip distribution). However, if we estimate the angular coefficients b 1 and b 2 for the first two iterations, we have b 2 \b 1 , meaning that the remaining unwanted trend tends to be attenuated. Iterating the procedure and replacing the parameter b in Eq. (9) by b 1 þ b 2 we get a new log-linear behavior having more gently steeping b 3 \b 2 . Finally, imposing a tolerance, a limited number of iterations m is always found such that b m $ 0: Therefore, replacing b in Eq. (9) by P mÀ1 l¼1 b l balances the mean of the slip defined in Eq. (7). Within the presented scheme, and for all the three subduction-zones, it was verified that after two iterations the value of b is reduced by at least an order of magnitude.
The stack computed with Eqs. (7) and (9) and using b ¼ P mÀ1 l¼1 b l is shown in Fig. 5c for the Calabrian Arc. The long-term seismic slip now quite satisfactorily matches the desired coupling. Figure 5d shows the finald H n as a function of the along-strike section average depth z when Eq. (9) is used to define the conditional probability given the magnitude of each slip distribution. For the sake of clarity,d H n is compared with the same quantity plotted in Fig. 5b, but they are now both plotted as a function of depth. Figure 5d evidences that the total balanced long-term slip matches the a priori imposed coupling both at shallower and intermediate depths.

From Slip Distributions to S-PTHA
This section describes how the balanced slip distributions are used for the S-PTHA.

Mean Annual Rates of Tsunami Hazard Intensity
Exceedance at a Point of Interest (POI) The total annual rate of exceedance of a given level of inundation height H 0 at each POI can be now computed as (e.g., Lorito et al. 2015): where P Sl i jM w ð Þis the balanced slip distribution Sl i conditional probability, given the magnitude M w , computed as in Eq. (9); P M w ð Þ is the cumulative tapered Pareto as in Eq. (7); k j is the mean annual rate for earthquakes with M w ! 6; for the j-th subduction zone (j ¼ 1; 2; 3); finally, P POI H [ H 0 jSl i ð Þ ½ is the conditional probability of exceedance of the tsunami intensity threshold H 0 at a given POI.
The mean annual rates k j are inherited from TSUMAPS-NEAM. The TSUMAPS-NEAM model considers epistemic uncertainty and the uncertainty on hazard curves is quantified as an ensemble distribution (Marzocchi et al. 2015;Selva et al. 2016). Here, for the sake of simplicity, we always consider only the mean of the epistemic uncertainty.
is evaluated starting with the computation of each individual slip distribution (see the example in Fig. 6a). The sea-bottom coseismic displacement generated from a slip distribution is 1508 A. Scala et al. Pure Appl. Geophys. computed by using dislocations on triangular subfaults in a homogeneous Poisson's solid half space (Meade 2007). The water column acts as a low-pass filter when the sea-bottom displacement is transferred to the sea-surface. This attenuation is considered by applying a two-dimensional filter of the form 1= cosh kH ð Þ, where k is the wavenumber and H the effective height of the water column (Kajiura 1963). The greater importance of the filtering with respect to other approximations, like the linear combination described below, was shown for example by Løvholt et al. (2012). The sea surface displacement obtained in this way from the slip distribution of Fig. 6a is shown in Fig. 6b.
To produce virtual mareograms at the POIs (which lie approximately on the 50 m isobath), the sea surface elevation is used as the initial condition for pre-computed tsunami Green's functions. The Green's functions are the elementary mareograms produced at the POIs by Gaussian-shaped, of * 4 km standard deviation (* 20 km base width) and spacing * 7 km, elementary sea surface elevations . These mareograms were simulated with Tsunami-HySEA, which is a nonlinear shallow water GPU-optimized and NTHMP benchmarked code (de la Asunción et al. 2013;Macías et al. 2016Macías et al. , 2017. The simulation time was 8 h on a spatial domain enclosing the entire Mediterranean from the Gibraltar Strait (with a small buffer in the Atlantic) to the Eastern Mediterranean including the Aegean and Marmara Seas. The topobathymetry employed is SRTM30 ? , which has a resolution of 30 arc-seconds (* 900 m) and is available at http://topex.ucsd.edu/WWW_html/ srtm30_plus.html. The coefficients for linearly combining the mareograms produced by the elementary sources are those allowing the optimal reconstruction of the initial sea level displacement as a linear combination of Gaussian elementary displacements (Fig. 6c, which is the reconstruction of the displacement in Fig. 6b). This reconstruction is based on the potential energy of the displacement field, following Molinari et al. (2016). An example of the virtual mareogram obtained as a linear combination using the coefficients determined in this way is presented in Fig. 6d. The approximations introduced by this . a A slip distribution on the Calabrian Arc. b The corresponding initial sea level elevation by using the dislocations on triangular subfaults and the low-pass filter. c The reconstructed initial sea level elevation from the linear combination of Gaussian elementary sources d The synthetic mareogram at the POI. e From the analysis of the dominant wave period and polarity and the application of the corresponding local amplification factor, the tsunami log-normal PDF is computed Vol. 177, (2020) Effect of Shallow Slip Amplification Uncertainty technique were addressed by Molinari et al. (2016); in particular, it was noted that the non-linearity of tsunami propagation, in the framework of the present linear combination scheme, did not introduce a significant bias but just some dispersion of the residuals of the reconstructed mareograms. Noting this is important since unwanted dispersion or non-linearities related to the numerical scheme might be potential issues arising from using small elementary sources for approximating large-scale tsunamis, which the superposition may indeed practically fix (e.g., Baba and Cummins, 2005). Here, we illustrate this approximation for two earthquakes on the Calabrian Arc (the same as in Fig. 6 and Fig. ESM5) and on the Hellenic Arc (Fig. ESM6), by showing the residuals of the reconstructed field, and of several mareograms at different POIs. We also point out that dispersion of the residuals of the linear combinations with respect to directly simulated mareograms ) is addressed by an error propagation technique, within the Glimsdal et al. (2019) scheme described here below. Further, and finally, we need to convey simulated offshore tsunami heights into amplified heights at the coastline using approximated amplification factors (Løvholt et al. , 2015Glimsdal et al. 2019). For a specific point on the coastline, the amplified height acts as a proxy for the Maximum Inundation Height (MIH) on the coast beyond. As discussed by Løvholt et al. (2013), the amplification factor method gives exact estimates for the MIH under the special condition of non-breaking plane waves, but this method also assumes incident plane waves and hence neglects local effects such as focusing and refraction. It was recently shown that the amplification factor applied to the offshore value at the POI provides a good and almost unbiased estimator of the median of the whole MIH distribution for a set of onshore transects over a coastline stretch of a few kilometers behind the POI (Glimsdal et al. 2019). This MIH distribution is generally well approximated by a lognormal distribution (see also Davies et al. 2017). Glimsdal et al. (2019) estimated the amplification factors as a function of the offshore wave period and polarity for the TSUMAPS-NEAM POIs in the Mediterranean. They also took into account the coastal bathymetry around each POI. In this study, we use their amplification factors. To estimate the median MIH throughout the coast behind a POI, we take the product of the offshore maximum height at the POI, as estimated from the virtual mareogram, with the specific local amplification factor from Glimsdal et al. (2019). In doing this, we extract for each scenario the period and the polarity of the leading wave at the POI estimated, as illustrated in the example of Fig. 6d, from the reconstructed mareogram. Further examples of maxima, period and polarity extraction, following Glimsdal et al. (2019), are given in Figs. ESM5 and ESM6 in the Electronic Supplementary Material. Finally, following Davies et al. (2017), for a given slip distribution Sl i , we compute the conditional probability from a log-normal distribution having the computed MIH as median and a standard deviation r ¼ 0:3.

S-PTHA
In the classical hypothesis of earthquake occurrence as a Poissonian arrival time process, the probability of at least one exceedance of a threshold level of inundation height H 0 over an exposure time T is given by (e.g., Geist and Parsons 2006): The hazard curves are computed through Eq. (11) for each POI.
In the following S-PTHA examples for the casestudy in the Mediterranean, an exposure time T ¼ 50 year is adopted.

S-PTHA Sensitivity
In this section, we show a sensitivity analysis performed using the case study which considers three subduction zones in the Mediterranean as potential tsunami sources. In Fig. 7, we show the hazard curves at three POIs each one located nearby one of the subduction zones.
We compare the hazard curves obtained from the depth-independent and the depth-dependent slip distribution sets. In this latter case, we present the tsunami hazard obtained by either using or not using 1510 A. Scala et al. Pure Appl. Geophys. the balancing of the total long-term slip presented in Sect. 2.5. The probabilities of exceedance corresponding to the average return periods (ARP) of 500 year (*10% 2 50yr) and 2500 year (*2% 2 50yr) are also highlighted. We observe a repeating pattern in the hazard curves for the three models. This similarity at the three locations could perhaps have been expected since each of the three sites is relatively close and landward of the neighboring subduction zone.
Our model from the balanced depth-dependent set, as compared to the classically-used depth-independent case, features a lower probability of smaller intensities, to which a decreased probability of occurrence of the shallow lower magnitude events may contribute, and exhibits a larger probability for higher intensities, likely due to the shallow slip amplification associated with the largest events rupturing almost everywhere over the subduction interface. The cross-over point between the two hazard curves slightly oscillates between MIH = 0.75-1.2 m; it also always occurs for ARPs shorter than 500 year.
It is also worth noting that the unbalanced model would overestimate the tsunami hazard at all the evaluated ARPs for the Calabrian site, or at least up to the 500 year ARP in the other cases, due to the accumulated shallow slip excess. The balanced and unbalanced models tend to provide more and more similar results for longer ARPs/larger intensities. Again, this is likely because larger magnitude events, producing larger slip and overall larger tsunamis, also feature wider ruptures along-dip, on which the balancing effect is less pronounced.
The S-PTHA sensitivity can also be illustrated by directly comparing the tsunami hazard maps. Tsunami hazard maps are obtained from hazard curves by plotting on a map view the MIH values corresponding to a fixed probability/ARP level. Defining MIH DÀD and MIH DÀI as the balanced depth-dependent and depth-independent MIH at each POI corresponding to a given ARP, respectively, we show in Fig. 8 the difference MIH DÀD À MIH DÀI for the two ARPs of 500 year (Fig. 8a) and 2500 year (Fig. 8b) mentioned above. For the first case, on the south-west coasts of the Peloponnesus, on the western coasts of Libya and Cyprus, and on the Ionian coast of Calabria, the depth-dependent case provides larger MIH estimates compared to the depth-independent case. Elsewhere, we found MIH DÀD \MIH DÀI . However, as the ARP increases (corresponding to smaller probability of exceedance and larger expected maximum inundation heights), the balanced depth-dependent set tends to provide larger MIH estimates over all the coastlines in the vicinity of the subduction zones, as well as relatively far and perpendicularly to the source (ideally along the main tsunami energy propagation direction). As an example, for ARP = 2500 year, we found that the depthindependent slip distributions may lead to MIH underestimation from * 0.5 to * 6 m compared to the depth-dependent distributions case (Fig. 8b). Two further sensitivity tests were performed to address how the S-PTHA depends on the slip distribution features.
The first one is a sensitivity test to the variation of the rigidity profile. As shown in Sect. 2, the earthquake occurrence probability and the slip distributions are constrained by the choice on the rigidity/coupling profiles. Hence, we repeated the analysis for the two extreme rigidity profiles of Fig. 1a (Bilek and Lay 1999;and the PREM model, Dziewonski and Anderson 1981). We found that, even for these end-member cases, the depth-dependent probability of occurrence (Eq. (9)) ensures estimates of the balanced mean slip per earthquaked n similar to the one shown in Fig. 5c. The results of this sensitivity test are summarized in Fig. 9, where the tsunami hazard curves are computed at the same POIs of Fig. 7. For the sake of comparison, the depth-dependent and the depth-independent tsunami hazard curves of Fig. 7 are also plotted in Fig. 9. At all three POIs, the expected overall inverse dependence of the tsunami hazard with the rigidity is obtained. Regardless, these two extreme cases still feature a smaller hazard at lower levels of inundation and a larger hazard for the higher MIH as compared to the depth-independent case.
The second and final sensitivity concerns the minimum magnitude for which the stochastic slip is modelled. To extend stochastic k À2 slip distributions to the events with smaller magnitude (7:9 M w [ 8:6) we reduced the minimum size R min of the stochastic slip asperities from 5Dx to Dx (See Sect. 2.3). It is worth stressing that in all the cases, below the imposed limit magnitude, for the depthdependent distributions a slip value proportional to the SWF n is assigned to each mesh cell (see Eq. (6) and Sect. 2.4), while above this limit the stochastic slip also allows spatial slip heterogeneity (e.g., to have small concentrated patches of large slip). The results of this sensitivity test are shown in Fig. 10 where the hazard curves of Fig. 7 are again also reported. Compared to the original results, a crossover between the two sets is still present, while an overall increase of the hazard occurs. This means that in a real application and for coastlines in the nearfield of the tsunami source, the stochastic slip should be applied as much as possible even at relatively low magnitudes, for which a finer discretization might be necessary.
The last two analyses also evidenced a higher hazard sensitivity at the Cyprus POI with respect to that at the other POIs. This might be due to the shallower minimum seismogenic depth combined with the vicinity of the Cyprus coast to the trench. Such configurations appear to be characterized by a larger model (epistemic) uncertainty, therefore at least a more detailed description (i.e., a finer discretization) of the seismic scenarios is recommended to partly reduce this uncertainty.

Discussion
We proposed a method for S-PTHA that allows for the exploration of the expected natural variability  Fig. 7 for the sensitivity test extending stochastic slip distributions down to M w ¼ 7:9. The hazard curves obtained for ''depth-dependent'' and ''depth-independent'' case (red and blue dashed lines respectively) are compared with the similar cases already shown in the Fig. 7 (red and blue solid lines respectively) Vol. 177, (2020) Effect of Shallow Slip Amplification Uncertainty of the seismic slip on a subduction interface to a significant extent. Future possible improvements may include further aspects that are not fully addressed here, which we nevertheless discuss in this section. For example, here the rupture area was limited quite tightly around the expected (best guess) value from the earthquake scaling relations. Although some along-strike variability was imparted, more systematic sampling of different aspect ratios would be perhaps desirable (e.g., Davies and Griffin 2018). Moreover, we quite subjectively judged the sampling of slip variability (five slip distributions per each rupture area) to sufficiently represents the tsunami hazard of our case study, also given the very dense spatial sampling of the rupture positions, which effectively makes the number of slip samples larger than five. However, this is known to be a challenging issue to manage (e.g., LeVeque et al. 2016;Sepúlveda et al. 2017) and a quantitative hazard convergence testing with respect to the sample dimension could be performed. Note that the proposed approach is in principle suitable for including a larger variability of all the seismic parameters, certainly including rupture size and slip distributions. Moreover, our model considers only the endmember case of depth-dependent rigidity and uniform stress drop, whereas a more realistic model explaining observed earthquake durations should also include a variable stress drop with depth (Bilek and Lay 1999;Saloor and Okal 2018). In fact, the rigidity values implied by the constant stress drop endmember case are too low to explain some tsunami observations (Geist and Bilek 2001). Hence, we employed higher rigidity values than those of Bilek and Lay (1999), which however slightly underpredict the observed durations. To avoid this inconsistency, in a future update, we might then impart a decrease of the stress drop toward the surface, along with an increase of the rupture length, to compensate for the shorter duration, while preserving the rigidity value at a given depth. This addition would imply a modification of the SWF n definition, which would then become: SWF n ¼ CK n = l n L 2 n À Á , with L n depending on the average depth of the cell. The consequent combined effect of a reduced shallow slip amplification and of a narrower rupture aspect ratio would have a complex impact on the resulting tsunami and on the tsunami hazard in the near field, which certainly deserves further studies.
In our model the shallow slip amplifications are all only roughly represented through the rigidity variability as a proxy for the fault conditions in the broader sense. As a consequence, our model generates systematic larger slip where the rigidity decreases. Other improvements should be oriented to quantify and consider shallow slip amplifications due to geometrical, frictional, and structural features, as they emerge from rupture dynamics modelling (Ma and Beroza 2008;Murphy et al. 2016Murphy et al. , 2018Scala et al. 2017Scala et al. , 2019, in a more thorough way. This would address for which tectonic settings and to what extent the modelled seismogenic zones are expected to feature shallow co-seismic slip amplification.
We employed a simplified 1-D coupling model that is considered an acceptable assumption in the absence of specific local coupling models. However, the procedure can be readily extended, possibly including a lateral variation of the coupling, if such a model becomes available. It is also worth noting that the seismic coupling only appears in the definition of the SWF n . Its net effect is to reduce the slip amount at very shallow depths for a single event. Alternatively, the coupling could be included in the definition of the conditional probability of occurrence as: where K z i ð Þ is the coupling computed at the average depth z i of the slipping area. The introduction of this term would ensure that an earthquake breaking a limited portion of the less-coupled zone is more likely than an event slipping only within the lesscoupled zone. We verified through preliminary tests that this approach only slightly modifies the hazard curves without changing the qualitative comparison between depth-independent and depth-dependent curves of Figs. 7,9 and 10.
The seismic slip distributions here presented are based on the k -2 paradigm that is widely used by the seismological community due to its ability to reproduce several macroscale direct observations. However, more efforts are needed to compare our slip distributions with real observations systematically.
Recently, the analysis of Source Time Functions is giving important answers about some macroscopic source properties, such as stress drop, duration, and rupture velocity for subduction events as compared to crustal earthquakes and when different subduction zones are considered . Provided enough resolution, a similar approach could be attempted to look for evidence of macroscopic differences between subduction events occurring at different depths, as already highlighted in terms of moment normalized radiated energy (Newman and Okal 1998;Saloor and Okal 2018). Seismic inversion catalogs (e.g., SRCMOD; USGS and others; Mai and Thingbaijam 2014;Ye et al. 2016) could be used to build ''datadriven'' SWF n as a basis for stochastic slip distributions (Mai and Beroza 2002;Goda et al. 2014). Finally, further consistency testing, such as with mareographic and runup tsunami data (Davies and Griffin 2018), would be also desirable.
We also point out that given the, on the average, limited sea depths in the Mediterranean and the relatively short distances and propagation times, the dispersion may be considered perhaps negligible (e.g. Glimsdal et al. 2013). Hence, the results produced by the Tsunami-HySEA code used here may be considered accurate enough. Otherwise, the dispersion would influence the final hazard results, combining with the effects introduced by our depth-dependent model.

Conclusions
In this work, we proposed a methodology to define stochastic slip distributions for moderate-tolarge magnitude earthquakes in a subduction zone, accounting for possible shallow slip amplification. These sets of events are made compatible with the convergence rate and depth-dependent coupling along the subduction interface. Depth-dependent seismicity features have been already investigated or reviewed as described in several recent papers (see Lay et al. 2012 or discussions in Lay 2018 and references therein). Such features were for example interpreted as controlled by the variability of either geometrical or structural and thermal factors (e.g., Satriano et al. 2014;Bletery et al. 2016) and deserve further investigation.
For illustrative and sensitivity testing purposes, we performed a simplified S-PTHA using the proposed approach for exploring the earthquake slip aleatory variability of three subduction zones in the Mediterranean. The proposed method, however, is completely general and it can thus be applied to any other subduction zone.
The shallow slip amplification is included in the single-event distributions through the definition of a depth-dependent Slip Weight Function, directly proportional to the coupling and inversely proportional to the rigidity variation. The k -2 slip distributions thus obtained are characterized by larger patches of higher amplitude slip as compared to the depth-independent case.
To make these single-event slip distributions compatible with the expected long-term slip reflecting convergence rate and the coupling, a depthdependent probability of occurrence must be defined. Since this probability of occurrence is imposed to increase with increasing average rigidity of the rupture area, it is evident that the largest magnitude seismic events, rupturing almost everywhere on the surface fault, can be considered spatially equiprobable for a fixed magnitude value. Therefore, to balance the systematic slip amplification due to those largest events, the probability of occurrence of relatively smaller events (6:0 M w 8:5) needs to be changed accordingly, with the deeper ones more probable than the shallower ones. By comparing our depth-dependent S-PTHA approach for subduction zone earthquakes to the more common depth-independent one, we found a lower probability for smaller hazard intensities and a higher probability for larger hazard intensities.
Some possible improvements of this approach have been extensively discussed, including consideration of variable stress drop together with more systematic exploration of the earthquake magnitudesize relations, and lateral along-strike seismic coupling variations. Vol. 177, (2020) Effect of Shallow Slip Amplification Uncertainty 1515