Probabilistic Tsunami Hazard Assessment in Meso and Macro Tidal Areas. Application to the Cádiz Bay, Spain

Introduction

The mitigation of marine hazards on threaten coasts is a challenging scientific and engineering topic. Addressing such a challenge requires extensive preparation through the development of hazard and risk assessment methods. Tsunami is one of these marine hazards, for which the assessment methods have significantly progressed over the past 2-3 decades. However, as the tsunami historical catalogs are often inherently incomplete, obtaining a purely empirical hazard assessment remains difficult to achieve (Davies et al., 2018). Alternatively, methods combining data and numerical modeling have been developed (Geist and Lynett, 2014), gaining increasing improvement in the aftermath of the 2004 Indian Ocean tsunami (Grezio et al., 2017).

Existing tsunami hazard assessment methods for earthquake sources are basically developed in an analogous way to the seismic hazard assessment. They include the deterministic tsunami hazard assessment (DTHA) and the probabilistic tsunami hazard assessment (PTHA, or SPTHA from seismic probabilistic tsunami hazard assessment). Prior to the 2004 Indian Ocean tsunami, the DTHA method was the widely used among the tsunami scientific community (e.g., Tinti and Armigliato, 2003). Nevertheless, some few PTHA studies could also be found in the literature (Lin and Tung, 1982; Rikitake and Aida, 1988). The DTHA employs particular source scenarios (i.e., worst-case scenario) to numerically predict and map selected tsunami intensity measures (e.g., maximum wave height, flow depth, current velocity) (Tinti et al., 2005; Lorito et al., 2008). Although the DTHA is considered simple to conduct and provides easily interpreted results, it is very sensitive to the selection of the source scenario (Geist and Lynett, 2014). The PTHA, on the other hand, determines the likelihood of exceeding a tsunami intensity measure (e.g., wave height, runup height) at a particular location within a given time period by integrating the modeled hazard from all considered events (Geist and Parsons, 2006; Power et al., 2007; Grezio et al., 2017). In comparison to DTHA, PTHA is relatively complex and computationally demanding. The complexity of PTHA relies on the treatment of uncertainties that often leads to generate a large number of possible events. Annaka et al. (2007) distinguished two categories of PTHA uncertainties: aleatoric referring to the random nature of tsunami generation and its effects, and epistemic related to the insufficient knowledge on data and modeling accuracy. Covering all sources of uncertainties in PTHA remains highly challenging being the incorporation of uncertainties on the tsunami source one of the main difficulties (Basili et al., 2013; Lorito et al., 2015; Lorito et al., 2015; Sepúlveda et al., 2019; Behrens et al., 2021).

A framework for modern PTHA was first presented in the pioneering work of Geist and Parsons (2006). This framework formed the basis for developing numerous PTHA research works (e.g., Burbidge et al., 2008; Otero, 2008; Taubenböck et al., 2008; Suppasri et al., 2012; Yadav et al., 2013; Horspool et al., 2014; Omira et al., 2015; Shin et al., 2015; El-Hussain et al., 2016; Griffin et al., 2016; Hoechner et al., 2016; Zamora and Babeyko, 2020). 7 years after the 2004 event, the massive 2011 Japan tsunami raised the need to consider the complexity of the earthquake rupture in the PTHA. Subsequently, PTHA studies incorporating heterogeneous slip distribution on the earthquake fault (Li et al., 2016; Sepúlveda et al., 2019) or rigidity depth-dependence (Davies and Griffin, 2019) have been conducted.

Probabilistic tsunami hazard analysis applies to different geographic scales. These scales include global-scale PTHA (Davies et al., 2018), regional-scale PTHA (Thio et al., 2007; Sorensen et al., 2012; Power et al., 2013; Lorito et al., 2015; Omira et al., 2015; Li et al., 2016; Zamora and Babeyko, 2020; Basili et al., 2021), national-scale PTHA (Grezio et al., 2012; Suppasri et al., 2012; Horspool et al., 2014; De Risi and Goda, 2016; El-Hussain et al., 2016; Davies and Griffin, 2019; Kotani et al., 2020), and local-scale PTHA (González et al., 2009; Omira et al., 2016; Volpe et al., 2019). Detailed PTHA of a specific coastal segment community (local-scale) often involves high-resolution inundation modeling and treatment of uncertainties from local effects, such as the tidal stage (Omira et al., 2016). Commonly, the sea level variation and the currents forced by the tidal motion are either disregarded or considered and treated as aleatory uncertainties in the numerical models (González et al., 2010; Omira et al., 2016). In some other cases, tidal elevations are added linearly to the maximum tsunami heights, without considering the non-linear shoaling and topo-bathymetry local effects. However, in zones with a macro and meso tidal regime, the effect of the tides on the probability distribution of tsunami hazard can be very important.

In this study, a PTHA methodology and its application to a costal site of the NE Atlantic, Cádiz Bay in Spain, is presented. Cádiz is characterized by a meso tidal regime (higher than 2 m and lower than 4 m, as defined in Davis et al., 1972) and, therefore, the study considers a combined statistical treatment of both tsunami hazard and sea level variation improving the existing mentioned approaches (sea level add linearly or a constant sea level). As a result, inundation maps expressing the likelihood of exceeding tsunami intensity measures (maximum wave elevation, maximum current speed and maximum wave forces) are derived at given return periods.

Test Site Selection

The November 1st, 1755 Lisbon event is probably the earliest well documented tsunami, with references to all over the Atlantic basin, and one of the largest in European historical times. The tsunami was possibly generated by one of the many faults present in the area of the Gulf of Cádiz (Medialdea et al., 2004; Terrinha et al., 2009; Zitellini et al., 2009). This area is tectonically deformed by a slow (3.5 - 4 mm/year) ESE-WNW oblique convergence between the African and Eurasian lithospheric plates (Argus et al., 1989; Stich et al., 2006), which is translated into a relatively low seismic activity distributed over a wide deformation area. Although there are several seismic sources capable of generating major earthquakes, the slow strain rate causes a low frequency of large earthquakes. The historical record of tsunamis, therefore, is very limited in the area, although there are records of paleotsunamis with similar characteristics to that of 1755 (Luque et al., 2001, 2002; Ruiz et al., 2005; Morales et al., 2008; Baptista and Miranda, 2009).

The Bay of Cádiz is geographically located between the longitudes 6° W and 6° 25′ W and latitudes 36° 20′ N and 36° 40′ N, in the Southwest of Iberian Peninsula (Figure 1). It faces West to the Gulf of Cádiz and is landlocked around its Southwestern, Southern and Eastern margins by the mainland. This area is a natural protection zone with large tidal flats, tidal channels (mean spring tidal range ∼ 3.7 m) and several beaches. Cádiz city constitutes the municipality with the highest population in this area, with approximately 130 000 inhabitants. This city (Figure 1) occupies a small peninsula connected by a tombolo (a sandy isthmus or spit that connects the island to the mainland or to another island) to the Leon Island, where the San Fernando municipality is located. The main economic activities in Cádiz are related to tourism, naval construction and harbor activities. Due to its long and nice beaches, the warm weather and its historic heritage, tourism has increased considerably in the last decades.

FIGURE 1

Figure 1. Geographical localization of the Cádiz Bay and its main municipalities. Pale yellow polygons represent urban areas; red lines are the main roads and the gray line is the railroad.

Probabilistic Methodology

Probabilistic tsunami hazard analysis methods generally follow the PSHA (Probabilistic Seismic Hazard Analysis) method, originally developed by Cornell (1968). The PSHA is a widely used method for assessing the seismic hazard and has been described and extensively reviewed, as stated in various reports and documents (National Research Council, 1988, 1997; Senior Seismic Hazard Analysis Comittee, 1997), and several authors have discussed various aspects of it, such as the distinction between aleatory and epistemic uncertainties, the use of synthetic earthquake catalogs (Ward, 1991, 1996, 2000), the Monte Carlo methods (Savage, 1991, 1992; Cramer et al., 1996; Ebel and Kafka, 1999), and the use of Logic Trees (Coppersmith and Youngs, 1986).

Since the seismicity, and consequently the seismic generation of the tsunami is a stochastic phenomenon, it is necessary to use probabilistic analysis to estimate the hazard of these phenomena and their impact on the coast. The approach we used here is based on Monte Carlo techniques, which comprises the following steps:

• Building a database of numerical simulations of tsunamis for different variables: faults, seismic magnitudes, epicenter locations and sea levels.

• Generating a set of scenarios from synthetic seismic catalogs and tidal conditions based on the probabilistic distributions of the involved variables.

• Mapping the maximum values of the selected tsunami intensity measures to assess the hazard for a given probability.

Methods based on Monte Carlo techniques allow the processing of the work variables and their respective uncertainties in a robust way. The variables are treated by their probability density functions (PDF) and different values can be sampled from the distribution function for each simulation.

In the framework of the European FP6 TRANSFER project¹, we developed a PTHA methodology based on the generation of synthetic seismic catalogs and the incorporation of the sea level variation into a Monte Carlo simulation. In this work, we present the application of this methodology step-by-step to one of the project test-sites, the Bay of Cádiz area in Spain. The main novelty of the work is the combined statistical treatment of both tsunami hazard and tidal variation by incorporating the tidal state into the PTHA. The tide is included using a discrete set of constant sea level elevations, and the tsunami tidal current interaction is neglected.

This paper is structured as follows: first, the characterization of the seismic sources and the statistical analysis are presented. Second, the sea level characterization and statistical analysis are performed. Third, the methodology based on the Monte Carlo approach, which has been applied to obtain the synthetic seismic catalog, is described in detail. Fourth, the numerical database and the methodology to generate the probabilistic thematic maps are presented. Finally, the probabilistic thematic maps for the Bay of Cádiz area are derived from the application of the developed methodology.

Tsunamigenic Sources Characterization

The tsunami origin accounts for a major source of uncertainty in PTHA. The importance of the seismic source characteristics into tsunami hazard has been shown in a number of works (Okal and Synolakis, 2004; Gica et al., 2007; Basili et al., 2013; Knighton and Bastidas, 2015; Davies and Griffin, 2019) and some approximations to deal with them from a probabilistic point of view have been proposed (e.g., Sørensen et al., 2012; Knighton and Bastidas, 2015; Lorito et al., 2015). The incorporation of variable or stochastic slip distribution functions (e.g., Goda et al., 2015; Li et al., 2016), as well as improved 3D geological rupture models (e.g., Yamada et al., 2013), into probabilistic analysis are part of the current development of the PTHA methodology.

Here, the characterization of tsunamigenic sources basically consists of identifying the seismic faults that can affect the studied region, and the probabilistic analysis of their parameters. This analysis is based on the historical earthquake data in the zone and can be divided into the following elements:

(1) Seismic zonation and identification of the main tectonic structures based on historic seismicity and geology. This zoning can include far-field and local sources.

(2) Determination of statistical parameters based on the Gutenberg-Richter Law.

(3) Determination of the geometrical parameters for the main fault ruptures and of the scaling relationships between the seismic moment and the source parameters.

Tectonic Structures

Tectonically, the study area can be defined as a diffuse plate boundary (Sartori et al., 1994), where the deformation is taking place by means of different structures distributed along a broad area. This area comprises the eastern edge of the Gloria and Tydeman transforms (where the deformation is mainly concentrated in these shear corridors), the Gorringe Bank, the Horseshoe Abyssal plain, the Portimao and Guadalquivir banks, and the western termination of the arcuated Gibraltar Arc (Medialdea et al., 2004; Zitellini et al., 2009). This deformation zone is the eastern edge of the Azores – Gibraltar seismic zone, being the present-day boundary between the Eurasian and African plates. Here, the strain is mainly compressive or transpressive, and changes to almost pure strike-slip along the Gloria Fault (Laughton and Whitmarsh, 1974; Borges et al., 2001). The relative motion between the two plates in this area is of the order of 3.5 - 4 mm/year (Argus et al., 1989; Stich et al., 2006).

Most damaging earthquakes and tsunamis that affected the coasts of Portugal, Morocco and Spain were probably generated in the SWIT (SW Iberian Transpressive Domain) zone, being the main active tectonic zone and plate boundary in the area (Zitellini et al., 2009; Martínez-Loriente et al., 2013), including the November 1st, 1755 massive tsunami. In the present study only the tsunamigenic areas in the SWIT region have been considered. This area of tectonic deformation is responsible for the offshore seismicity of Southwestern Iberia (Figure 2). The largest instrumental earthquake was the February 28th, 1969 Horseshoe Ms 8.0 (M_w 7.8) event (Fukao, 1973; Grandin et al., 2007) and the respective tsunami was recorded all over the North Atlantic (Gjevik et al., 1997; Baptista, 2019). Moreover, the SWIT was the place of several moderate size events that have been recorded in the last decades. Focal mechanisms and slip vectors show the coexistence of two different processes of deformation, on one hand the presence of mainly reverse and reverse-oblique faulting with N-S to NNW-SSE compression and on the other hand the occurrence of mainly strike-slip faulting with an associated direction of compression approximately NW-SE (Geissler et al., 2010; Stich et al., 2010). These characteristics are explained by the occurrence of strain partitioning in this area of diffuse deformation (Rosas et al., 2009; Terrinha et al., 2009).

FIGURE 2

Figure 2. Map of seismicity and geological structures used as sources in the present study. GBF, Gorringe Bank Fault; MPF, Marques de Pombal Fault; HSF, Horseshoe Fault; PBF, Portimao Bank Fault; CWF, Cádiz Wedge Fault. Circles show the epicenters of the Spanish National Seismic Network instrumental catalog, with the size proportional to the event magnitude. The colormap shows the SWIM high resolution bathymetry (Zitellini et al., 2009).

As depicted in Figure 2, the seismicity can be grouped in clusters associated with the main geological structures. These structures are the main tsunamigenic sources in the Gulf of Cádiz and they are used in this work as maximum earthquake ruptures. These fault traces are parameterized in Álvarez-Gómez et al. (2008); Omira et al. (2009), Lima et al. (2010) and are described in greater detail in Zitellini et al. (2009). The seismicity associated with these structures is used to define its statistical parameters, basically the Gutenberg-Richter Law, as shown in the next section.

Statistical Analysis of the Seismicity

The Gutenberg-Richter Law relates the earthquake magnitude with its frequency:

$$lo{g}_{10}N=a-b{M}_w(1)$$

where N represents the accumulated frequency of the earthquakes with a magnitude higher or equal to M_w. The a coefficient establishes the seismic activity in the zone, while the parameter b is related with the difference between the frequency of small and big earthquakes. The annual rate of earthquakes with a magnitude higher or equal to M_w (λM_W (N/T), where T refers to the time),

is given by the following equation:

$${\mathrm{\lambda }}_{M_W}=e^{(\mathrm{\alpha }-\mathrm{\beta }{M}_W)}(2)$$

where α = a ln(10) and β = b ln(10).

From the Gutenberg-Richter Law the following distribution function F_M0 (M₀) can be derived for the seismic moment M₀:

$$F_{Mo}\left(M_0\right)={\left[1-{\left(\frac{m_{0l}}{m_{0u}}\right)}^{\mathrm{\beta }}\right]}^{-1}\left[1-{\left(\frac{m_{0l}}{M_0}\right)}^{\mathrm{\beta }}\right]m_{0l}\le M_0\le m_{0u}(3)$$

The associated density function f_M0 (M₀) is the following:

$$f_{Mo}\left(M_0\right)={\left[1-{\left(\frac{m_{0l}}{m_{0u}}\right)}^{\mathrm{\beta }}\right]}^{-1}\mathrm{\beta }\left(\frac{m_{0l}^{\mathrm{\beta }}}{M_0^{\mathrm{\beta }+1}}\right)m_{0l}\le M_0\le m_{0u}(4)$$

where m_0l represents the lower limit of the considered earthquake seismic moment and m_0u represents the upper limit of the considered earthquake seismic moment. For each fault m_0u should be equal to the worst-case scenario seismic moment.

We have selected the time period between 1900 and 2002 in the IGN seismic catalog IGN (2011). In this interval the used magnitude m_b Lg (Mezcua and Martínez-Solares, 1983) is homogeneous. After 2002, the method to obtain the magnitude in the Spanish seismic network changed, and no direct relationship between them has been calculated. This magnitude has been transformed to M_w using the relationship published by Rueda and Mezcua (2002) for Iberian earthquakes:

$$M_w=0.311+0.637m_{b}Lg+0.061m_{b}L{g}^2(5)$$

The number of events increases exponentially with time, denoting the increase in sensitivity of the seismographs. To avoid this effect on the statistical framework, we have calculated the magnitude of completeness evolution through time using Gutenberg-Richter fits with a mobile window of 200 events. We obtained a magnitude of completeness of M_W > 3.2 since 1975, which was tested analyzing the cumulative distribution of events. To compute a synthetic catalog assuming a Poissonian model, all the “non-independent” events must be eliminated. We have used the declustering methodology of Gardner and Knopoff (1974), with time and space windows predefined as function of the event magnitude to eliminate the aftershocks. With the catalog filtered and declustered, we calculated the values a and b by means of a maximum likelihood method. The obtained values were b = 1.05 ± 0.07, and a = 4.31 (Figure 3A). We used the software ZMAP (Wiemer, 2001) to perform these statistical analyses.

FIGURE 3

Figure 3. (A) Gutenberg-Richter fit of the instrumental catalog after declustering of the catalog shown in Figure 2 (see details in the text). Triangles show discrete frequencies while the squares show the cumulative distribution. (B) Plot comparing the Gutenberg-Richter Law fitted with instrumental data (thin line) and the historical data. Squares are the mean annual frequencies, light-gray downward triangles are minimum frequencies and dark-gray upward triangles are maximum frequencies. The fitted Gutenberg-Richter Law to the historical data is shown as a thick line. The intersection between both fits takes place on the magnitude value M_W = 5.28.

Although the fit is quite good, the magnitude interval in which the values are computed is very limited and we have no control over the fit for bigger earthquakes. To evaluate the validity of this relation for tsunamigenic earthquakes, we have estimated the annual frequency of great earthquakes from the historical data. We have used magnitudes estimated from the bibliography (Luque et al., 2001, 2002; Ruiz et al., 2005; Garcia-Orellana et al., 2006; IGN, 2009, 2011) and the minimum, maximum and mean annual frequencies have been computed.

In Figure 3B historical data are plotted over the previously calculated Gutenberg-Richter fit. As can be seen, the distribution of the mean frequencies can be fitted confidently by a Gutenberg-Richter Law (R² = 0.998) with values: a = 1.72, b = 0.56; intersecting the instrumental fit in the magnitude M_W = 5.28. This type of two-branch distribution has been proposed in different environments regarding the distribution of great earthquakes and implying for some authors a characteristic earthquake model (Schwartz and Coppersmith, 1984; Davison and Scholz, 1985; Wesnousky, 1994; Main, 2000; Zielke and Arrowsmith, 2008). In our approximation to the PTHA we are interested in earthquakes capable of generating tsunamis, in big and great earthquakes; and in our case, the second branch better describes the occurrence of the tsunamigenic earthquakes in the Gulf of Cádiz.

Once these parameters are known, the annual rate of earthquakes greater or equal to a given value can be obtained from the Gutenberg-Richter Law. For the SWIT zone we are constraining the minimum seismic moment magnitude for tsunamigenic earthquakes to 6.5. Therefore, earthquakes with smaller magnitudes have not been considered.

Geometrical Parameters and Scaling Relations

To define a set of potential earthquake ruptures with different magnitudes on the same fault, we have obtained first the aspect ratio of the largest rupture scenario in the synthetic catalog for each source and from its characteristics and the scaling relations we have derived those of the smaller events. The worst cases, as well as the fault traces that we used for the Cádiz Bay area, are similar to the ones described in Lima et al. (2010); but accounting for maximum earthquake magnitudes of M_w = 8.6 in the Cadiz accretionary Wedge Fault (Figure 2). For the rest of the events, we have maintained constant the strike, dip and rake of the fault, as well as the length (L) - width (W) relation in order to have enough constraints to define the rest of the parameters from the event magnitude (M_w).

The scalar seismic moment equation (Aki, 1966) relates the size of the earthquake with the elastic properties of the media where the event is taking place, the dimensions of the rupture and the amount of slip on the fault plane:

$$M_0=\mathrm{\mu }S\overline{D}(6)$$

where μ is the shear modulus of the upper lithosphere involved in the fault rupture, S is the area of the rupture that generates the earthquake, and D is the average displacement over the fault surface. In this study we have used a shear modulus of 3 × 10¹⁰ Pa, which is widely used for this area in the literature (Grandin et al., 2007; Gutscher et al., 2009; Omira et al., 2011; Roger et al., 2011), although higher values have also been proposed. To obtain the seismic moment we used the Kanamori (1977) relation:

$$M_w=\frac{2}{3}\left(lo{g}_{10}{M}_o-9.1\right)(7)$$

To obtain the source dimensions from the earthquake magnitude the common practice is the use of an empirical relation adequate for the tectonic environment (Stirling et al., 2013). Some of the most commonly used empirical relationships are mainly developed for earthquakes occurred in continental lithosphere (e.g., Wells and Coppersmith, 1994; Leonard, 2010) or subduction context (Blaser et al., 2010; Strasser et al., 2010; Murotani et al., 2013), while for oceanic lithosphere the empirical relations are scarce and related to normal faulting events (Strasser et al., 2010; Álvarez-Gómez et al., 2012). Scholz (2002) proposed a relation between the length of the fault (L) and the average slip of the rupture (D) as a function of the ratio L/W and the tectonic environment. We have maintained this ratio constant for each fault and assumed an oceanic intraplate environment, which is coherent for a distributed deformation boundary. The relation proposed by Scholz (2002) takes the form:

$$\overline{D}\left(m\right)=6.5x{\mathrm{\hspace{0.33em}10}}^{-5}L(m)(8)$$

The average displacement for the maximum magnitude event is obtained once the length of the fault is defined. The maximum length of the fault is constrained from the sea bottom mapping (Figure 2) while the width is limited by the maximum seismogenic depth. This depth varies from 50 km in the western area of the Gulf of Cádiz, where the oceanic lithosphere is coupled (Stich et al., 2005; Cunha et al., 2010; Geissler et al., 2010), to 25 km on its eastern edge where detachment levels exist at shallower depths in transitional or continental lithosphere (Medialdea et al., 2004; Pereira et al., 2011; Silva et al., 2017).

Once the L/W ratio is defined for each fault, combining eq. 8, L/W = constant for a fault and S = W.L in eq. 6, W can be obtained. All the parameters scales in relation with the source area defined from the earthquake magnitude (eq. 6).

Sea Level Statistical Analysis

In this study the statistical treatment of the tidal variation is incorporated towards a PTHA. Based on the data measured by the tidal gauge closest to the study zone, the sea level probability distribution function is obtained. The net sea level variation is due to atmospheric processes, such as atmospheric pressure changes, the effect of wind and due to the astronomical effects. The former contribution has a random nature, while the astronomical tides are deterministic. Therefore, from the tide gauge data, conducting a harmonic analysis both contributions can be separated and treated as different processes.

The astronomical tide is characterized once the main tidal components are identified and their amplitudes and phases are known. The longest frequency that we are able to consider is the one corresponding to the nodal cycle, which has a period of 18.6 years. With the information derived from the harmonic analysis, 18.6 years of astronomical tide elevations (covering a nodal cycle) with a 1 h resolution are reconstructed.

The 20-year sea level data to perform the statistical analysis has been obtained from a tidal gauge located in Bonanza (Figure 1). The Bonanza tidal gage is part of the REDMAR network from Puertos del Estado, it is located just at the end of harbor breakwater in the mouth of the Guadalquivir river. The Bonanza station is selected for several reasons: (1) it is the closest station to Cadiz (30 km away), it has records of more than 25 years, it is a high-quality station taking measurements of sea level every 1 minute, it is a station included in the Spanish Tsunami Early Warning System in the IGN for tsunami detection, this is why its high temporal sampling in taking sea level data. The other two nearest tide stations are the one in Huelva, 90 km away from Cádiz, located also in a mouth of the Odiel river, just in the directional dike, and the Tarifa station at 85 km located in the Strait of Gibraltar, where it is already affected by the tidal wave of the Mediterranean and the Atlantic; (2) Given the large scale of spatial and temporal variation of the tidal wave, the tide is almost the same in Huelva as in Bonanza, there is only a slight temporal lag; (3) In regards the meteorological tide (storm surge) slight variations in the storm surge regime was obtained between the three gauges; and (4) The effect of Guadalquivir river floods in the sea level is really small, the river is being regulated by upstream dams, the average flow discharge is less than 5 m^3/s, and in floods in rain storms the flow discharge is controlled with maximal 40 m^3/s, being the effect on the sea level very small at the mouth. Therefore, the collected sea level records at the Bonanza station have high quality to characterize the harmonics of the tidal wave and the storm surge regime.

The main components arising from the harmonic analysis, after applying the Tide software, are shown in Table 1.

TABLE 1

Table 1. Main components (amplitude and phase) in Bonanza.

On the other hand, the meteorological tide needs a statistical treatment in order to apply this information in a Monte Carlo method. The distribution function of the meteorological tide can be fitted to different theoretical distribution functions, such as the Gumbel distribution, given by the following expression:

$$F_p=e^{-e^{\left(\frac{S_{mm}-\text{Λ}}{\text{Δ}}\right)}}(9)$$

where F_p represents the exceedance probability of a given meteorological tide S_mm, Λ represents the localization parameter and Δ the scale parameter.

Once the localization and scale parameters of the selected distribution function are obtained, all the information needed in the Monte Carlo simulations is known.

To derive the meteorological tide, the astronomical tide has been reconstructed along the same tidal gauge data series period. Reconstructed data has been subtracted from the total sea level variation measured by the tidal gauge. The residue corresponds to the meteorological tide. The distribution function of the meteorological tide and the fit to a Gumbel distribution function is shown in Figure 4. As can be seen, two different tendencies are apparent in the meteorological distribution function. Because of that, two different Gumbel distribution functions have been fitted, one in the upper part representing the extremal regime and the other in the mean and lower part, representing the mean regime. In Table 2, the values corresponding to the location parameters and to the scale parameters for both regimes are shown.

FIGURE 4

Figure 4. Meteorological tide distribution function and its fit to a Gumbel distribution function. Open circles indicate the empirical distribution function of the meteorological tide and straight lines represent the two fitted Gumbel distribution functions.

TABLE 2

Table 2. Gumbel distribution function, localization, and scale parameters.

Monte Carlo Simulation Procedure

Based on the probabilistic characteristics of the sources, and by applying the Monte Carlo method, synthetic earthquake series are obtained for each of the considered faults. The time series corresponding to each fault are integrated to generate a global tsunamigenic earthquake and sea level catalog for the study region. The length of the time series in the present study was selected to 10000 years. To derive confidence bands, for each thematic map, we performed 1000 Monte Carlo simulations for each fault.

The steps followed in the Monte Carlo simulation for each fault (i) are the following (Figure 5):

FIGURE 5

Figure 5. Summary of the algorithm used to generate the synthetic catalog of tsunamis.

(a) Generation of a random number for the probability of the seismic moment and computation of M₀ (Eq. 3).

(b) Computation of the fault dimensions, rupture width (Wi), length (Li) and Area (Si = WixLi) using the selected scaling relations (Eqs. 6, 7, 8).

(c) Computation of the focal depth (h_f) assuming surface rupture:

$$h_f=\frac{W_i}{2}sin\left(\mathrm{\delta }\right)(10)$$

where δ is the fault dip and W_i the rupture width.

(d) Computation of the epicenter. For a given rupture length Li a random number of locations is generated along the rupture zone, leaving a distance of Li /2 on both sides of the rupture (Figure 6). Once the point along the rupture zone has been generated, and assuming that the rupture reaches the surface, the location of the epicenter is computed with the following equations:

FIGURE 6

Figure 6. Fault parameters. (A) Map view; (B) Cross section B - B′.

$$\begin{array}{c}{X}_{0i}=X_{fi}+\frac{L_i}{2}sin\left(\mathrm{\theta }\right)+\mathrm{\sigma }\left(L_{max}-L_i\right)sin\left(\mathrm{\theta }\right)\hfill \\ Y_{0i}=L_{fi}+\frac{L_i}{2}cos\left(\mathrm{\theta }\right)+\mathrm{\sigma }\left(L_{max}-L_i\right)cos\left(\mathrm{\theta }\right)\hfill \end{array}(11)$$

$$\begin{array}{c}{X}_{ei}=X_{0i}+\frac{W_i}{2}cos\left(\mathrm{\delta }\right)\cdot \text{cos}(-\mathrm{\theta })\hfill \\ Y_{ei}=Y_{0i}+\frac{W_i}{2}cos\left(\mathrm{\delta }\right)\cdot \text{sin}(-\mathrm{\theta })\hfill \end{array}(12)$$

where X_0i and Y_0i represent the longitude and the latitude of the middle point along the fault trace with length L_i; X_ei; and Y_ei represent the longitude and latitude of the epicenter for an earthquake with seismic moment Moi (Figure 6); σ is a random number between 0 and 1 (with uniform distribution); δ is the dip angle and θ is the strike angle and L_max is the maximum rupture length for the fault. Thus, l_e (l_e = L_max - L_i) is the effective distance along the trace (see Figure 6), where the point (X_0i,Y_0i) associated to the epicenter can be located in a random way, for a given rupture length (L_i). If the fault length L_i = L_max then only one point in the fault trace is possible and therefore only one epicenter location.

(e) Computation of the time interval from the last tsunamigenic earthquake event to the next. In the present methodology, earthquakes are assumed to be Poissonian processes; therefore, the probability P(T) of having the next earthquake in a time interval of t or higher is given by the following equation:

$$P\left(T\right)=1-e^{\left(-\mathrm{\lambda }t\right)}(13)$$

Generating a random number for the probability, the time of the interval between two consecutive earthquakes can be obtained. _λ is the annual rate of tsunamigenic earthquakes of magnitude M_Wi or higher. Eq. 13 can be expressed as a function of P(T) to obtain the time t between two consecutive events:

$$t=\frac{\text{ln}(1-P\left(T\right))}{-\mathrm{\lambda }}(14)$$

(f) Generation of a random number for the time along a period of the 18.6 year astronomical sea level series. The astronomical tidal level corresponding to this random moment is considered as the tidal level coinciding with this earthquake event.

(g) Generation of a random number between 0 and 1 for the probability of the storm surge. With this random number, the storm surge S_mm can be derived according to the following expression (see Eq. 9):

$$S_{mm}=\text{Λ}-\text{Δ}\cdot \text{ln}(-\text{ln}(F(x)))(15)$$

where F(x) is the cumulative probability.

(h) Computation of the net tidal level. In this step, the astronomic tidal level and the storm surge are linearly added.

We believe that the assumption of rupture to the surface in steps (c) and (d) concerns magnitudes higher than 6.5 is reasonable taking into account the mechanics of fracture of the schizosphere (the elasto-frictional upper layer of the crust). It is well recognized that major earthquakes nucleate at depth, where the strength of the rocks is higher due to thermo-mechanical conditions (Sibson, 1982). When an earthquake ruptures at depth, overcoming the high frictional strength of the faults, the amount of energy liberated during the earthquake is enough to allow the lateral and vertical propagation of the rupture, reaching the surface and propagating laterally (Scholz, 1998, 2002). Nevertheless, some of the modeled events, especially those with smaller magnitudes, where the surface rupture assumption could be more arguable, do not necessarily reach the surface in nature. Consequently this assumption could be a conservative bias in these cases. Nevertheless, being these events the smaller ones, and as we are using average uniform slip over the fault, this potential bias is reduced and is counterbalanced by the threat posed by the bigger events.

Steps “a” to “h” are repeated for each fault until a simulation of an established number of years is obtained (10000 years in our simulations). Once the synthetic earthquake series are derived for each fault, they all are joined in a single global synthetic catalog.

1000 different 10000-year long tsunamigenic seismic catalogs have been created for each fault. As an example, Figure 7A shows the evolution of one time series (10000-year) of the earthquakes (seismic moment magnitude) considering all the source regions. In Figure 7B, the sea levels corresponding to each of the simulated earthquake scenarios are represented. This figure shows that there are earthquakes that coincide with high sea levels, while others occur during the mean or low tide. These scenarios can create very different impacts on the coastal region.

FIGURE 7

Figure 7. Example of one time series of Moment magnitude (10000 years) for seismic events greater than 6.5 in all source regions (A) and associated sea levels (B).

The distribution of the total number of synthetic events from all sources, in the 1000-year series is shown in Figure 8A; it is worth noting the maximum generated magnitude M_W = 8.55. The epicenter location for each scenario is represented in Figure 8B.

FIGURE 8

Figure 8. (A) Distribution of all synthetic generated magnitudes M_w in discrete (squares) and cumulative(circles) form, and (B) location of the synthetically generated epicenters for all the tsunamigenic earthquake scenarios.

Numerical Database

Once a synthetic catalog is generated (magnitude and epicenter location) based on (1) the synthetic earthquake series (magnitude and epicenter location) for each of the considered faults; and (2) the synthetic sea level series, it is necessary to assign for each synthetic tsunamigenic event in the catalog a flooding scenario on the coast. Thus, a numerical database of scenarios is needed. These scenarios should cover the entire range of possible tsunami events in the synthetic catalog (earthquakes and sea levels). Each tsunami scenario is propagated using a tsunami numerical model, from the source region to the coast (study zone). In this section the numerical model and the simulated scenarios are described.

Numerical Model

We used the C3 (Cantabria-Comcot-Tsunami-Claw) numerical model to simulate the wave generation, propagation and coastal flooding (Olabarrieta et al., 2011). This model integrates finite difference and finite volumes model schemes, and solves the conservative form of the Non-Linear Shallow Water equations (NLSW), considering a Cartesian coordinate system.

These equations constitute a hyperbolic system that can be expressed by the following general conservation equation:

$$\frac{\partial q}{\partial t}+\frac{\partial f}{\partial x}+\frac{\partial g}{\partial y}=s(16)$$

where q represents the vector of the conserved variables, f and g are the numerical flux vectors and s represents the source terms. These vectors are expressed in function of variables such as the water depth H = (h + ξ) and the vertically integrated flow velocities (u, v)

$$\begin{array}{c}q=\left[\begin{array}{c}\hfill H\hfill \\ \hfill uH\hfill \\ \hfill vH\hfill \end{array}\right]\hfill \\ f=\left[\begin{array}{c}\hfill Hu\hfill \\ \hfill H{u}^2+\frac{1}{2}g{H}^2\hfill \\ \hfill Huv\hfill \end{array}\right]\hfill \\ g=\left[\begin{array}{c}\hfill Hv\hfill \\ \hfill Huv\hfill \\ \hfill H{v}^2+\frac{1}{2}g{H}^2\hfill \end{array}\right]\hfill \\ s=\left[\begin{array}{c}\hfill 0\hfill \\ \hfill gh\frac{\partial h}{\partial x}-\frac{{\mathrm{\tau }}_x}{\mathrm{\rho }}\hfill \\ \hfill gh\frac{\partial h}{\partial y}-\frac{{\mathrm{\tau }}_y}{\mathrm{\rho }}\hfill \end{array}\right]\hfill \end{array}(17)$$

where ξ is the free surface elevation; x, y represent both components of the Cartesian coordinate system; g is the gravitational acceleration; ρ is the water density; h is the water depth at rest and τx, τy are the terms that represent the bottom shear stress.

The model solves the differential form of the mass and momentum conservation equations, applying a leap-frog scheme, based on the COMCOT numerical model of Liu et al. (1995). In the inner (detailed) nested grids a finite volume method (George and LeVeque, 2006) is applied to solve the integral form of the conservative equations. The C3 model (Olabarrieta et al., 2011) has been validated in the framework of the European Project TRANSFER using the standards and guidelines proposed by Synolakis et al. (2007).

The datasets used for the bathymetry and topography reconstruction are the Weatherall et al. (2015) database for the regional grids improved with the SWIM database (Zitellini et al., 2009). For the internal, more detailed bathymetric grids, we digitized the nautical charts of the Navy Hydrographic Institute of Spain. The topography and the geographical data used is from the National Geographic Institute of Spain (Figure 9).

FIGURE 9

Figure 9. Topo-bathymetric grids. GEBCO data, SWIM: Zitellini et al. (2009) bathymetric data (The red line is the limit of the high resolution SWIM bathymetry), IHM: Instituto Hidrográfico de la Marina (Navy Hydrographic Institute) data (the black square is the limit of the high resolution local bathymetry IHM), IGN: Instituto Geográfico Nacional (National Geographic Institute) data. Labels in italic font type show the locations commented in the text: R, Rota; RNS, Rota Naval Station; PSM, Puerto de Santa María; V, Valdelagrana; C, Cádiz; CP, Cádiz Port; LCB, La Caleta Beach; SMB, Santa María Beach; LVB, La Victoria Beach; CB, Cortadura Beach; PR, Puerto Real; SF, San Fernando; CH, Chiclana; SP, Sancti Petri.

Three different nested grids have been used to run the numerical model. From the regional bathymetry used to propagate the wave from the source, with a resolution of 864 m, to the detailed grid that includes the topography and local bathymetry with a cell size of 50 m. Between them we used an intermediate grid with a resolution of 216 m.

The submerged profile of the beaches in the area is characterized by the existence of a rocky slab at a water depth of about 5 m and several sandy shoals. In front of the Cádiz urban area there are several elongated shoals and rocky bars (Figure 9). These morphological features play an important role in the local tsunami wave propagation.

Simulated Scenarios

A set of tsunami scenarios are considered for tsunami propagation and for the creation of the tsunami datasets. These cases cover the entire range of possible tsunami scenarios in each zone (potential sources, earthquake magnitudes, and epicenter locations), together with all the possible sea levels.

To establish the tsunamigenic earthquake scenarios in each fault, the relation between the seismic moment magnitude and the rupture area has been taken into account. The width of the fault (W) is constrained by the L/W constant ratio for each fault and the momentum magnitude (eqs. 6-8). The relation between the momentum magnitude and the rupture area for each fault is shown as example for CWF Fault. A sensitivity analysis was carried out to define the minimal area of rupture that generates the earthquake and it was defined as 500 km². Based on a geological analysis of each fault, the maximal rupture area in this case for CWF fault was defined as 15.000 Km². Thus, the definition of moment magnitudes was obtained for areas between 500 km² and 15.000 km² each ΔS = 500 km², with minimal moment magnitude of 6.5 and maximal of 8.4, generating 31 magnitude events. As the magnitude increases, the area also increases exponentially. It is noteworthy that moment magnitudes selected are not equidistant, a change in the magnitude generates a higher variation in the area when the magnitude is large. Therefore, a larger number of scenarios have been considered for the big magnitudes.

On the other hand, in each fault, up to 4 different epicenter locations have been considered for each magnitude. The number of epicenter locations for a specific magnitude event, in a specific fault, depends on the moment magnitude (it defines the source size) and the maximal size of each fault. In the case of smaller magnitude events the size of the source is also small and more epicenter locations can be distributed along the fault. In the case of CWF fault, the largest fault in the area, 4 possible epicenters distributed along the fault were selected for magnitudes between 6.5 and 7.2, 3 for magnitudes 7.3 to 7.6, 2 for magnitudes 7.7 to 8.0 and just 1 epicenter location is possible for magnitudes larger than 8.0. The combination of magnitudes and epicenter locations on CWF fault led to 53 possible earthquake scenarios.

A sensitivity analysis was carried out to evaluate the sea level ranges based on numerical simulations and taking into account the computational time constraints. This analysis allowed considering 4 sea level ranges (Δsl = 1.0 m). In Table 3 the selected tidal levels are indicated where the reference level corresponds to the mean low spring tidal level.

TABLE 3

Table 3. Considered tidal level scenarios (the reference level corresponds to the mean low spring tidal level).

Finally, the combination of magnitudes, epicenter location and sea levels for the CWF fault requires the simulation of 212 scenarios.

A total of 619 scenarios have been modeled (Table 4) using the C3 numerical model. The numerical simulations were performed in the Altamira Supercomputer which is part of the Spanish Supercomputation network. This supercomputer is operated by IFCA (Cantabrian Physics Institute) and it has 512 processors able to manage 4.5 billion of operations per second.

TABLE 4

Table 4. Number of considered scenarios in each fault.

Probabilistic Thematic Hazard Maps

Once the 1000 time series of tsunamigenic earthquake events, coinciding with a given sea level, are created in the Monte Carlo computation, the effect of each tsunami on the coast is obtained using the numerical database generated in the previous steps. An algorithm written in Matlab^® was used for this purpose. For each scenario in the synthetic catalog, the closest case (sea level + earthquake) in the numerical database was identified. For the interpolation of the events of the synthetic catalog, we start from the basis of a 3D matrix of numerically modeled scenarios, with three inputs: magnitude, location of the epicenter and sea level. For a synthetic event of a given magnitude, the closest modeled magnitude is searched in the matrix, for this closest magnitude, the closest modeled epicenter to the synthetic catalog is searched among the modeled epicenters, and finally, the corresponding event with the closest sea level is searched in the matrix. Thus, the event of the synthetic series is associated with the closest simulated event in magnitude, location of the epicenter and sea level. The tsunami intensity measures corresponding to this numerical case (maximum current speed, maximum water depth, maximum wave elevation, maximum Froude number, and maximum force) are assigned to the tsunamigenic earthquake in the synthetic catalog. Total forces, consider hydrostatic and hydrodynamic forces as defined in FEMA (2000), and are associated with a drag coefficient CD, which represents the shear coefficient which depends on the bottom roughness, on the Reynolds number and on the geometry of the structures in the flooding area. CD is based on FEMA (2000); Synolakis (2003); Yeh (2006), which provided recommended shear coefficient values for various obstacle geometries.

The time series of the tsunami intensity measures can be derived for the whole coastal region following the described procedure. Therefore, these time series of the maximum value can be used to perform a statistical analysis of the selected tsunami intensity measures and derive maps with an associated return period.

Once the time series are reconstructed in each high-resolution grid node, these data are sorted in an ascendant order. The cumulative probability is assigned to each organized scenario. This is given by its position, i:

$$P_i=\frac{i}{(m+1)}(18)$$

where i represents the data position and m represents the total number of data. The distribution function of the analyzed parameter vi is given by the following expression:

$$F{Z}_i\left(v\right)=1-e^{-\mathrm{\lambda }(1-P_i)}(19)$$

where FZi (v) is the annual exceedance probability and λ is the tsunamigenic earthquakes annual rate. For the computation of return period T_r the following expression is used:

$$T_r=\frac{1}{FZ(v_r)}(20)$$

where T_r represents the mean number of years in which the variable v_r exceeds a given value.

For each return period and each synthetically generated catalog, a value of the study variable is obtained in all the grid nodes. Therefore, for a given return period, the analyzed variable can be mapped in the high-resolution grid. If several synthetically generated earthquake and sea level catalogs are used, the value of the analyzed parameter for a given return period will present a mean value and a variance. Computing the distribution function of the values obtained for the given return period, the value for a given return period with a given confidence can be obtained.

For a given return period it is possible in each grid point of the detailed mesh (the flooding area) to obtain the desired parameter (e.g., maximum water depth) using the grid point distribution function. The representation on a map of this parameter for all the grid points allow us to produce thematic maps for specific return periods.

In the present study, return periods of 500, 1000, 5000, and 10000 years with a confidence band of 99% have been considered. This confidence band is representative of the mean extreme events. Other confidence bands can easily be derived.

Results

The probabilistic methodology applied to the Cádiz test-site enables the generation of different thematic maps in the study region. These maps do not represent the impact generated by a single tsunamigenic earthquake but the combined effect of different earthquakes in the area. For example, it is possible that an earthquake scenario generates a big impact in La Victoria beach but has a minor impact in Rota. However, another earthquake can generate a big impact in Rota and yet has a small effect in La Victoria. In these maps both Rota and La Victoria will appear to have a great impact, merging the effects of both earthquakes.

The maximum tsunami wave elevation is plotted only in those areas covering the intertidal zones and in those areas above this zone that are flooded due to tsunamis (Figure 10). For a 500-year return period the effects due to tsunamis are very low, and they are restricted to very specific areas. The area showing the highest sea surface elevations corresponds to the harbor of Cádiz and Rota. As expected, as the return period increases, the flooded area also increases, especially in the inner zone of Cádiz Bay as well as in the Sancti Petri zone. In Cádiz city, several beaches such as La Caleta, Santa María del Mar, Victoria, and Cortadura, are affected by the combined effect of sea level (tides) and tsunamis.

FIGURE 10

Figure 10. Maximum wave elevation for different return periods. 99% confidence interval.

In Figure 11, the maximum current speed for different return periods has been plotted. For a return period of 500 years only the zone of Rota presents velocities exceeding the value of 0.5 m/s. For a return period of 1000 years, velocities higher than 4 m/s are achieved in the urban beaches of Cádiz and Rota. For higher return periods most of the coastal area present maximum current speeds higher than 9 m/s.

FIGURE 11

Figure 11. Maximum current speed for different return periods. 99% confidence interval.

Figure 12 the maximum total forces, including the hydrostatic and hydrodynamic contributions with a drag coefficient of 2, have been plotted. For a 500-year return period the net force does not exceed the value of 1 ton/m in the study area. For the 1000-year return period the maximum observed net force increases to 4 tons/m. Rota and the urban beaches in Cádiz (Caleta, Santa María del Mar, Victoria and Cortadura) are the zones with the highest net forces. For higher return periods values of 5 tons/m are reached in large areas. In Rota and in the urban beaches of Cádiz, values of up to 10 tons/m are possible.

FIGURE 12

Figure 12. Maximum forces for different return periods, per unit area of exposed element. 99% confidence interval.

Discussion

Method and Limitations

The PTHA methodology proposed here treats aleatoric uncertainties on the seismic source and tidal stage through a Monte Carlo type framework. It, however, overlooks the epistemic uncertainties such as those related to the modeling of tsunami generation, propagation, and inundation and also those related to the accuracy of the bathymetric-topographic data. Our methodology allows deriving probabilistic coastal hazard maps (wave heights, current speed, and total force) for different return periods. It deals only with tsunamis of tectonic origin and can be applied to any tsunami-threatened coast where high-resolution DEM and records of tidal variations are available.

The PTHA developed here adopts a simple approximation to the source definition fixing some parameters of the earthquake faults (strike, dip, and rake angles) to constrain the variability of the models. This approximation is common in regional PTHA (e.g., Cheung et al., 2011; Heidarzadeh and Kijko, 2011; Power et al., 2012; Davies et al., 2018) but our methodological approximation can be extended to fully incorporate the aleatory uncertainties on the source definition by means of probability density functions using the Monte Carlo approach. Although the incorporation of these branches into the Monte Carlo approach would increase exponentially the number of simulations to compute, the use of GPU-optimized models can overcome this difficulty.

An important source of aleatoric uncertainty in the PTHA concerns the sea level when the tsunami waves reach the nearshore areas. However, this source of uncertainty is commonly disregarded in tsunami hazard assessments. On the other hand, tsunami-tide interaction has been a discussion topic in the literature. Several studies have been carried out to understand whether tide-tsunami interactions enhance/reduce tsunami elevation, flow speed, and arrival times. Based on observations and/or numerical experiments, various authors have found that tsunami inundation in rivers, inlets, harbors, and estuaries is strongly affected by tidal conditions (e.g., Kowalik et al., 2006; Kowalik and Proshutinsky, 2010; Zhang et al., 2011; Lee et al., 2015; Ayca and Lynett, 2016; Shelby et al., 2016; Tolkova, 2016; Ayca et al., 2017). Causes of tide-tsunami interaction are attributed to tidally induced currents and changes in the depth altering the background conditions during the propagation of the tsunami from its source (Weisz and Winter, 2005), both effects are small in the open ocean but increase as the tsunami shoals, mainly in coastal and bathymetric particular configurations as energetic tidal channels communicating large bodies of water or coastal configuration that can induce resonance effects between large scale shallow shelfs, narrow channels, islands, etc. The hydrodynamic pattern of tsunamis occurring under strong tidal regimes can be assessed using different approaches with increasing levels of complexity, depending on the scope of the study and local conditions of the site. Our PTHA approach accounts for the uncertainty of the tide level when the tsunami reaches the coast off Cádiz by incorporating four tidal stages (0.5 m, 1.5 m, 2.5 m, and 3.5 m) into the Monte Carlo simulations. The sensitivity related to these tides levels has not been tested, and it is subject to future studies. The reference for the selected tidal levels corresponds to the mean low spring tidal level. While this approach allows reducing the computational time, it has the limitation of overlooking the non-linear interaction between the tsunami wave, the tidal wave and their induced currents, which could be important on a large spatial and temporal scales. However, in the city of Cádiz that is exposed to the open sea, there are no relevant bathymetry features to generate tsunami-tide interactions, that together with the very large dimensions and periods of the tide (1000 kms and hours) and the tsunami (100 kms and minutes), we assume a small tsunami-tide interaction. Thus, the tide is a boundary condition that can be assumed as instantaneous due to its low speed of variation (in hours), compared to the tsunami speed (in minutes). On the other hand, inside the bay a tsunami-tidal interaction could occur mainly in current speeds that have not been considered here. The study of this tsunami-tide interaction is a challenge that requires further investigation beyond the scope of this work.

Coastal Tsunami Impact

The impact of tsunamis generated in the SW Iberian Transpressive Domain is significant on the Bay of Cádiz and increases, as expected, with the increase of the return period as shown in the results. It is also expected that the incorporation of the sea level variation due to the local tide highly influences the tsunami hazard results on the Bay of Cádiz. To evidence the relevance and the influence of the sea level incorporation, numerical simulations conducted for an example of a large tsunami generated in CWF (Cadiz Wedge Fault) with a M_w = 8.4 in the study area of Cadiz are presented (Figures 13A-C). Figure 13A presents the maximum tsunami wave elevation for a simulation in a low tide sea level. Figure 13B shows the result of the same scenario but incorporating the sea level to the propagation (spring tide 3.5 m, characteristic of the area). And finally, Figure 13C depicts the resulting flooded area for the same scenario executed in low tide condition, and then adding linearly the tide sea level (3.5 m) after the numerical simulation, which is a common approximation to include the tidal effect in PTHA analysis. Two main remarks can be drawn from these results. First, the incorporation of the tide is relevant because the flooded area is larger, and second, this incorporation was made using different initial levels, not just adding tidal levels linearly. Thus, in order to reduce the uncertainty of the sea level in meso and macro tidal areas in a probabilistic analysis, it is required to include the sea level as another variable taking into account its probability of occurrence.

FIGURE 13

Figure 13. Importance of using tides in the simulations. (A) low tide propagation. (B) high tide propagation. (C) low tide propagation and tide added linearly.

Another important issue concerns the use of deterministic and probabilistic thematic tsunami hazard maps. Whether deterministic (conservative scenario) or probabilistic analysis serves better increasing the tsunami resilience of coastal areas depends on the scale and scope of the study. Figure 14 depicts a comparison of the aggregated results (envelope of the worst credible cases in each fault) for the high tide level obtained in a deterministic analysis with those obtained for the probabilistic results for a return period of 10000 years (confidence interval 99%). The comparison clearly indicates that the deterministic approach is very conservative. It is also noticeable that the aggregated results from the worst-case scenarios for each source of the deterministic approach (parameterized in Álvarez-Gómez et al., 2008; Omira et al., 2009; Lima et al., 2010) is associated with a return period longer than (Tr = 10000 years), as expected. Therefore, deterministic maps could be used to define those zones with no probability of being affected by a tsunami.

FIGURE 14

Figure 14. Comparison of the aggregated maximum wave elevation obtained for the high tide condition in the deterministic and in the probabilistic approaches (Tr = 10000 years and a confidence interval of 99%).

Both approaches can address different objectives in the analysis of the affected zone in a tsunami–prone area. Deterministic maps are more subjective but they are useful to establish inundation limits in areas where special infrastructures (nuclear plants, inflammable materials, pipelines) could be built or to assign the evacuation zones. On the other hand, the probabilistic maps are more valuable for vulnerability and risk assessments as well as for insurance purposes, as they provide a more complete picture of the hazard for a given zone with an associated return period. Additionally, building and infrastructure engineering design include probabilistic risk scenarios and risk return periods in the design.

To summarize, despite the limitations of the proposed PTHA methodology in dealing with limited sources of aleatoric uncertainty and disregarding the epistemic uncertainties it interestingly allows combining uncertainties in the generation parameters with the those of sea level in meso and maro tidal areas and performs high-resolution probabilistic tsunami hazard assessment. We believe that the method could highly benefit from further improvements. Particularly, the implementation in a specific case can be improved, with better probabilistic description of the different involved variables (definition of sources, geometrical parameters and scaling relations, Mw distribution, epicenter location in each fault, sea level distributions), and also, with higher resolution of the number of numerically modeled scenarios and higher grid resolutions inland taking into account the different source and sea level parameters, using more powerful and faster models using for example GPU techniques. However, the proposed methodology allows obtaining a reduction of the uncertainties for probabilistic flooding maps, achieving a distribution function of different tsunami intensity measures (max. water depth, max. wave elevation, max. current sped, max. forces, etc.) in each nodal grid point inland.

Conclusion

Deterministic or scenario-based approaches are based on tsunami events that could happen in the future, but have occurrence probabilities that are not well understood. In most deterministic assessments the tsunami scenario represents the worst tsunami case that can take place in the study region, often inferred from historical data. The most difficult task in this approach is to assign a return period to the considered scenario, which sometimes is impossible because of the time span of the seismic catalogs. Nevertheless, the results obtained from scenario-based approaches have been used for several applications, such as the establishment of inundation zones, evacuation maps (Geist and Parsons, 2006) and urban planning.

In this study, a new probabilistic method for the tsunami hazard estimation has been presented and its application to the Cádiz area has been described. This probabilistic method is combined with the application of high-resolution numerical simulations for tsunami propagation and inundation (in the present case the C3 model, Olabarrieta et al., 2011) to derive inundation and hazard maps for given return period and confidence interval. A novel characteristic of the present methodology, based on Monte Carlo simulations, is the statistical treatment of both tsunami hazard and sea level variation, including both astronomical and meteorological effects. Therefore, the present method does not assume that the tsunamigenic earthquake occurs at a given tidal level condition.

The application of the PTHA method described here requires a priori zoning of the source region and the identification of the different faults in each source region. In the specific case of the Cádiz area, only the tsunamigenic areas in the SWIT region were considered. The main tsunamigenic sources in the Gulf of Cádiz were identified as clusters of seismicity in the SWIT region, which represent the maximum earthquake ruptures. GBF (Gorringe Bank Fault), HSF (Horseshoe Fault), MPF (Marques de Pombal Fault), PBF (Portimao Bank Fault), and CWF (Cádiz Wedge Fault) were the five main earthquake ruptures identified in the analysis region. The a and b parameters of the Gutenberg-Ritcher Law for the SWIT region were obtained from the study of the instrumental and historical seismic catalog of the IGN. The tidal level probability distribution function and the characteristics of the main tidal harmonic constituents were obtained from a tidal gauge located near Cádiz. 1000 different time series of combined tsunamigenic earthquake and tidal levels were synthetically generated using the Monte Carlo technique. Each time series had a 10000-year duration.

The C3 model was used to propagate the tsunami from the source region to the coast. The tsunami characteristics were analyzed statistically to derive different thematic maps for the return periods of 500, 1000, 5000, and 10000 years. The derived thematic maps include the maximum wave elevation, the maximum current velocity, the maximum Froude number, and the maximum total forces. The thematic maps obtained for a return period of 10000-year (confidence interval 99%) were compared to those derived from the aggregated deterministic worst-case scenarios for the Cádiz region. This deterministic scenario assumed that the worst-case scenarios for each considered source coincided with the maximum water level measured by the tidal gauge, which means that the probability of occurrence of this extreme event is very low. Thus, the deterministic approach was very conservative. Deterministic maps could be used to define those zones with a very low probability or with no probability of being affected by a tsunami. This is a useful information for visualizing the hazard and the establishment of the maximum limits of flooding, which is necessary for planning the evacuation zones and for defining those areas in which critical infrastructures (e.g., nuclear plants, flammable materials, pipelines) can be built. Ideally, deterministic approach should be linked to the probabilistic study. The probabilistic maps are more valuable for vulnerability and risk assessment and also for insurance proposes.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

MG conceived the manuscript, did the probabilistic approximation, and managed the process. JÁ-G conducted the characterization of tsunamigenic sources. ÍA-Q computed the Numerical simulations. LO and MO compiled the Montecarlo procedure and created the database. RJ, EK, and JB provided some data and internal discussion on the probabilistic model and revised the manuscript. MB, RO, and MM wrote parts of the manuscript and participated in internal discussion. ÍA-A and AL reviewed the probabilistic model. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor declared a past co-authorship with several of the authors MB and RO.

Funding

This research was part of the TRANSFER project (http://www.transferproject.eu/). AL acknowledges the support of the Spanish Grant MTM2008-00759. Some of the figures have been produced with the software GMT (Wessel and Smith, 1998). This is a contribution from the Environmental Hydraulics Institute (Universidad de Cantabria) and the Active Tectonics, Paleoseismology, and Associated Hazards research group (UCM-910368).

Acknowledgments

We thank to the Environmental Ministry of Spain (Dirección General de Costas de Andalucia) for providing us the local bathymetry. We also thank to the Spanish Supercomputing Network the use of the Altamira IFCA supercomputer. We would like to thank to the TSUSY Project (IHCantabria).

Footnotes

1. ^ http://www.transferproject.eu/

References

Aki, K. (1966). Generation and propagation of G waves from Niigata earthquake of June 16, 1964. Part 2. Estimation of earthquake moment, released energy, and stress-strain drop from the G wave spectrum. Bull. Earthq. Res. Inst. 44, 73–88.

Google Scholar

Álvarez-Gómez, J. A., Gutiérrez, O. Q. G., and Aniel-Quiroga, Í, and González, M. (2012). Tsunamigenic potential of outer-rise normal faults at the Middle America trench in Central America. Tectonophysics 574, 133–143. doi: 10.1016/j.tecto.2012.08.014

CrossRef Full Text | Google Scholar

Álvarez-Gómez, J. A., Miranda, J. M., Carreño, E., Matías, L., Martínez-Solares, J. M., and Baptista, M. A. (2008). Review of tsunamigenic faults in the Gulf of Cádiz. Geotemas 10, 1425–1428.

Google Scholar

Annaka, T., Satake, K., Sakakiyama, T., Yanagisawa, K., and Shuto, N. (2007). Logic-tree approach for probabilistic tsunami hazard analysis and its applications to the japanese coasts. Pure Appl. Geophys. 164, 577–592. doi: 10.1007/s00024-006-0174-3

CrossRef Full Text | Google Scholar

Argus, D. F., Gordon, R. G., DeMets, C., and Stein, S. (1989). Closure of the Africa-Eurasia-North America plate motion circuit and tecton- ics of the gloria fault. J. Geophys. Res. 94, 5585–5602. doi: 10.1029/JB094iB05p05585

CrossRef Full Text | Google Scholar

Ayca, A., and Lynett, P. J. (2016). Effect of tides and source location on nearshore tsunami-induced currents. J. Geophys. Res. Oceans 121, 8807–8820. doi: 10.1002/2016jc012435

CrossRef Full Text | Google Scholar

Ayca, A., Lynett, P. J., and Wilson, R. (2017). Effects of tide and wave directionality on localized tsunami-induced currents in port and harbors. Coast. Eng. Proc. 1:8. doi: 10.9753/icce.v35.currents.8

CrossRef Full Text | Google Scholar

Baptista, M. A. (2019). Tsunamis along the azores gibraltar plate boundary. Pure Appl. Geophys. 177, 1713–1724. doi: 10.1007/s00024-019-02344-8

CrossRef Full Text | Google Scholar

Baptista, M. A., and Miranda, J. M. (2009). Revision of the Portuguese catalog of tsunamis. Nat. Haz. Earth Syst. Sci. 9, 25–42. doi: 10.5194/nhess-9-25-2009

CrossRef Full Text | Google Scholar

Basili, R., Brizuela, B., Herrero, A., Iqbal, S., Lorito, S., Maesano, F. E., et al. (2021). The Making of the NEAM Tsunami Hazard Model 2018 (NEAMTHM18). Front. Earth Sci. 8:616594. doi: 10.3389/feart.2020.616594

CrossRef Full Text | Google Scholar

Basili, R., Tiberti, M. M., Kastelic, V., Romano, F., Piatanesi, A., Selva, J., et al. (2013). Integrating geologic fault data into tsunami hazard studies. Nat. Haz. Earth Syst. Sci. 13, 1025–1050. doi: 10.5194/nhess-13-1025-2013

CrossRef Full Text | Google Scholar

Behrens, J., Løvholt, F., Jalayer, F., Lorito, S., and Salgado-Gálvez, M. (2021). Probabilistic tsunami hazard and risk analysis -a review of research gaps. Front. Earth Sci. doi: 10.3389/feart.2021.628772

CrossRef Full Text | Google Scholar

Blaser, L., Krüger, F., Ohrnberger, M., and Scherbaum, F. (2010). Scaling relations of earthquake source parameter estimates with special focus on subduction environment. Bull. Seismol. Soc. Am. 100, 2914–2926. doi: 10.1785/0120100111

CrossRef Full Text | Google Scholar

Borges, J. F., Fitas, A. J., Bezzeghoud, M., and Teves-Costa, P. (2001). Seismotectonics of Portugal and its adjacent Atlantic area. Tectonophysics 331, 373–387. doi: 10.1016/s0040-1951(00)00291-2

CrossRef Full Text | Google Scholar

Burbidge, D., Cummins, P. R., Mleczko, R., and Thio, H. K. (2008). A probabilistic tsunami hazard assessment for Western Australia. Pure Appl. Geophys. 165, 2059–2088. doi: 10.1007/978-3-0346-0057-6_6

CrossRef Full Text | Google Scholar

Cheung, K. F., Wei, Y., Yamazaki, Y., and Yim, S. C. (2011). Modeling of 500-year tsunamis for probabilistic design of coastal infrastructure in the Pacific Northwest. Coast. Eng. 58, 970–985. doi: 10.1016/j.coastaleng.2011.05.003

CrossRef Full Text | Google Scholar

Coppersmith, K. J., and Youngs, R. R. (1986). “Capturing uncertainty in probabilistic seismic hazard assessments within intraplate tectonic environments,” in Proceedings of the Third U.S. National Conference on Earthquake Engineering, Charleston, SC, 301–312.

Google Scholar

Cornell, C. A. (1968). Engineering seismic risk analysis. Bull. Seismol. Soc. Am. 58, 1583–1606.

Google Scholar

Cramer, C. H., Petersen, M. D., and Reichle, M. S. (1996). A Monte Carlo approach in estimating uncertainty for a seismic hazard assessment of Los Angeles, Ventura, and Orange Counties, California. Bull. Seismol. Soc. Am. 86, 1681–1691.

Google Scholar

Cunha, T., Watts, A., Pinheiro, L., and Myklebust, R. (2010). Seismic and gravity anomaly evidence of large-scale compressional deformation off SW Portugal. Earth Planet. Sci. Lett. 293, 171–179. doi: 10.1016/j.epsl.2010.01.047

CrossRef Full Text | Google Scholar

Davies, G., and Griffin, J. (2019). Sensitivity of probabilistic tsunami hazard assessment to far-field earthquake slip complexity and rigidity depth-dependence: case study of Australia. Pure Appl. Geophys. 177, 1521–1548. doi: 10.1007/s00024-019-02299-w

CrossRef Full Text | Google Scholar

Davies, G., Griffin, J., Løvholt, F., Glimsdal, S., Harbitz, C., Thio, H. K., et al. (2018). A global probabilistic tsunami hazard assessment from earthquake sources. Geol. Soc. Spec. Publ. 456, 219–244. doi: 10.1144/sp456.5

CrossRef Full Text | Google Scholar

Davis, R. A., Fox, W. T., Hayas, M. O., and Boothroyd, J. C. (1972). Comparison of ridge and runnel systems in tidal and non-tidal environments. J. Sediment. Petrol. 42, 413–421.

Google Scholar

Davison, F. C., and Scholz, C. H. (1985). Frequency-moment distribution of earthquakes in the Aleutian Arc: a test of the characteristic earthquake model. Bull. Seismol. Soc. Am. 75, 1349–1361.

Google Scholar

De Risi, R., and Goda, K. (2016). Probabilistic earthquake-tsunami multi-hazard analysis: application to the Tohoku region, Japan. Front. Built Environ. 2, 2–25. doi: 10.3389/fbuil.2016.00025

CrossRef Full Text | Google Scholar

Ebel, J. E., and Kafka, A. L. (1999). A Monte Carlo approach to seismic hazard analysis. Bull. Seismol. Soc. Am. 89, 854–866.

Google Scholar

El-Hussain, I., Omira, R., Deif, A., Al-Habsi, Z., Al-Rawas, G., Mohamad, A., et al. (2016). Probabilistic tsunami hazard assessment along Oman coast from submarine earthquakes in the Makran subduction zone. Arab. J. Geosci. 9:668. doi: 10.1007/s12517-016-2687-0

CrossRef Full Text | Google Scholar

FEMA (2000). Coastal Construction Manual, Tech. Rep. 55. Washington, DC: Federal Emergency Management Agency.

Google Scholar

Fukao, Y. (1973). Thrust faulting at a lithospheric plate boundary the portugal earthquake of 1969. Earth Planet. Sci. Lett. 18, 205–216. doi: 10.1016/0012-821x(73)90058-7

CrossRef Full Text | Google Scholar

Garcia-Orellana, J., Gràcia, E., Vizcaino, A., Masqué, P., Olid, C., Martínez- Ruiz, F., et al. (2006). Identifying instrumental and historical earthquake records in the SW Iberian margin using 210Pb turbidite chronology. Geophys. Res. Lett. 33:L24601. doi: 10.1029/2006GL028417

CrossRef Full Text | Google Scholar

Gardner, J. K., and Knopoff, L. (1974). Is the sequence of earthquakes in Southern California, with aftershocks removed, Poissonian? Bull. Seismol. Soc. Am. 64, 1363–1367.

Google Scholar

Geissler, W. H., Matias, L., Stich, D., Carrilho, F., Jokat, W., Monna, S., et al. (2010). Focal mechanisms for sub-crustal earthquakes in the Gulf of Cadiz from a dense OBS deployment. Geophys. Res. Lett. 37:L18309. doi: 10.1029/2010GL044289

CrossRef Full Text | Google Scholar

Geist, E., and Parsons, T. (2006). Probabilistic analysis of tsunami hazards. Nat. Hazards 37, 277–314. doi: 10.1007/s11069-005-4646-z

CrossRef Full Text | Google Scholar

Geist, E. L., and Lynett, P. J. (2014). Source processes for the probabilistic assessment of tsunami hazards. Oceanography 27, 86–93.

Google Scholar

George, D. L., and LeVeque, R. J. (2006). Finite volume methods and adaptive refinement for global tsunami propagation and local inundation. Sci. Tsunami Hazards 24, 319–328.

Google Scholar

Gica, E., Teng, M. H., Liu, P. L. F., Titov, V., and Zhou, H. (2007). Sensitivity analysis of source parameters for earthquake-generated distant tsunamis. J. Waterw. Port Coast. Ocean Eng. 133, 429–441. doi: 10.1061/(asce)0733-950x(2007)133:6(429)

CrossRef Full Text | Google Scholar

Gjevik, B., Pedersen, G., Dybesland, E., Harbitz, C. B., Miranda, P. M. A., Baptista, M. A., et al. (1997). Modeling tsunamis from earthquake sources near Gorringe Bank southwest of Portugal. J. Geophys. Res. Oceans 102, 27931–27949. doi: 10.1029/97jc02179

CrossRef Full Text | Google Scholar

Goda, K., Yasuda, T., Mori, N., and Mai, P. M. (2015). Variability of tsunami inundation footprints considering stochastic scenarios based on a single rupture model: application to the 2011 T ohoku earthquake. J. Geophys. Res. Oceans 120, 4552–4575. doi: 10.1002/2014jc010626

CrossRef Full Text | Google Scholar

González, F. I., Geist, E. L., Jaffe, B., Kânoglu, U., Mofjeld, H., Synolakis, C. E., et al. (2009). Probabilistic tsunami hazard assessment at seaside, Ore- gon, for near- and far-field seismic sources. J. Geophys. Res. 114:C11023. doi: 10.1029/2008JC005132

CrossRef Full Text | Google Scholar

González, M., Medina, R., Olabarrieta, M., and Otero, L. (2010). Tsunami hazard assesment in the Southern coast of Spain. Turk. J. Earth Sci. 19, 351–366.

Google Scholar

Grandin, R., Borges, J. F., Bezzeghoud, M., Caldeira, B., and Carrilho, F. (2007). Simulations of strong ground motion in SW Iberia for the 1969 February 28 (Ms= 8.0) and the 1755 November 1 (M ∼ 8.5) earthquakes – II. Strong ground motion simulations. Geophys. J. Int. 171, 807–822. doi: 10.1111/j.1365-246x.2007.03571.x

CrossRef Full Text | Google Scholar

Grezio, A., Babeyko, A., Baptista, M. A., Behrens, J., Costa, A., Davies, G., et al. (2017). Probabilistic tsunami hazard analysis: multiple sources and global applications. Rev. Geophys. 55, 1158–1198. doi: 10.1002/2017RG000579

CrossRef Full Text | Google Scholar

Grezio, A., Sandri, L., Marzocchi, W., Argnani, A., and Gasparini, P. (2012). Probabilistic tsunami hazard assessment for messina strait area (Sicily, Italy). Nat. Hazards 64, 329–358. doi: 10.1007/s11069-012-0246-x

CrossRef Full Text | Google Scholar

Griffin, J., Pranantyo, I. R., Kongko, W., Haunan, A., Robiana, R., Miller, V., et al. (2016). Assessing tsunami hazard using heterogeneous slip models in the Mentawai Islands, Indonesia. Geol. Soc. Spec. Publ. 441, 47–70. doi: 10.1144/SP441.3

CrossRef Full Text | Google Scholar

Gutscher, M.-A., Dominguez, S., Westbrook, G. K., and Leroy, P. (2009). Deep structure, recent deformation and analog modeling of the Gulf of Cadiz accretionary wedge: implications for the 1755 Lisbon earthquake. Tectonophysics 475, 85–97. doi: 10.1016/j.tecto.2008.11.031

CrossRef Full Text | Google Scholar

Heidarzadeh, M., and Kijko, A. (2011). A probabilistic tsunami hazard assessment for the Makran subduction zone at the northwestern Indian Ocean. Nat. Hazards 56, 577–593. doi: 10.1007/s11069-010-9574-x

CrossRef Full Text | Google Scholar

Hoechner, A., Babeyko, A. Y., and Zamora, N. (2016). Probabilistic tsunami hazard assessment for the Makran region with focus on maximum magnitude assumption. Nat. Hazards Earth Syst. Sci. 16, 1339–1350. doi: 10.5194/nhess-16-1339-2016

CrossRef Full Text | Google Scholar

Horspool, N., Pranantyo, I., Griffin, J., Latief, H., Natawidjaja, D. H., Kongko, W., et al. (2014). A probabilistic tsunami hazard assessment for Indonesia. Nat. Hazards Earth Syst. Sci. 14, 3105–3122. doi: 10.5194/Nhess-14-3105-2014

CrossRef Full Text | Google Scholar

IGN (2009). Catálogo de Tsunamis en las Costas Españolas. Madrid: Instituto Geográfico Nacional.

Google Scholar

IGN (2011). Catálogo Sísmico. Madrid: Instituto Geográfico Nacional.

Google Scholar

Kanamori, H. (1977). The energy release in great earthquakes. J. Geophys. Res. 82, 2981–2987. doi: 10.1029/jb082i020p02981

CrossRef Full Text | Google Scholar

Knighton, J., and Bastidas, L. A. (2015). A proposed probabilistic seismic tsunami hazard analysis methodology. Nat. Hazards 78, 699–723. doi: 10.1007/s11069-015-1741-7

CrossRef Full Text | Google Scholar

Kotani, T., Tozato, K., Takase, S., Moriguchi, S., Terada, K., Fukutani, Y., et al. (2020). Probabilistic tsunami hazard assessment with simulation-based response surfaces. Coast. Eng. 160:103719. doi: 10.1016/j.coastaleng.2020.103719

CrossRef Full Text | Google Scholar

Kowalik, Z., and Proshutinsky, A. (2010). Tsunami–tide interactions: a cook inlet case study. Cont. Shelf Res. 30, 633–642. doi: 10.1016/j.csr.2009.10.004

CrossRef Full Text | Google Scholar

Kowalik, Z., Proshutinsky, T., and Proshutinsky, A. (2006). Tide-tsunami interactions. Sci. Tsunami Hazards 24, 242–256.

Google Scholar

Laughton, A. S., and Whitmarsh, R. B. (1974). The Azores-Gibraltar Plate Boundary. Dordrecht: Reidel Publishing, 63–81.

Google Scholar

Lee, H. S., Shimoyama, T., and Popinet, S. (2015). Impacts of tides on tsunami propagation due to potential Nankai Trough earthquakes in the Seto Inland Sea, Japan. J. Geophys. Res. Oceans 120, 6865–6883. doi: 10.1002/2015jc010995

CrossRef Full Text | Google Scholar

Leonard, M. (2010). Earthquake fault scaling: self-consistent relating of rupture length, width, average displacement, and moment release. Bull. Seismol. Soc. Am. 100, 1971–1988. doi: 10.1785/0120090189

CrossRef Full Text | Google Scholar

Li, L., Switzer, A. D., Chan, C.-H., Wang, Y., Weiss, R., and Qiu, Q. (2016). How heterogeneous coseismic slip affects regional probabilistic tsunami hazard assessment: a case study in the South China Sea. J. Geophys. Res. 121, 6250–6272. doi: 10.1002/2016jb013111

CrossRef Full Text | Google Scholar

Lin, I.-C., and Tung, C. C. (1982). A preliminary investigation of tsunami hazard. Bull. Seismol. Soc. America Am. 72, 2323–2337.

Google Scholar

Lima, V. V., Miranda, J. M., Baptista, M. A., Catalão, J., Gonzalez, M., Otero, L., et al. (2010). Impact of a 1755-like tsunami in Huelva, Spain. Nat. Hazards Earth Syst. Sci. 10, 139–148. doi: 10.5194/nhess-10-139-2010

CrossRef Full Text | Google Scholar

Liu, P. L.-F., Cho, Y.-S., Briggs, M. J., Kanoglu, U., and Synolakis, C. E. (1995). Runup of solitary waves on a circular island. J. Fluid Mech. 302, 259–285. doi: 10.1017/s0022112095004095

CrossRef Full Text | Google Scholar

Lorito, S., Selva, J., Basili, R., Romano, F., Tiberti, M. M., and Piatanesi, A. (2015). Probabilistic hazard for seismically induced tsunamis: accuracy and feasibility of inundation maps. Geophys. J. Int. 200, 574–588. doi: 10.1093/gji/ggu408

CrossRef Full Text | Google Scholar

Lorito, S., Tiberti, M. M., Basili, R., Piatanesi, A., and Valensise, G. (2008). Earthquake generated tsunamis in the Mediterranean Sea: scenarios of potential threats to southern Italy. J. Geophys. Res. Solid Earth 113, 1978–2012.

Google Scholar

Luque, L., Lario, J., Civis, J., Silva, P. G., Zazo, C., Goy, J. L., et al. (2002). Sedimentary record of a tsunami during Roman times, Bay of Cadiz, Spain. J. Quat. Sci. 17, 623–631. doi: 10.1002/jqs.711

CrossRef Full Text | Google Scholar

Luque, L., Lario, J., Zazo, C., Goy, J. L., Dabrio, C. J., and Silva, P. G. (2001). Tsunami deposits as paleoseismic indicators: examples from the Spanish coast. Acta Geol. Hisp. 36, 197–211.

Google Scholar

Main, I. (2000). Apparent breaks in scaling in the earthquake cumulative frequency- magnitude distribution: fact or artifact? Bull. Seismol. Soc. Am., 90, 86–97. doi: 10.1785/0119990086

CrossRef Full Text | Google Scholar

Martínez-Loriente, S., Gràcia, E., Bartolome, R., Sallarès, V., Connors, C., Perea, H., et al. (2013). «Active deformation in old oceanic lithosphere and significance for earthquake hazard: seismic imaging of the coral patch ridge area and neighboring abyssal plains (SW Iberian Margin)». Geochem. Geophys. Geosy. 14, 2206–2231. doi: 10.1002/ggge.20173

CrossRef Full Text | Google Scholar

Medialdea, T., Vegas, R., Somoza, L., Vázquez, J., Maldonado, A., del Río, V. D., et al. (2004). Structure and evolution of the “Olistostrome” complex of the Gibraltar Arc in the Gulf of Cádiz (east- ern Central Atlantic): evidence from two long seismic cross-sections. Mar. Geol. 209, 173–198. doi: 10.1016/j.margeo.2004.05.029

CrossRef Full Text | Google Scholar

Mezcua, J., and Martínez-Solares, J. M. (1983). Sismicidad en el Área Ibero-Mogrebí, Publicación Técnica 203. Madrid: Instituto Geográfico Nacional.

Google Scholar

Morales, J. A., Borrego, J., Miguel, E. G. San, López-González, N., and Carro, B. (2008). Sedimentary record of recent tsunamis in the Huelva Estuary (southwestern Spain). Quat. Sci. Rev. 27, 734–746. doi: 10.1016/j.quascirev.2007.12.002

CrossRef Full Text | Google Scholar

Murotani, S., Satake, K., and Fujii, Y. (2013). Scaling relations of seismic moment, rupture area, average slip, and asperity size for M˜ 9 subduction-zone earthquakes. Geophys. Res. Lett. 40, 5070–5074. doi: 10.1002/grl.50976

CrossRef Full Text | Google Scholar

National Research Council (1988). Probabilistic Seismic Hazard Analysis, Tech. Rep. Washington, DC: National Academy press.

Google Scholar

National Research Council (1997). Review of Recommendations for Probabilistic Seismic Hazard Analysis: Guidance on Uncertainty and Use of Experts, Tech. rep. Washington, DC: National Academy press.

Google Scholar

Okal, E. A., and Synolakis, C. E. (2004). Source discriminants for near-field tsunamis. Geophys. J. Int. 158, 899–912. doi: 10.1111/j.1365-246x.2004.02347.x

CrossRef Full Text | Google Scholar

Olabarrieta, M., Medina, R., Gonzalez, M., and Otero, L. (2011). C3: a finite volume- finite difference hybrid model for tsunami propagation and runup, computers andamp. Geosciences 37, 1003–1014. doi: 10.1016/j.cageo.2010.09.016

CrossRef Full Text | Google Scholar

Omira, R., Baptista, M. A., and Matias, L. (2015). Probabilistic tsunami hazard in the northeast atlantic from near-and far-field tectonic sources. Pure Appl. Geophys. 172, 901–920. doi: 10.1007/s00024-014-0949-x

CrossRef Full Text | Google Scholar

Omira, R., Baptista, M. A., Matias, L., Miranda, J. M., Catita, C., Carrilho, F., et al. (2009). Design of a sea-level tsunami detection network for the Gulf of Cadiz. Nat. Hazards Earth Syst. Sci. 9, 1327–1338. doi: 10.5194/nhess-9-1327-2009

CrossRef Full Text | Google Scholar

Omira, R., Baptista, M. A., and Miranda, J. M. (2011). Evaluating tsunami impact on the Gulf of Cadiz Coast (Northeast Atlantic). Pure Appl. Geophys. 168, 1033–1043. doi: 10.1007/s00024-010-0217-7

CrossRef Full Text | Google Scholar

Omira, R., Matias, L., and Baptista, M. A. (2016). Developing an event-tree probabilistic tsunami inundation model for NE Atlantic coasts: application to a case study. Pure Appl. Geophy. 173, 3775–3794. doi: 10.1007/s00024-016-1367-z

CrossRef Full Text | Google Scholar

Otero, L. (2008). Metodología Para Evaluar la Peligrosidad Debido a Tsunamis en Zonas Costeras. Ph.D. thesis. Santander: University of Cantabria.

Google Scholar

Pereira, R., Alves, T. M., and Cartwright, J. (2011). Post-rift compression on the SW Iberian margin (eastern North Atlantic): a case for prolonged inversion in the ocean–continent transition zone. J. Geol. Soc. 168, 1249–1263. doi: 10.1144/0016-76492010-151

CrossRef Full Text | Google Scholar

Power, W., Wallace, L., Wang, X., and Reyners, M. (2012). Tsunami hazard posed to New Zealand by the Kermadec and southern New Hebrides subduction margins: an assessment based on plate boundary kinematics, interseismic coupling, and historical seismicity. Pure Appl. Geophys. 169, 1–36.

Google Scholar

Power, W., Downes, G., and Stirling, M. (2007). Estimation of tsunami hazard in New Zealand due to South American earthquakes. Pure Appl. Geophys. 164, 547–564.

Google Scholar

Power, W., Wang, X., Lane, E. M., and Gillibrand, P. A. (2013). A probabilistic tsunami hazard study of the Auckland region, Part I: propagation modelling and tsunami hazard assessment at the shoreline. Pure Appl. Geophys. 170, 16–21. doi: 10.1007/s00024-012-0543-z

CrossRef Full Text | Google Scholar

Rikitake, B. T., and Aida, I. (1988). Tsunami hazard probability in Japan. Bull. Seismol. Soc. Am. 78, 1268–1278.

Google Scholar

Roger, J., Baptista, M. A., Sahal, A., Accary, F., Allgeyer, S., and Hébert, H. (2011). The Transoceanic 1755 Lisbon tsunami in Martinique. Pure Appl. Geophys. 168, 1015–1031.

Google Scholar

Rosas, F., Duarte, J., Terrinha, P., Valadares, V., and Matias, L. (2009). Morphotectonic characterization of major bathymetric lineaments in Gulf of Cadiz (Africa– Iberia plate boundary): insights from analogue modelling experiments. Mar. Geol. 261, 33–47.

Google Scholar

Rueda, J., and Mezcua, J. (2002). Estudio del terremoto de 23 de septiembre de 2001 en Pego (Alicante). Obtención de una relación m_bLg - M_W para la Península Ibérica. Rev. Soc. Geol. Esp. 15, 159–173.

Google Scholar

Ruiz, F., Rodríguez-Ramírez, A., Cáceres, L. M., Vidal, J. R., Carretero, M. I., Abad, M., et al. (2005). Evidence of high-energy events in the geological record: mid-holocene evolution of the southwestern Doñana National Park (SW Spain). Palaeogeogr. Palaeoclimatol. Palaeoecol. 229, 212–229.

Google Scholar

Sartori, R., Torelli, L., Zitellini, N., Peis, D., and Lodolo, E. (1994). Eastern segment of the Azores-Gibraltar line (central-eastern Atlantic): an oceanic plate boundary with diffuse compressional deformation. Geology 22, 555–558.

Google Scholar

Savage, J. C. (1991). Criticism of some forecasts of the National Earthquake Prediction Evaluation Council. Bull. Seismol. Soc. Am. 81, 862–881.

Google Scholar

Savage, J. C. (1992). The uncertainty in earthquake conditional probabilities. Geophys. Res. Lett. 19, 709–712.

Google Scholar

Scholz, C. H. (1998). «Earthquakes and friction laws». Nature 391, 37–42.

Google Scholar

Scholz, C. H. (2002). The Mechanics of Earthquakes and Faulting, 2nd Edn. Cambridge: Cambridge University Press, 471.

Google Scholar

Schwartz, D. P., and Coppersmith, K. J. (1984). Fault behavior and characteristic earthquakes: examples from the wasatch and san andreas fault zones. J. Geophys. Res. 89, 5681–5698. doi: 10.1029/JB089iB07p05681

CrossRef Full Text | Google Scholar

Senior Seismic Hazard Analysis Comittee (1997). Recommendations for Probabilistic Seismic Hazard Analysis: Guidance on Uncertainty and Use of Experts, Tech. rep. Rockville, MD: U. S. Nuclear Regulatory Commission.

Google Scholar

Sepúlveda, I., Liu, P. L. F., and Grigoriu, M. (2019). Probabilistic tsunami hazard assessment in South China Sea with consideration of uncertain earthquake characteristics. J. Geophys. Res. Solid Earth 124, 658–688. doi: 10.1029/2018JB016620

CrossRef Full Text | Google Scholar

Shelby, M., Grilli, S. T., and Grilli, A. R. (2016). Tsunami hazard assessment in the Hudson River Estuary based on dynamic tsunami-tide simulations. Pure Appl. Geophys. 173, 3999–4037.

Google Scholar

Shin, J. Y., Chen, S., and Kim, T.-W. (2015). Application of Bayesian Markov Chain Monte Carlo method with mixed gumbel distribution to estimate extreme magnitude of tsunamigenic earthquake. J. Civ. Eng. 19, 366–375. doi: 10.1007/s12205-015-0430-0

CrossRef Full Text | Google Scholar

Sibson, R. H. (1982). «Fault zone models, heat flow, and the depth distribution of earthquakes in the continental crust of the United States». Bull. Seismol. Soc. Am. 72, 151–163.

Google Scholar

Silva, S., Terrinha, P., Matias, L., Duarte, J. C., Roque, C., Ranero, C. R., et al. (2017). Micro-seismicity in the Gulf of Cadiz: is there a link between micro-seismicity, high magnitude earthquakes and active faults? Tectonophysics 717, 226–241. doi: 10.1016/j.tecto.2017.07.026

CrossRef Full Text | Google Scholar

Sørensen, M. B., Spada, M., Babeyko, A., and Wiemer, S., and Grunthal, G. (2012). Probabilistic tsunami hazard in the Mediterranean Sea. J. Geophys. Res. 117:B01305. doi: 10.1029/2010JB008169

CrossRef Full Text | Google Scholar

Sorensen, M. B., Spada, M., Babeyko, A., Wiemer, S., and Grunthal, G. (2012). Probabilistic tsunami hazard in the Mediterranean Sea. J. Geophys. Res. 117:B01305. doi: 10.1029/2010JB008169

CrossRef Full Text | Google Scholar

Stich, D., de Lis Mancilla, F., and Morales, J. (2005). Crust-mantle coupling in the Gulf of Cadiz (SW-Iberia). Geophys. Res. Lett. 32:L13306.

Google Scholar

Stich, D., Martín, R., and Morales, J. (2010). Moment tensor inversion for Iberia– Maghreb earthquakes 2005–2008. Tectonophysics 483, 390–398.

Google Scholar

Stich, D., Serpelloni, E., de Lis Mancilla, F., and Morales, J. (2006). Kinematics of the Iberia–Maghreb plate contact from seismic moment tensors and GPS observations. Tectonophysics 426, 295–317.

Google Scholar

Stirling, M., Goded, T., Berryman, K., and Litchfield, N. (2013). Selection of earthquake scaling relationships for seismic-hazard analysis. Bull. Seismol. Soc. Am. 103, 2993–3011.

Google Scholar

Strasser, F. O., Arango, M. C., and Bommer, J. J. (2010). Scaling of the source dimensions of interface and intraslab subduction-zone earthquakes with moment magnitude. Seismol. Res. Lett. 81, 941–950. doi: 10.1785/gssrl.81.6.941

CrossRef Full Text | Google Scholar

Suppasri, A., Imamura, F., and Koshimura, S. (2012). Probabilistic tsunami hazard analysis and risk to coastal populations in Thailand. J. Earthq. Tsunami 06:1250011. doi: 10.1142/S179343111250011x

CrossRef Full Text | Google Scholar

Synolakis, C. E. (2003). Tsunami and Seiche, pp. 9.1–9.99, Earthquake Engineering Handbook. New York, NY: CRC Press.

Google Scholar

Synolakis, C. E., Bernard, E. N., Titov, V. V., Kanoglu, U., and GonzáLez, F. I. (2007). Standard, Criteria, and Procedures for NOAA Evaluation Of Tsunami Numerical Models, Technical Memorandum PMEL-135. Washington, DC: NOAA.

Google Scholar

Terrinha, P., Matias, L., Vicente, J., Duarte, J., Luís, J., Pinheiro, L., et al. (2009). Morphotectonics and strain partitioning at the Iberia–Africa plate boundary from multibeam and seismic reflection data. Mar. Geol. 267, 156–174.

Google Scholar

Taubenböck, H., Post, J., Roth, A., Zosseder, K., Strunz, G., and Dech, S. (2008). A conceptual vulnerability and risk framework as outline to identify capabilities of remote sensing. Nat. Hazard. Earth Syst. Sci. 8, 409–420. doi: 10.5194/nhess-8-409-2008

CrossRef Full Text | Google Scholar

Thio, H. K., Somerville, P., Polet, J., and Ichinose, G. (2007). “Probabilistic tsunami hazard analysis for ports and harbors. Ports 2007: 30 Years of Sharing Ideas 1977-2007,” in Proceedings of the Eleventh Triennial International Conference. doi: 10.1061/40834(238)103

CrossRef Full Text | Google Scholar

Tinti, S., and Armigliato, A. (2003). The use of scenarios to evaluate the tsunami impact in southern Italy. Mar. Geol. 199, 221–243.

Google Scholar

Tinti, S., Armigliato, A., Pagnoni, G., and Zaniboni, F. (2005). Scenarios of giant tsunamis of tectonic origin in the Mediterranean. ISET J. Earthq. Technol. 42, 171–188.

Google Scholar

Tolkova, E. (2016). Tsunami penetration in tidal rivers, with observations of the Chile 2015 tsunami in rivers in Japan. Pure Appl. Geophys. 173, 389–409.

Google Scholar

Volpe, M., Lorito, S., Selva, J., Tonini, R., Romano, F., and Brizuela, B. (2019). From regional to local SPTHA: efficient computation of probabilistic tsunami inundation maps addressing near-field sources. Nat. Hazard. Earth Syst. Sci. 19, 455–469.

Google Scholar

Ward, S. N. (1991). A synthetic seismicity model for the middle America trench. J. Geophys. Res. 96, 433–421.

Google Scholar

Ward, S. N. (1996). A synthetic seismicity model for southern California: cycles, probabilities, and hazard. J. Geophys. Res. 101, 393–322.

Google Scholar

Ward, S. N. (2000). San Francisco Bay Area earthquake simulations: a step toward a standard physical earthquake model. Bull. Seismol. Soc. Am. 90, 370–386.

Google Scholar

Weisz, R., and Winter, C. (2005). “Tsunami, tides and run-up: a numerical study,” in Proceedings of the International Tsunami Symposium, eds G. A. Papadopoulos and K. Satake, Chania, 322.

Google Scholar

Wells, D. L., and Coppersmith, K. J. (1994). New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. Bull. Seismol. Soc. Am. 84, 974–1002.

Google Scholar

Wesnousky, S. G. (1994). The Gutenberg-Richter or characteristic earthquake distribution, which is it? Bull. Seismol. Soc. Am. 84, 1940–1959.

Google Scholar

Wessel, P., and Smith, W. H. F. (1998). New improved version of generic mapping tools released. EOS Trans. Am. Geophys. Union 79:579.

Google Scholar

Weatherall, P., Marks, K. M., Jakobsson, M., Schmitt, T., Tani, S., Arndt, J. E., et al. (2015). A new digital bathymetric model of the world’s oceans. Earth Space Sci. 2, 331–345. doi: 10.1002/2015EA000107

CrossRef Full Text | Google Scholar

Wiemer, S. (2001). A software package to analyze seismicity: ZMAP. Seismol. Res. Lett. 72, 373–382.

Google Scholar

Yadav, R. B. S., Tsapanos, T. M., Bayrak, Y., and Koravos, G. C. (2013). Probabilistic appraisal of earthquake hazard parameters deduced from a bayesian approach in the northwest frontier of the himalayas. Pure Appl. Geophys. 170, 283–297. doi: 10.1007/s00024-012-0488-2

CrossRef Full Text | Google Scholar

Yamada, Y., Masui, R., and Tsuji, T. (2013). Characteristics of a tsunamigenic megasplay fault in the Nankai Trough. Geophys. Res. Lett. 40, 4594–4598.

Google Scholar

Yeh, H. (2006). “Tsunami forces in the runup zone,” in NSF Caribbean Tsunami Workshop, eds A. Mercado-Irizarry and P. L.-F. Liu (Singapore: World Scientific).

Google Scholar

Zamora, N., and Babeyko, A. Y. (2020). Probabilistic tsunami hazard assessment for local and regional seismic sources along the Pacific coast of Central America with emphasis on the role of selected uncertainties. Pure Appl. Geophys. 177, 1471–1495.

Google Scholar

Zhang, Y. J., Witter, R. C., and Priest, G. R. (2011). Tsunami–tide interaction in 1964 Prince William Sound tsunami. Ocean Model. 40, 246–259.

Google Scholar

Zielke, O., and Arrowsmith, J. R. (2008). Depth variation of coseismic stress drop explains bimodal earthquake magnitude-frequency distribution. Geophys. Res. Lett. 35, L24301. doi: 10.1029/2008GL036249

CrossRef Full Text | Google Scholar

Zitellini, N., Gràcia, E., Matias, L., Terrinha, P., Abreu, M., DeAlteriis, G., et al. (2009). The quest for the Africa–Eurasia plate boundary west of the Strait of Gibraltar. Earth Planet. Sci. Lett. 280, 13–50.

Google Scholar