Modeling tsunami propagation in the Iberia–Africa plate boundary: Historical events, regional exposure and the case-study of the former Gulf of Tartessos

Modeling tsunami propagation in the Iberia–Africa plate boundary: Historical events, regional exposure and the case-study of the former Gulf of Tartessos R. Periáñez ⁎, J.M. Abril Dpto. Física Aplicada I, ETSIA, Universidad de Sevilla, Ctra. Utrera km 1, 41013-Sevilla, Spain abstractarticle info Article history: Received 17 September 2012 Received in revised form 12 November 2012 Accepted 14 November 2012 Available online 27 November 2012 Keywords: Tsunami modeling Iberia Risk assessment Gulf of Tartessos A numerical model to simulate tsunami propagation in south of Iberia waters has been developed. It is based on the 2D non-linear hydrodynamic equations and allows calculating tsunami run-ups. The model has been validated through the simulation of historical tsunamis. Then it is applied to a risk assessment study to evaluate tsunami flooding along the Spanish and Moroccan coasts. A spectral analysis of tsunami waves has also been carried out. Finally, the model has been applied to simulate tsunami propagation in the former Gulf of Tartessos, where the city of the same name could have been located. The objective consists of evaluating if a large tsunami could have destroyed it, as some historians suggest. © 2012 Elsevier B.V. All rights reserved. 1. Introduction The interest in modeling tsunamis has increasedin recentyears after the terrible consequences of the 2004 Indian Ocean tsunami and of the 11th March 2011 tsunami in Japan. Tsunamis can produce primary catastrophic effects through strong water currents and waves of large amplitude which break on the coastline, being able to flood large flat areas. But catastrophic destruction may also occur when the frequencies of the arriving tsunami waves match the resonant frequencies of a harbor or bay (Rabinovich, 2009). Thus, the Spanish International Cooperation Agency has recently funded a project whose objective consists of the development of a tsunami propagation numerical model for southern Iberia and Moroccan waters (Fig. 1). In fact, the north-African margin and south of Iberia are areas seismically active since African and Eurasian plates meet here. As a consequence several notable tsunamis have been triggered in historic times (Rodríguez-Vidal et al., 2011)aswellasrecently:we may consider the 1755 Lisbon tsunami (Baptista et al., 2003), the 1856 and 2003 tsunamis in Algeria (Roger and Hébert, 2008; Sahal et al., 2009) and the 1969 tsunami in the Atlantic (Guesmia et al., 1998). Some models for predicting tsunami propagation and arrival times have already been developed in the past (see for instance Choi et al., 2003; Ioualalen et al., 2010; Kowalik et al., 2007; Ma and Lee, 1997). In recent years, an intensive search has been conducted in the Gulf of Cádiz to identify active faults being able to generate large earthquakes and tsunamis (Duarte et al., 2011; Zitellini et al., 2009). Although it cannot be considered a closed question, five candidate sources for tsunamis have been identified (Lima et al., 2010). These kinds of studies have also been carried out for the US East Coast (Barkan et al., 2009). Also, spectral analysis has been applied, among other uses, to compare recorded time series of water elevations induced by a tsunami against synthetic ones, from numerical simulation (Miranda et al., 1996), and for studying the evolution of the periodic components during the tsunami propagation from its origin (Lobkovsky et al., 2010; Sahal et al., 2009). The objective of this paper consists of describing the tsunami propagation numerical model which has been developed in the frame of the research project mentioned above. The model domain (Fig. 1) covers the Gulf of Cádiz (Atlantic Ocean) and the western Mediterranean, which are connected through the Strait of Gibraltar. Differently to some other previous models, it is fully non-linear and also allows run-up calculations. The model has been tested through the simulation of the 1755 Lisbon tsunami and the 1856 and 2003 Algeria tsunamis, for which historical data or previous model simulations exist, to compare with the results of our model. Once tested, the model is applied to three case studies: the first risk assessment carried out along the Spanish–Moroccan coasts, spectral analysis of tsunami waves (on account of its interest) and, finally, the model has been applied to simulate tsunami propagation in the former Gulf of Tartessos, in SW Spain (see Fig. 1 for its location). It was a large marine gulf which extended over the nowadays Guadalquivir Marsh since the last Holocene sea level rise (some 6900 years BP) until it was finally colmated. It is speculated that the city of Tartessos, located in this area, could have been destroyed by a tsunami occurring during the sixth century BC. The bathymetry of the Gulf of Tartessos has been reconstructed and the effects of a large tsunami have been Journal of Marine Systems 111–112 (2013) 223–234 ⁎Corresponding author. Fax: +34 954486436. E-mail address: [email protected] (R. Periáñez). 0924-7963/$ –see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jmarsys.2012.11.005 Contents lists available at SciVerse ScienceDirect Journal of Marine Systems journal homepage: www.elsevier.com/locate/jmarsys simulated by the model so as to evaluate if such a hypothesis is reasonable. The model is described in Section 2. Next, the results are presented and discussed. Initially (Section 3), historical tsunamis are simulated and results compared with previous simulations and/or historical data, to evaluate the performance of our model. Then, model applications (risk assessment, spectral analysis and Gulf of Tartessos studies) are presented in Section 4. 2. Model description The 2D depth-averaged barotropic hydrodynamic equations, which describe the propagation of surface shallow water gravity waves, are (see for instance Kowalik and Murty, 1993): ∂ζ ∂tþ∂ ∂xHuðÞþ ∂ ∂yHvðÞ¼0ð1Þ ∂u ∂tþu∂u ∂xþv∂u ∂yþg∂ζ ∂x−Ωvþτu ρH¼A∂2u ∂x2þ∂2u ∂y2 ! ð2Þ ∂v ∂tþu∂v ∂xþv∂v ∂yþg∂ζ ∂yþΩuþτv ρH¼A∂2v ∂x2þ∂2v ∂y2 ! ð3Þ where uand vare the depth averaged water velocities along the xand y axes, his the depth of water below the mean sea level, ζis the displacement of thewater surface above the mean sea level measured upwards, H=h+ζis the total water depth, Ωis the Coriolis parameter (Ω= 2wsinλ, where wis the Earth rotational angular velocity and λis latitude), gis acceleration due to gravity, ρis a mean value of water density and Ais the horizontal eddy viscosity. τ u and τ v are friction stresses which have been written in terms of a quadratic law: τu¼kρuffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2þv2 p τv¼kρvffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2þv2 pð4Þ where kis the bed friction coefficient. Essentially, these equations express mass and momentum conservation. They have been written in Cartesian coordinates given the relatively small model domain. At oceanic scale spherical coordinates would be required (Kowalik and Murty, 1993). All the equations are solved using explicit finite difference schemes (Kowalik and Murty, 1993) with second order accuracy. In particular, the MSOU (Monotonic Second Order Upstream) is used for the advective non-linear terms in the momentum equations. The computational domain is presented in Fig. 1. Water depths have been obtained from the GEBCO08 digital atlas, available on-line, with a resolution of 30 seconds of arc both in longitude and latitude. Thus, the domain is covered by a total number of 2,052,000 grid cells. Due to the CFL stability condition (Kowalik and Murty, 1993) time step for model integration was fixed as Δt=1 s. Boundary conditions have to be specified. There is no water flow towards a dry grid cell. Wet and dry cells are identified through the flood/dry algorithm described below. A gravity wave radiation condition (Herzfeld et al., 2011)isusedfor sea surface elevation along open boundaries. In this approximation it is assumed that processes near the boundaries are always governed by outgoing waves. This is a valid assumption since tsunamis will propagate outside from the interior of the domain. The boundary condition is implemented in an implicit form, thus: ζxB;tþ1ðÞ¼ ζxB;tðÞþμζ xB−1;tþ1ðÞ 1þμð5Þ where t+1 is the updated time step, x B designs the boundary location and x B−1 means the position of the nearest interior point to the open boundary. Also, μ¼cΔt Δxð6Þ where c¼ffiffiffiffiffiffiffi gH pis the shallow water barotropic phase speed. This boundary condition allows the transit of information across the open boundary without spurious reflection and/or amplification. Note that tsunamis may be treated as shallow water waves since wavelengths are significantly larger than water depths. Aflood/dry algorithm is required since when the tsunami reaches the coast new wet or dry grid cells may be generated due to run-up or rundown. The numerical scheme described in Kampf (2009) has been Fig. 1. Model domain. Water depths (GEBCO08 bathymetry, 30 seconds of arc resolution) are given in m. Places mentioned in the text are located. The gray box indicates the area of the former Gulf of Tartessos. 224 R. Periáñez, J.M. Abril / Journal of Marine Systems 111–112 (2013) 223–234 adopted. Wet grid cells are defined as those with a total water depth Hlarger than a threshold value H min typically set as a few centimeters. Dry cells are defined as cells where H≤H min . Flooding and drying are implemented in the code via the calculation of the water velocity normal to the interface between wet and dry cells. The calculation is performed when the pressure gradient force is directed towards the dry cell. Otherwise velocity is set to zero at this point. In the case of a non-zero velocity, water level in the dry cell will increase and the cell turns into a wet one once the water depth is larger than H min .In the present study, H min =0.10 m has been fixed. Tsunami initial conditions also have to be specified. Owing to the incompressibility of water, sea surface is considered to be displaced from the mean level the same magnitude as the sea bottom is displaced from its original position due to the earthquake (Choi et al., 2003; Ma and Lee, 1997;Sahal et al., 2009). We also assume that this deformation is instantaneous: rupture duration (a few tens of seconds) is much smaller than tsunami periods (Alasset et al., 2006; Harbitz et al., 2012; Ma and Lee, 1997). Therefore, the vertical sea-floor deformation is considered as the initial condition for the tsunami calculation, together with zero flow over all the domain (Cecioni and Bellotti, 2010; Choi et al., 2003; Ioualalen et al., 2010;Ma and Lee, 1997; Sahal et al., 2009). The sea-floor deformation produced by the earthquake is computed using the classical Okada formulae (Okada, 1985). Inputs to these formulae are fault plane location, rake, strike, dip, slip, length and width, as well as seismic moment and rigidity. This also allows the estimation of the total potential energy related to the bottom deformation which is transferred to the sea surface (Kowalik et al., 2007). AFORTRAN code has been written to implement model equations. The hydrodynamic model and numericalschemes have been tested through simulation of barotropic tides in these waters (Periáñez, 2007, 2009), which allowed selecting optimum values for the bed friction coefficient and horizontal viscosity. Their values are, respectively, k=0.0025 and A=10 m 2 /s. The code runs on a PC with two Intel X6800 2.93 GHz processors and Ubuntu operating system. The Intel FORTRAN compiler has been used. Computation times are long due to the large number of grid cells and small time step. One hour of simulation time takes 80 min. 3. Model testing Three historical tsunamis have been simulated. The objective of these simulations is to carry out an intercomparison exercise with previous models. Once we know that model results are comparable with previous modeling studies, some applications are described in Section 4. Table 1 Fault parameters used in the simulations. Geographical coordinates correspond to the fault center. Rake is 90° in all cases. [1] Barkan et al. (2009), [2] Alasset et al. (2006), [3] Roger and Hébert (2008) and [4] Lima et al. (2010) and Zitellini et al. (2009). The five fault resulting seafloor deformations are drawn in Fig. 6. Tsunami Lon Lat Length (km) Width (km) Slip (m) Strike (deg) Dip (deg) Reference Lisbon (source 5 in reference [1]) −10.753 36.042 200 80 19 345 40 [1] Lisbon (source 8 in reference [1]) −11.467 36.015 200 80 19 345 40 [1] Algeria 2003 3.635 36.893 60 24 0.71 70 45 [2] Algeria 1856 5.736 37.080 110 20 1.5 60 40 [3] GBF (Gorringe Bank Fault) −11.281 36.948 137 60 8.3 233 25 [4] HSF (Horseshoe Fault) −10.000 36.133 106 70 10.7 222.1 25 [4] MPF (Marqués de Pombal Fault) −10.067 36.895 86 70 8.0 200.0 25 [4] PBF (Portimao Bank Fault) −8.664 36.105 100 55 7.2 266.3 25 [4] CWF (Cádiz Wedge Fault) −9.329 34.833 133 200 11.1 346.3 6 [4] Fig. 2. Computed maximum wave amplitudes (m) after the 1755 Lisbon tsunami. Simulation time 3 h. 225R. Periáñez, J.M. Abril / Journal of Marine Systems 111–112 (2013) 223–234 3.1. The 1755 Lisbon tsunami The tsunami generated by the Lisbon earthquake of November 1st, 1755, is well documented in historical accounts. It reached not only Spain, Morocco and Portugal, but also Ireland, England and the Caribbean. In spite of the extensive work carried out to determine the earthquake epicenter (Baptista et al., 1998, 2003; Barkan et al., 2009), it still remains an open problem. The most recent study by Barkan et al. (2009) suggests that the most likely epicenter was located in the Horseshoe Plain area, and not in the previously proposed locations: active faults denoted as GBF (Gorringe Bank Fault), MPF (Marqués de Pombal Fault) and GCF (Gulf of Cádiz Fault). We have used source number 5 in Barkan et al. (2009) to carry out a simulation with our model. These authors found that two sources, denoted as 5 and 8 in their paper, produced the best agreement with historical data. After some test runs, slightly better results are produced by our model with source 5. Their characteristics are summarized in Table 1. The computed maximum amplitude of the tsunami with this source is presented in Fig. 2 after a simulation time of 3 h. Most energy is directed towards the open ocean, as in the simulations by Barkan et al. (2009). However, significant wave amplitudes are apparent in the area of Cape San Vicente, which focuses waves due to bathymetry, and along the southwestern Spanish and Portuguese coasts. Waves along Morocco are smaller than along SW Spain. Finally, the tsunami hardly penetrates the Mediterranean due to the presence of the narrow Strait of Gibraltar. Indeed, taking a characteristic depth of some 500 m in the Strait, the wave speed ffiffiffiffiffiffi gh pis70m/s.Thewaveperiod in Gibraltar (see Section 2) is about 1800 s. Thus, the corresponding wave length is about 126 km, much larger than the Strait of Gibraltar width. Consequently wave diffraction is likely to occur, with the corresponding amplitude decrease with increasing distance to the Strait. A comparison of model results with historical data is presented in Table 2. Generally speaking, model results are in acceptable agreement with historical data. Amplitudes are well predicted in Spanish and Portuguese locations, although underestimated in Ceuta and Gibraltar. Arrival times are consistent with historical reports, although tend to be overestimated in some locations. Nevertheless, the overall quality of results is similar to that of previous modeling studies (Baptista et al., 2003; Barkan et al., 2009). 3.2. The 1856 Algeria tsunami On August 1856 an earthquake with an estimated magnitude of 5.7 occurred about 300 km east of Algiers. There are no historical data on the estimated arrival times or on tsunami run-ups, but only some reports on the damage that occurred in the Balearic Islands and along the Algerian coasts. This tsunami and its possible source have been studied in detail by Roger and Hébert (2008). We have used the same source proposed by these authors (Table 1) and compared our modeling results with those presented by them. The computed maximum wave amplitude, aftera 2 hour simulation, may be seen in Fig. 3. Results are in very good agreement with the previous calculations of Roger and Hébert (2008) (Fig. 3B in their paper). Most energy is directed towards the southeast of Majorca Island and, to smaller degree, towards Menorca Island. Computed maximum amplitudes in these areas, about 0.80 m, are also in close agreement with those presented by Roger and Hébert (2008). 3.3. The 2003 Algeria tsunami More information is available for the tsunami triggered on May 21, 2003, by a 6.9 magnitude earthquake whose epicenter was 50 km east of the city of Algiers. This tsunami has been studied in detail Table 2 Run-up reports (m) and arrival times after the initial condition (min) from historical data and modeling results for the 1755 Lisbon tsunami. Arrival 1 refersto the computed arrival time of a 10 cm amplitude signal and arrival 2 to the computed arrival time of the first wave or subsidence maximum. Points are located in Fig. 1. Location Run-up (hist.) Arrival (hist.) Run-up (model) Arrival 1 Arrival 2 S. Vicente >10 16±7 17.4 8 22 Lisbon 5 25±10 4.3 28 36 Huelva 50±10 5.4 64 69 Cádiz 15 78±15 8.9 57 65 Gibraltar 2 52–68 0.53 76 79 Ceuta 2 0.38 76 79 Figueira 45±10 1.9 57 69 Sale 33–50 4.7 37 44 Fig. 3. Computed maximum wave amplitudes (m) after the 1856 Algeria tsunami. Simulation time 2 h. 226 R. Periáñez, J.M. Abril / Journal of Marine Systems 111–112 (2013) 223–234 before (Alasset et al., 2006; Sahal et al., 2009). These authors found that with the source parameters described by Delouis (details on it may be seen in the paper by Alasset and coworkers, see also Table 1), their computed travel times were in agreement with observations. The computed maximum wave amplitude, after a 2 hour simulation, may be seen in Fig. 4. This map is in excellent agreement with the computations of Alasset et al. (2006), using the Delouis source (Fig. 6a in their paper), and those of Sahal et al. (2009) using the same source (Fig. 5 in their paper). Energy is again mainly directed towards the Balearic Islands but, differently to the 1856 tsunami, Ibiza is more affected than Menorca and Majorca. Nevertheless, as in Alasset et al. (2006) and Sahal et al. (2009), maximum amplitudes tend to be underestimated by the models with respect to tidal gauge measurements. This has been attributed by Sahal et al. (2009) to the fact that resonance plays an essential role in wave amplification in bays, specially for small magnitude events. Such resonance effect is not totally described in a relatively coarse (some 1 km) model. Also, it could be linked to an underestimation of the initial deformation at the source (Sahal et al., 2009). Computed travel times have been compared with those obtained from tide gauges and previous models. Results may be seen in Table 3; there is a generally good agreement between both set of data. 4. Model applications Some model applications are described in this section. Risk assessment (Section 1) and spectral analysis of tsunami wave (Section 2) studies are carried out. Finally, the possibility that the capital of Tartessos was destroyed by a large tsunami has been investigated in Section 3. 4.1. Risk assessment The coasts of the Gulf of Cádiz have been exposed to devastating tsunamis in the past, as the 1755 event. Thus, there is a need for quantitative tsunami hazard assessment to support planning in the case of an emergency. This work has already been very carefully carried out in the case of Huelva town, in SW Spain (Lima et al., 2010). We have now extended such work to all SW Spanish shore, south Portugal coasts and Moroccancoasts, which have neverbeen studied before.Historically, tsunamis with origin in the western Mediterranean Sea have had much more smaller amplitudes than those generated in the Gulf of Cádiz. Consequently, the study has been limited to this last region. The five main active faults in theGulf of Cádiz have been selected for this study (Table 1). Tsunamis originating in each of the faults have been simulated and run-ups along the coasts calculated. Asan example, maximum wave amplitudes after 3 h of simulation are presented in Fig. 5 for four of the faults. Ascould be expected, the fault strike is essential, since for a 90° rake thetsunamimaximum amplitudes are produced in the direction perpendicular to the fault (Roger and Hébert, 2008). Thus, in the case of a tsunami originating in HSF, most of the energy is directed towards Moroccan coasts. For PBF and CWF faults, both the south of Portugal and Morocco receive large waves, specially in the case of CWF (see the color scale). A tsunami originating in MPF mainly affects the area of Lisbon and a small portion of the north Moroccan coast. A tsunami produced in GBF (not shown) mainly impacts the Moroccan coast. Nevertheless, the area around Cape San Vicente is, at least, moderately affected by tsunamis generated in all faults. This is due to the fact that, because of bathymetry, this region focuses waves. Fig. 4. Computed maximum wave amplitudes (m) after the 2003 Algeria tsunami. Simulation time 2 h. Table 3 Comparison of calculated arrival times with observations for the 2003 Algeria tsunami. Numbers in parenthesis indicate results from the model by Sahal et al. (2009). The calculated results correspond to the arrival time of a 2 cm amplitude signal. Points are located in Fig. 1. Location Arrival time (min) Calculated arrival time (min) Argel 4 7 Palma 51 53 Ibiza 50 37 Valencia (81) 81 Alicante (70) 69 Málaga (100) 104 Almeria (75) 76 227R. Periáñez, J.M. Abril / Journal of Marine Systems 111–112 (2013) 223–234 Our results are consistent with those of Lima et al. (2010), who found that Huelva town (SW Spain, Fig. 1)issignificantly affected by a tsunami only if its origin is in CWF (see Fig. 5). Run-ups larger than 3 m for the five simulated tsunamis are summarized in Fig. 6, together with the seabed deformation produced by the Okada formulation for each source. It seems that the most dangerous tsunamis, leading to extensive flooding along the Iberian coast, are those originating in CWF. Tsunamis generated in PBF and MPF produce run-ups larger than 3 m along the Portuguese coasts. In the case of Morocco, there is no a significant run-up if the tsunami is generated in MPF, although its south coast is exposed to PBF originated events. CWF tsunamis only affect significantly the south coasts of Morocco. On the other hand, virtually all the Atlantic coast of Morocco is affected by tsunamis generated in GBF and HSF. A risk factor, RF j , has been defined for each location jto summarize all these results. Thus, for a given coastal location: RFj¼∑5 n¼1rnjðÞ Mð7Þ where r n (j) is run-up produced in location jby source nand Mis the maximum of all sums in the numerator: M¼max X 5 n¼1 rnkðÞ "# ð8Þ with the index krunning for all locations as j.RF is normalized, thus it is a number between 0 and 1. An RF closer to 1 indicates that such location is more exposed to be affected by a significant run-up after a potential tsunami originating in any of the active faults than a location with a low RF. Risk factors are presented in Fig. 7. It may be seen that virtually all the coasts of SW Spain, south of Portugal and Morocco may suffer a considerable flooding. However, the Portuguese region to the east of Cape San Vicente is the most potentially exposed to suffer severe flooding. There are also areas with large RF in the south of Morocco and north from Lisbon. It must be commented that run-up results should be interpreted with care, given the relatively low spatial resolution of the model. Nevertheless, they constitute a valuable first order approximation. Fig. 5. Maximum wave amplitudes computed after 3 h of simulation for four of the source parameters listed in Table 1. The color scale provides elevations in m. 228 R. Periáñez, J.M. Abril / Journal of Marine Systems 111–112 (2013) 223–234 Fig. 6. Calculated run-up for each fault. The red circles indicate run-ups larger than 3 m. The color scale provides seabed deformation in m, from Okada formulation. 229R. Periáñez, J.M. Abril / Journal of Marine Systems 111–112 (2013) 223–234 4.2. Spectral analysis It has been shown (Rabinovich, 2009) that the spectral signature of tsunamis depends not only on the fault parameters, but also on the eigenmodes of the medium where they propagate. The Atlantic coasts of SW Spain and Morocco can be roughly seen as forming a semi-circular or a semi-elliptic basin (length ∼210 km and wide ∼400 km), with a semi-parabolic depth profile (∼3100 m at the open boundary). This allows estimating a fundamental period T o ∼5400 s (Rabinovich, 2009), and decreasing values with factors ∼0.7, ∼0.6 and ∼0.5 for the three consecutive harmonic modes (Rabinovich, 2009). A similar approach can be applied to the Huelva– Cádiz arc, in SW Spain, and to the Alboran Sea, leading to fundamental periods of some 4000 s and 10,800 s respectively. From this crude approach, it may be expected that an important fraction of the tsunami energy will be associated to the low frequency spectrum. To test this hypothesis, several synthetic tide gauges have been placed along the coastline as shown in Fig. 1, in addition to gauges in the epicenter of each tsunami and in two fixed control points in the open sea. Time series of water elevations have been then generated with Δt=1 s and length of 3 h. The power spectral density (proportional tothe energy spectraldensity) was then evaluated from FFT analysis (Sahal et al., 2009). Fig. 8 shows the results for a set of four selected locations. 02000 4000 6000 8000 10000 12000 0 200 400 600 800 1000 1200 1400 1600 1800 2000 HUELVA Period (s) Power spectral density HSF GBF PBF LISBON CWF MPF 02000 4000 6000 8000 10000 12000 0 200 400 600 800 1000 1200 1400 CASABLANCA Period (s) Power spectral density HSF GBF PBF LISBON CWF MPF 0 2000 4000 6000 8000 10000 12000 0 500 1000 1500 2000 2500 LISBON Period (s) Power spectral density HSF GBF PBF LISBON CWF MPF 02000 4000 6000 8000 10000 12000 0 10 20 30 40 50 60 70 80 90 100 ALGECIRAS Period (s) Power spectral density HSF GBF PBF LISBON CWF MPF Fig. 8. Computed power spectral density (using FFT) from time series of water elevations(Δt=1 s and length 3 h)at four selectedlocations (see Fig.1) and for six different tsunamis with epicenter in the Atlantic (see Table 1). Note that Algeciras is the name of the bay where Gibraltar is located. −13 −12 −11 −10 −9 −8 −7 −6 −5 33 34 35 36 37 38 39 40 Longitude Latitude RF<0.33 0.33<RF<0.66 RF>0.66 Fig. 7. Computed risk factors. They have been plotted in three intervals for a clearer presentation. 230 R. Periáñez, J.M. Abril / Journal of Marine Systems 111–112 (2013) 223–234 From Fig. 8 (note that not all the coastal locations where the analysis has been carried out are shown in this figure), and as expected, most of the energy appears in the low-frequency spectrum. In some coastal areas (as Gibraltar –located in Algeciras Bay–and Ceuta, see Fig. 1) the tsunami signal is almost “monochromatic”, while in others several main constituents may be distinguished. For each location the main frequencies are excited under most of the tsunami conditions, which suggests that they are linked to eigenmodes. Furthermore, in places with polychromatic spectra, the main frequencies closely follow the harmonic sequence. Thus, in Casablanca the main period appears around 1800 s, and the most energetic tsunamis –as that of 1755–are able to excite the 900 s period. Other examples can be found in Málaga, where the following periods are excited: ∼5400 s, ∼3600 s, ∼2700 s ∼2160 s, ∼1800 s and ∼1350 s (which correspond, respectively, to 1/2, 1/3, 1/4, 1/5, 1/6 and 1/8 of the fundamental period of ∼10800 s estimated from the geometry of the Alboran Sea —according to Rabinovich, 2009). These results have been tested with a simulation time of 10 h. Also we have the cases of Faro, with periods ∼2700 s, ∼1800 s and ∼1200 s; or Lisbon, with main periods around 1800–2160 s and 900–1080 s. 4.3. Tsunami propagation in the former Gulf of Tartessos 4.3.1. Historical context With the last Holocene sea level rise, about 6900 BP, a large marine gulf, the so-called Gulf of Tartessos, extended over the south-western area of the nowadays Guadalquivir Valley. Since then, a series of depositional littoral landforms has been generated, along with the progressive colmation of the lagoon (Rodríguez-Ramírez et al., 1996). This way, 2000 years BP, the area was an inland lagoon, the so-called Lacus Ligustinus, with a double outlet to the sea as described by Estrabon in his Geographica. The kingdom of Tartessos flourished in this area, based upon its mineral resources and in the early trade ties with the Phoenician sailors. The first written references to the kingdom of Tarshish, contained in the Bible, are dated from the tenth century BC. Several human settlements from this period have been identified around the former coastline of the Gulf of Tartessos, being of particular relevance the finding of the Carambolo Treasure, in the vicinity of the city of Seville. By the sixth century BC this culture vanished, perhaps due to the loss of eastern markets for their products. Alternatively, it has been speculated about the possibility of a catastrophic end by some natural hazards; particularly, a tsunami could destroy the main city of this kingdom. Furthermore, it has been suggested that the Atlantis of Plato could refer to the culture of Tartessos (Kühne, 2004). In 1912–1930 the German archeologist Adolf Schulten excavated in the area of Doñana Natural Park, in a fruitless search for the main town of Tartessos. Other authors have suggested that the lost city of Tartessos could have been located in the vicinity of the ancient Caura (nowadays Coria del Río, about 12 km to the south from Seville), or by the mouth of the Guadalquivir River (followers of Adolf Schulten, 1924). Between 218 and 209 BC, this area suffered the impact of a historical tsunami, with epicenter likely located near the Cape San Vicente and producing a run-up of ~5 m (Rodríguez-Vidal et al., 2011). As reported by these authors, this tsunami had a significant morphological effect, eroding the previous littoral spits and foredunes, and generating overwash deposits within the estuary. Earlier tsunamis affected this area, but this is the first one properly documented. Kühne (2004) claimed, based on the analysis of aerial photographs, to have identified at geographical coordinates 36°57′25″±6″and − 6°22′58″±8″a“quadratic”structure of size 280 m times 240 m which could equate to the temple of Cleito and Poseidon, in the lost city of Atlantis–Tartessos. According to the work by Rodríguez- Vidal et al. (2011), this location (depicted in Fig. 9)wasaflooded marshland at ~220 BC. Fig. 9. Computed maximum current magnitude after 3 h of propagation of a 1755-like tsunami (units for the color scale are m/s). The continuous line corresponds to the nowadays coastline, and the main cities have been depicted as a reference, including the one proposed by Kühne (2004) for the lost city of Tartessos. 231R. Periáñez, J.M. Abril / Journal of Marine Systems 111–112 (2013) 223–234