A numerical modeling study on oceanographic conditions in the former Gulf of Tartessos (SW Iberia): Tides and tsunami propagation

A numerical modeling study on oceanographic conditions in the former Gulf of Tartessos (SW Iberia): Tides and tsunami propagation R. Periáñez 1 ,J.M.Abril Dpto. Física Aplicada I, ETSIA, Universidad de Sevilla, Ctra. Utrera km 1, 41013 Sevilla, Spain abstract Keywords: Numerical modeling Iberia Tsunamis Tides Gulf of Tartessos At least five catastrophic tsunami events have affected the Iberian Atlantic coasts during the last 7000 years. During this time, the former Gulf of Tartessos evolved towards the present marshland area, in the lower Guadalquivir valley (SW Spain). Ancient cultures flourished and vanished around this waterbody which, with its tidal dynamics, was an essential part of their living environment. A numerical modeling study on oceano-graphic conditions (tide and tsunami propagation) in the former Gulf of Tartessos has been carried out, in order to improve our insight on how they could have influenced the human activities in this area. The model solves the 2D depth-averaged hydrodynamic equations with appropriate initial and boundary conditions for tides and tsunamis. The model has been tested through simulation of tides under present bathymetry and of past tsunamis for which historical records exist. Then the bathymetry of the Gulf of Tartessos has been reconstructed for 2200 and 4000 years BP. The two main tidal semidiurnal and diurnal constituents have been simulated. Results indicate that they suffer a strong attenuation in the Gulf for both bathymetries, with significant currents only in its entrance. The known main active faults in the Gulf of Cádiz and the 1755 Lisbon source have been selected as case studies of tsunami propagation in the former Gulf of Tartessos. The 1755 Lisbon tsunami represented the worst case situation, and its numerical simulation has been subjected to several sensitivity tests. Results reveal negligible impacts for the inner shoreline of the Gulf of Tartessos, but severe damage could have been produced along the Atlantic coasts of SW Spain, from Cádiz to Huelva. 1. Introduction With the last Holocene sea level rise, about 6900 years BP, a large marine gulf (the Gulf of Tartessos) extended over the south-western area of the present Guadalquivir Valley. Since then, the generation of littoral landforms and the progressive colmation of the inner gulf resulted in a inland lagoon (the so-called Lacus Ligustinus) described by Estrabo in his Geographica. Studies on the geological evolution of this area have been conducted, among others, by Gavala (1959), Menanteau (1984) and Arteaga and Roos (1995).Rodríguez-Ramírez et al. (1996) published a review on the recent coastal evolution of Doñana natural park. The transoceanic tsunami generated by the Lisbon earthquake (November 1st, 1755 AD) is well documented in historical accounts (e.g., the recorded run-up in Cádiz, Spain, was 15 m). But Lario et al. (2011) have identified at least five catastrophic tsunamievents generated by strong earthquakes affecting this area during the last 7000 years, with a recurrence interval between 1200 and 1500 years. Two of them seem to be isochrone with the collapse of ancient cultures in the area. Thus, at 4100 BP, a large extreme wave event was recorded in the Doñana marshlands as an interbedded marine layer within the inner marsh deposits of the Guadalquivir river estuary (Lario et al., 2011). Between ca. 4500 and ca. 4300 BP, the Valencina emporium collapsed. It was the firstknown political system in the Guadalquivir Basin, emerging at the beginning of the third millennium BC, and being centered on thesouth-western Pyrite Belt andtheLowerGuadalquivir Basin (Nocete et al., 2005, 2010). The site of Valencina was its main human settlement, and it became an important trade center for raw materials and products from regional to transcontinental scales (Nocete et al., 2010). Evidences of a widespread catastrophic tsunami have been reported with a radiocarbon data range from 2700 to 2200 cal BP (Lario et al., 2011; Luque, 2008). Thus, it is not clear if they correspond to the same historical documented tsunami of 218–209 BC (Rodríguez-Vidal et al., 2011), or to two distinct events. Particularly, this range includes the collapse of the Tartessian culture around the sixth century BC. It had flourished around the first millennium BC in this area under the Phoenicianinfluence. Thequest of its lost capital, the city of Tartessos described by Avieno, has been one of the most exciting archeological enterprises in the past century. Escacena et al. (1996) reported ancient saltworks in the “Marismillas”site (37.237° N, 6.132° W), located in the inner coastline of the former Gulf of Tartessos and dated at ca. 5000 years BP. In all this time the water body of the Tartessian Gulf, with its tidal dynamics, E-mail address: [email protected] (R. Periáñez). 1 Fax: +34 954486436. became an essential part of the living environment of ancient cultures settled around it. It allowed for an important boat traffic carrying regional and transcontinental trade, fisheries, saltworks, and other human activities. It could have also been the means through which catastrophic tsunamis arrived. A reliable oceanographic reconstruction of this marine palaeoenvironment is then relevant to improve our insight on how it could have influenced the human activities of former cultures. Thus, computed tidal elevationsand currents can provide some insight on theancient trades for ship trafficandfisheries. The simulation of tsunami propagation allows one to estimate their potential hazardous effects on ancient coastal cities. Numerical modeling of tsunami propagation is a well stated methodology which has been validated against recorded data from historical events over the world (Alasset et al., 2006; Choi et al., 2003; Ioualalen et al., 2010; Kuo-Fong and Mon-Feng, 1997; among many others). 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). Some models for predicting tsunami propagation inthe Iberian Atlantic coasts have already been developed (Baptista et al., 2003; Lima et al., 2010). Similarly, numerical modeling of tidal dynamics in marine systems is a well stated discipline (Kowalik and Murty, 1993; Periáñez, 2007; Pugh, 1987). Applications of this kind of models for the Gulf of Cádiz, the Strait of Gibraltar and the Alborán Sea have been conducted, among others, by Tejedor et al. (1999),Periáñez (2007) and Quaresma and Pichon (2013). Periáñez and Abril (2013) developed a model for tsunami propagation in the Iberia-Africa plate boundary and tested it against some recorded historical events, including the 1755 Lisbon tsunami and the ones of 1856 and 2003 originated in the Algerian coasts. That work included a preliminary study, through numerical modeling, on tsunami propagation in the former Gulf of Tartessos as it was at ca. 2200 years BP. An extension of this study was recently published by Abril et al. (2013), in which the methodology was proposed as a tool for Archaeological Science. The work, more focussed in the archeological context, included a second former configuration of the Gulf at ca. 4000 years BP, a basic study on tidal dynamics, and used the 1755 Lisbon tsunami as the worst case situation. Results indicated that a tsunami could hardly have produced any severe damage in the gulf interior, due to the reduction in wave amplitude and water velocity occurring in its entrance. The aim of the present work is conductinga full oceanographic study on numerical simulation of tides and tsunami propagation inthe former Gulf of Tartessos as it was at ca. 4000 and 2200 years BP. The model is described in Section 2; indications on how tsunamis and tides are simulated are given in Sections 2.1 and 2.2 respectively. The tsunami model is fully non-linear and also allows accurate run-up calculations through a flooding and drying algorithm. Model validation against recorded tsunamis was presented in Periáñez and Abril (2013) and is briefly summarized here. Validation of the tidal model has been conducted through the simulation of surface tides under present conditions (i.e., the present bathymetry) and comparing model results against recorded data on water elevations and currents, as reported in detail in Section 3. Some indications on the reconstruction of the Gulf of Tartessos bathymetry for 2200 BP and 4000 BP are presented in Section 4.1. The tidal model is then applied to the former Gulf of Tartessos for the two main semidurnal and diurnal constituents and for both bathymetries. Simulation results are described in Section 4.2. Up to seven different tsunami sources (the known main active faults in the Gulf of Cádiz, and the source of the 1755 Lisbon tsunami) have been used to study the associate extreme wave events in the former Gulf of Tartessos, also handling the two selected bathymetries. The 1755 Lisbon tsunami has been selected as the worst case scenario, and it has been subjected to sensitivity tests for the effect of numerical (spatial and temporal) resolution, friction coefficient in the shallower areas, and bathymetry. Results are described and discussed in Section 4.3. 2. Model description The 2D depth-averaged barotropic shallow water equations, which describe the propagation of surface shallow water gravity waves (see for instance Kowalik and Murty, 1993; Periáñez and Abril, 2013)constitute the basis of the tsunami and tidal models. Simulations of tides and tsunamis essentially differ in initial and boundary conditions, which are briefly described below. 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 s of arc both in longitudeand 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= 0.6 s. Nevertheless, model sensitivity to grid resolution in the Gulf of Tartessos has been studied, as discussed later. 2.1. Tsunami computations The application of the 2D shallow water equations to simulate the propagation of tsunami waves has been described in detail in Periáñez and Abril (2013). Some indications on the methodology are given below. A gravity wave radiation condition (Herzfeld et al., 2011)isusedfor sea surface elevation along open boundaries, which is implemented in an implicit form. This condition avoids artificial wave reflections in open boundaries. A flood/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 adopted. 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; Kuo-Fong and Mon-Feng, 1997; Sahal et al., 2009 among others). The sea-floor deformation produced by the earthquake is computed using the classical Okada formulae (Okada, 1985). 2.2. Tide calculations For simulating tides, water elevations are specified along open boundaries (Schwiderski, 1980a, 1980b) during the simulation period. The model is started from rest. Once that stable oscillations are achieved, standard tidal analysis is carried out (Pugh, 1987). The model includes the two main semidiurnal, M 2 and S 2 , and diurnal, O 1 and K 1 , constituents, as well as the shallow water constituent M 4 . Thus, the hydrodynamic equations are solved for each constituent and tidal analysis is also carried out for each constituent separately. Tidal residual transports are also calculated. This is done for each constituent separately. Tidal residuals for each constituent are calculated from the equation: qr !¼ bHqt !N bHNð1Þ which corresponds to the Eulerian residual transport velocity (Delhez, 1996). In this equation bN is the time averaging operator, qr !is the tidal residual and qt !is the instantaneous tidal current. 3. Model testing 3.1. Tidal computations under present conditions The main semidiurnal tides M 2 and S 2 have been simulated using today bathymetry. Calculated tidal constants (tide amplitudes and phases) are compared with observations available at some locations. Good results are obtained with a bed friction coefficient k= 0.0025 and a horizontal viscosity A=10m 2 /s. A comparison between measured and computed tidal constants for both constituents may be seen in Table 1. Measured values are obtained from NOAA (1982),Tsimplis et al. (1995) and Manzella and Elliott (1991). It may be seen that the model is providing a realistic representation of tidal elevations over the model domain. Tidal currents are also compared with observations, taken from García-Lafuente and Cayo-Lucana (1994),inTable 2.Generally speaking, there is also a goodagreementbetween themeasured and calculated currents for both constituents. Indeed, a regression analysis of computed tide amplitude and current amplitude vs. measured data gives r 2 values of 0.988 and 0.990. 3.2. Simulation of historical tsunamis The tsunami of Lisbon in 1755, as well as Algeria tsunamis in 1856 and 2003, havebeen simulated.Computed run-up, wave maximum amplitude and arrivaltimes have been compared with historical records or previous simulations (when historical data do not exist). Full details have been described in Periáñez and Abril (2013).Asanexample,comparison of model results with historical data is presented in Table 3 for the 1755 Lisbon tsunami. The source is described in Table 4.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 they tend to be overestimated in some locations. Nevertheless, the overall quality of −8 −7.5 −7 −6.5 −6 −5.5 −5 −4.5 −4 35 35.5 36 36.5 37 37.5 38 Longitude Latitude • Seville Gulf of Cadiz • Chipiona Rota • • Cadiz Ayamonte •• Huelva • Gibraltar • A B• • C • D • E • F • G A B Fig. 1. Model domain.Water depths are given in m. The gray box gives the location of the former Gulf of Tartessos. A detailed view of the Gulf of Cádiz is presented in the lower panel, with locations mentioned in the text and Tables. results is similar to that of previous modeling studies (Baptista et al., 2003; Barkan et al., 2009). Results of all simulations mentioned above have indicated that the model is correctly reproducing the propagation of tsunami waves, given the generally good agreement of our results with previous simulations and/or historical records. 4. Model results for the Gulf of Tartessos The bathymetry of the Gulf of Tartessos has been first reconstructed for two situations: 2200 and 4000 years BP. This is very briefly described in the next section. Later, tide and tsunami propagation in the Gulf are presented in Sections 4.2 and 4.3 respectively. A model sensitivity study is finally described in Section 4.3.3. 4.1. Reconstruction of the bathymetry of the Gulf of Tartessos The Gulf of Tartessos was a rather shallow water body at 2200– 4000 years BP. Accretion over the last two millennia roughly accounted for 1 m of sediments, while the accretion of the Doñana spit-barrier was higher by one order of magnitude (Rodríguez-Vidal et al., 2011; Ruiz et al. 2004). The topography of the Gulf as it was at 2200 BP and 4000 BP has been reconstructed following the procedure described in detail in Periáñez and Abril (2013); consequently it is not repeated here. Nevertheless, for self-consistence of the present work, both topographies are presented in Fig. 2. Some present cities are indicated as a reference. These reconstructions are based on previous geological works in the area (Arteaga and Roos, 1995; Gavala, 1959; Rodríguez- Vidal et al., 2011). 4.2. Tides The two main tidal constituents, M 2 and S 2 ,havebeenfirst studied. Calculated tide amplitude, phase and current amplitude for the M 2 constituent arepresented in Fig. 3. There is a relevantattenuation of the tide as entering the Gulf of Tartessos, which is higher for the bathymetry of 2200 BP since the Gulf is closer and shallower than for 4000 BP. As a consequence, energy dissipation is more intense. Even for the 4000 BP bathymetry, tide amplitude is reduced by a factor slightly larger than 2 with respect to the situation in the open sea. Current amplitudes are again smaller for 2200 BP than for 4000 BP, but reaching some 0.6 m/s in the Gulf entrance for both bathymetries. In both cases there is a significant change in the tide elevation phase. The M 2 residual current and surface elevation are presented in Fig. 4 for the 2200 BP situation as an example. Currents inside the Gulf are extremely weak. Only in the connection to the sea they reach values of the order of a few cm per second. The water level is, on average, some 8–10 cm higher in the Gulf than in the sea. This is a typical situation occurring in upper reaches of estuaries (Pugh, 1987),where friction becomes the dominant process: the time taken for water to drain away under gravity becomes longer than the tidal periods themselves, resulting in a long-term pumping-up of water. Residual currents are even weaker for the 4000 BP situation, since the entrance to the Gulf is wider. A higher mean water level than in the open sea is still evident. The behavior on the S 2 tide is similar (Fig. 5). We again find a reduction in tide amplitude inside the Gulf, which is higher for the 2200 BP bathymetry than for the 4000 BP one, as expected. Nevertheless, although the amplitude of the S 2 tide is about 47% of the M 2 one (Pugh, 1987), it is interesting to notice that current amplitudes at the Gulf entrance are not so different for both constituents. Maximum currents are about 0.6 and 0.4 m/s for the M 2 and S 2 tides, respectively (Figs. 3 and 5). Consequently, strong tidal currents, of the order of 1 m/s, would exist during spring tides at the Gulf entrance. Navigation in this area could be difficult to ancient ships during spring tides and strong wind episodes. Asignificant change in the tide elevation phase in the Gulf, with respect to the open sea, is again found for the S 2 constituent. Residual currents are extremely weak even in the Gulf entrance, hardly reaching 1 cm/s in some locations in this area. They are presented,forthe 2200 BP bathymetry, in Fig. 6. Results for the diurnal constituents are presented in Fig. 7.Only the 4000 BP case is shown. It is interesting to note that now a slight amplification of the tide amplitude occurs inside the Gulf. Currents are weak, reaching only some 0.1 m/s in the entrance. The tidal phase change inside the Gulf with respect to the open sea is not so significant as for the semidiurnal constituents. Finally, the M 4 tide has also been simulated for the 4000 BP configuration. The period of this constituent is 6.2 h, comparable to the natural period of oscillation of a semi-enclosed bay with the dimensions of the former Gulf of Tartessos (see Section 4.3.3). This makes possible the Table 2 Computed and observed current ellipse parameters for the M 2 and S 2 tides at points indicated with letters in Fig. 1.Mand θare the magnitude (m/s) of the major semiaxes and its direction (deg), respectively, measured anticlockwise from east. M 2 S 2 Location M obs θ obs M comp θ comp M obs θ obs M comp θ comp Ab0.03 –0.041 125 ––– – B 0.68 10 0.65 15 0.24 15 0.24 11 C 0.25 20 0.40 20 0.12 20 0.14 20 D 0.086 54 0.067 73 ––– – E 1.15 9 0.88 5 0.31 9 0.30 4 F 0.073 22 0.062 1 ––– – G0.011−50.025 −2––– – Table 3 Run-up reports (m) and arrival times after the initial condition (min) from historical data and modeling results for the Lisbon 1755 tsunami. Arrival 1 refers to 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 N10 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 Rabat 33–50 4.7 37 44 Table 1 Established, index obs,andcomputed,indexcomp, amplitudes (A,cm)andphases(g,deg) of tidal elevations at several locations indicated in Fig. 1. M 2 S 2 Location A obs g obs A comp g comp A obs g obs A comp g comp Faro 93 94 99 67 32 125 35 90 Chipiona 102 54 103 62 41 82 39 86 Rota 105 52 103 60 37 78 37 84 Cádiz 100 87 101 59 37 110 37 83 Ayamonte 100 59 101 65 32 88 36 89 Huelva 102 56 105 64 38 82 38 88 Casablanca 99 56 93 54 35 81 36 77 Rabat 88 59 98 58 35 83 37 81 Tarifa 42 41 41 54 14 85 15 79 Ceuta 30 50 25 60 11 76 9.3 88 Malaga 17 56 19 47 7 72 7.2 77 Alhucemas 18 61 18 60 7 80 7.0 89 Almeria 9 51 9.5 48 4 78 4.2 77 Gibraltar 30 46 30 51 11 72 11 78 generation of standing waves in the whole Gulf or in particular subdomains. M 4 charts in the Gulf of Tartessos are presented in Fig. 8. Corange and cotidal lines (Figs. 8 A and C respectively) indicate the existence of a degenerate amphidrome close to the southwest coast of the Gulf. Thus, a pattern of standing waves is effectively generated for this constituent. Tidal currents are weak (Fig. 8B), reaching only some 10 cm/s in the entrance of the Gulf. 4.3. Tsunamis 4.3.1. Tsunami propagation Periáñez and Abril (2013) studied the potential hazardous impacts in the coasts of the Gulf of Cádiz due to tsunamis generated in the five known active faults around the Atlantic area of the Iberia-Africa plate boundary. The propagation of tsunamis generated in these faults inside theGulf of Tartessos has now been investigated. In addition to these five faults, a potential new source described in Matias et al. (2013) has been added (denoted as GC06G), which is located in the Guadalquivir Bank, very close to the entrance of the Tartessian Gulf. Although it has not been possible to identify here a clear large thrust fault (Matias et al., 2013), these authors considered that GC06G represents a tectonic source similar to the well identified Portimao Bank Fault. Given its proximity to the Gulf of Tartessos and, thus, its large potential impact, it has been also included in the study, assigning to it the extreme magnitude of 8.75 Mw. A summary of fault parameters is presented in Table 4. The source for the Lisbon tsunami used in Periáñez and Abril (2013) is also included in the Table. If the present coastline is considered, the highest run-ups along the southwest Iberia coast are obtained for a Lisbon-like tsunami. Thus, this source was considered as representative of the worst case tsunami affecting the former Gulf of Tartessos in the preliminary study in Periáñez and Abril (2013) –for the 2200 BP bathymetry –and in Abril Table 4 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]: Lima et al. (2010);Zitellini et al. (2009);[3]Matias et al. (2013). 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] GBF (Gorringe Bank Fault) −11.281 36.948 137 60 8.3 233 25 [2] HSF (Horseshoe Fault) −10.000 36.133 106 70 10.7 222.1 25 [2] MPF (Marqués de Pombal Fault) −10.067 36.895 86 70 8.0 200.0 25 [2] PBF (Portimao Bank Fault) −8.664 36.105 100 55 7.2 266.3 25 [2] CWF (Cádiz Wedge Fault) −9.329 34.833 133 200 11.1 346.3 6 [2] GC06G (assuming Mw = 8.75) −7.330 36.667 220 120 9.9 −100 35 [3] Fig. 2. Reconstructed bathymetries for the 2200 BP (up) and 4000 BP (down) situations. Color scales are in m. Fig. 3. Computed tidal charts for the M 2 constituent and for the two considered bathymetries. A: tide amplitude (m). B: current amplitude (m/s). C: tide phase (deg). et al. (2013) —4000 BP bathymetry. Results for a Lisbon-like tsunami may be seen in these references and are not repeated here. Instead, results for the six remaining sources (Table 4) are presented. Thus, it will be possible to know if a tsunami generated in any other active fault could have a significant effect inside the Gulf. Results for the 4000 BP bathymetry are presented in Fig. 9 for six tsunami sources. The greatest impact outside the Gulf is obtained for tsunamis generated in CWF and also GC06G, given its proximity to the shore. In these two cases amplitudes exceed 5 m along the coast (color scale is limited to 4 m to appreciate details inside the Gulf). However, amplitudes inside the Gulf only exceed 1 m close the entrance. For the four remaining faults, amplitudes inside the Gulf hardly exceed 0.5 m. Computed amplitudes are even lower for the 2200 BP bathymetry (not shown), given that the Gulf in this case is shallower than in 4000 BP and its connection to the open sea is more limited as well. Nevertheless, energy dissipation is linked to other harmful effects of the tsunami through the morphological changes resulting from the associated strong water currents. Thus, the computed maximum water current amplitude (Fig. 10, 4000 BP) is over 3.0 m/s in the mouth of the Gulf for the GC06G tsunami (the color scale is limited to 3 m/s). Currents over 2 m/s are obtained in the case of CWF and they are in the order of 1 m/s for the remaining faults. Thus, from these figures, the expected morphological impacts of these tsunamis (mainly CWF and GC06G, and, in a less extent HSF, MPF and GBF) would consist of the erosion of the coastal sediments Fig. 4. Residual current for the M 2 constituent in the Gulf of Tartessos for the 2200 BP situation. The color scale provides the mean water surface elevation in m. Fig. 5. Same as Fig. 3 but for the S 2 constituent. Fig. 6. Same as Fig. 4 but for the S 2 constituent. Fig. 7. Computed tidal charts for the O 1 and K 1 constituents and for the 4000 BP bathymetry. A: tide amplitude (m). B: current amplitude (m/s). C: tide phase (deg). and their subsequent spreading on the western area of the former Gulf of Tartessos, according to the work by Rodríguez-Vidal et al. (2011), which describes the effects of the historical tsunami in 218–209 BC. Time series of water elevations for all tsunamis and for a point located in the entrance of the Gulf and a point inside it (TG-1 and TG-4 in Fig. 2) are presented in Fig. 11. These results correspond to the 4000 BP bathymetry, as this is more favorable to tsunami propagation inside the Gulf. It may be appreciated that the highest impact in the Gulf entrance is dueto a Lisbon-like tsunami. A similarwave is produced by the GC06G tsunami, which arrives sooner due to the source proximity to the Gulf. It may be noticed that in this case an amplitude different from zero is apparent at t= 0. This is due to the fact that TG-1 and TG-4 are very close to the GC06G fault. Thus, Okada's equation produces an initial deformation in these points. A clear wave attenuation is apparent inside the Gulf. The first (and highest) wave in TG-4 is about 0.5 m for both the Lisbon and CWF tsunamis. All other faults produce waves below 0.4 m inside the Gulf. It is also interesting to notice that, due to the GC06G fault strike, a rundown is first recorded. It reaches approximately −1mattheGulfentrance. It is worth noting that stratigraphic studies in the mentioned archeological site of the Marismilla (Escacena et al., 1996) did not find any noticeable disturbance, discarding then the occurrence of any extreme wave event during the last 5000 years. Similarly, the tidal regimen in this place located in the inner coastline of the former Gulf of Tartessos was not vigorous enough as to print any signal in the stratigraphy. These findings are in good agreement with present results. 4.3.2. Spectral analysis The spectral signature of tsunami induced waves depends not only on the fault parameters, but also on the eigenmodes of the medium where they propagate (Periáñez and Abril, 2013; Rabinovich, 2009). Thus, 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 (with an almost uniform depth ~3100 m at the open boundary). This allows one to estimate a fundamental period T 0 ~5400s(Rabinovich, 2009). Periáñez and Abril (2013) found this period along with its associated harmonic sequence in the spectral analysis of a synthetic tide Fig. 8. Computed tidal charts for the M 4 constituent and for the 4000 BP bathymetry. A: tide amplitude (m). B: current amplitude (m/s). C: tide phase (deg). Fig. 9. Maximum wave amplitudes (m) in the Gulf of Tartessos for six tsunamis and the 4000 BP configuration after a simulated time of 5 h. Note the different scales in color bars. Fig. 10. Maximum currents (m/s) inthe Gulf of Tartessosfor six tsunamis and the 4000 BP configuration after a simulated time of 5 h. Note the different scales in color bars. gauge placed in Huelva, excited by tsunamis generated by different active faults. The former Gulf of Tartessos, as it was at 2200 BP, can be roughly approached as an almost enclosed basin of rectangular shape, with a length of ~45 km and a mean water depth of ~3.5 m. The fundamental period may be then estimated from the analytical formulae provided by Rabinovich (2009) as ~3.1 × 10 4 s. It corresponds to a quarterwave oscillator, with anodal lineat the open mouth of thebay. The configuration of the Gulf for 4000 BP can be roughly seen as a semi-closed basin of rectangular shape with a length of ~50 km (slightly longer than for 2200 BP), andwith a semi-parabolic bed profile with maximum water depth (at its entrance) of ~7 m. The resulting fundamental period (after Rabinovich, 2009) is ~2.9 × 10 4 s. In Fig. 8, which plots the computed tidal amplitude for the M 4 constituent in the former Gulf of Tartessos at 4000 years BP, the amphidrome corresponds to the nodal line. Thus, the effective position of the open mouth should be slightly displaced towards the inner Gulf. This implies a shorter length and a lower fundamental period, close to that of M 4 (2.2 × 10 4 s). To test these hypotheses, several synthetic tide gauges have been placed within the former Gulf of Tartessos for both configurations (TG-1 to TG-4 in Fig. 2 and TG-5 to TG-7 in Fig. S1, in electronic supplementary material). Five to twelve hour long time series of water elevations have been generated for these gauges with a time resolution of 0.6 s. They have been computed for the tsunamis generated by the Lisbon and the CWF sources under the two bathymetries (chronologies) of the former Gulf of Tartessos (group A), and using only the one of 4000 BP for the rest of tsunamis in Table 4 (group B). It is worth noting that the bathymetry of the Gulf has a negligible effect in the tsunami signal at locations outside it. The power spectral density (proportional to the energy spectral density) was then evaluated from FFT analysis (Sahal et al., 2009). Results for TG-1 to TG-4 are depicted in Fig. 12 and for TG5 to TG7 in Fig. S1 (in electronic suplementary material). Concerning the configuration of 4000 BP, there is a large number of frequencies which are excited by most of the tsunamis, and thus, they can be associated to eigenmodes. In all the cases, the higher determined period is 3.0 ± 1.0 h (solved as 2.5 ± 0.5 h in TG-1 and TG-3) —the range corresponds to half the time window resolved by the applied Matlab FFT algorithm. The most energetic tsunamis, as CWF and GC06G, could excite even higher periods (Fig. S1, in electronic supplementary material), but they cannot be properly determined. The sequence of decreasing eigen-periods is shown in Fig. S2 (electronic supplementary material) for the inner Gulf (TG-2, TG-4, TG-5, TG-6 and TG-7) and comprises (in hours): 1.75 ± 0.25, 1.35 ± 0.15, 1.01 ± 0.08, 0.86 ± 0.06…Considering the resonance found for M 4 (Fig. 8), with a period of 6.2 h, if we take this as the Helmholtz mode (a quarter wave oscillator, n= 0), then the sequence of normalized periods closely follows the theoretical one for a rectangular basin with a semi-parabolic depth profile along its length (Rabinovich, 2009): 1.000, 0.409, 0.259, 0.189…(as shown in Fig. S2). Then, the position of the amphidrome in Fig. 8 delimits the effective open mouth of the inner Gulf, whose periods of free oscillation have been described. Thus, TG-1 is beyond this nodal line, butwe can observe essentially the same eigen-period sequence (Fig. S2). TG-3 is also beyond the nodal line (see Figs. 2 and 8) close to a small island and bay. Thus, the sequence of periods significantly differs from the rest. For the 2200 BP configuration of the former Gulf of Tartessos, with a narrower open mouth, the tsunami signal isless energetic,and only few eigen-periods can be identified (Fig. 12a). They are (in hours): 3.3 ± 1.7, 0.86 ± 0.14, 0.56 ± 0.06 and 0.38 ± 0.04. 4.3.3. Model sensitivity Initially, the model sensitivity to grid spatial resolution has been analyzed. This is a relevant issue since the relatively coarse grid may lead to an artificial wave attenuation inside the Gulf of Tartessos. A sub-domain containing the Gulf of Tartessos was created. It extends from −6.8° to −5.8° in longitude and from 36.55° to 37.55° in latitude. Spatial resolution was increased in factors 2 and 3 in the sub-domain (to 15 and 10 s of arc respectively). Open boundary conditions were obtained from the large domain, for the Lisbon tsunami since this leads to the highest impact. The 4000 BP configuration has been used, as more favorable to tsunami propagation. Time series of water surface elevation in point TG-4 for the three resolutions are presented in Fig. 13. It may be seen that increasing spatial resolution leads to a slight increase of wave amplitude inside the Gulf. However, the maximum difference in wave height (between the original and triple resolution) is 9.5%, and conclusions derived from this study are not changed. 0 50 100 150 200 250 300 −2 −1 0 1 2 3 4 Elevation (m) TG−1 Lisbon CWF GBF HSF PBF MPF GC06G 0 50 100 150 200 250 300 −0.4 −0.2 0 0.2 0.4 0.6 Elevation (m) Time (min) TG−4 Lisbon CWF GBF HSF PBF MPF GC06G Fig. 11. Time series of water elevations at points indicated in Fig. 2 for all simulated tsunamis (including Lisbon) and the 4000 BP bathymetry. The bed friction coefficient was calibrated when tides were calculated for today's bathymetry. Nevertheless, the sensitivity of model output in the very shallow waters of the Gulf of Tartessos to the friction coefficient has been investigated as well. Results are presented in Fig. 13. A reduction of friction produces a wave height increase inside the Gulf. The first wave actually increases from some 0.5 m with the nominal friction to some 0.75 m with friction reduced in a factor 2. Conversely, a bed friction increase leads to a wave height reduction, as may be seen in Fig. 13.Nevertheless, general conclusions of the present study do not change due to slight variations of the bed friction coefficient around its nominal value. Results from both configurations (two chronologies) can also be interpreted as a sensitivity test to small-scale bathymetric settings. These slight bathymetry changes do not modify the main conclusions derived from this work, i.e., tsunamis do not significantly propagate into the Gulf of Tartessos. 5. Summary and conclusions The propagation of tides and tsunamis in the former Gulf of Tartessos (SW Spain) has been studied through numerical modeling. The model is based on the 2D shallow water hydrodynamic equations with appropriate initial and boundary conditionsfor tides and tsunamis. It has been first tested through the simulation of M 2 and S 2 tides for the present day bathymetry. Computed tidal constants and current ellipse parameters are in agreement with observations. Also, the model has been applied to past tsunamis for which historical data or previous simulations exist. Again, our calculations agree with such data and simulations. The former Gulf of Tartessos has been reconstructed for two past situations: as it was 2200 and 4000 years BP. Then tides and tsunamis have been simulated with these two new bathymetries. 0 0.5 1 1.5 2 x 104 0 500 1000 1500 2000 2500 3000 3500 Gauge TG−1 Power spectral density Lisbon 4000 BP Lisbon 2200 BP CWF 4000 BP CWF 2200 BP GC06G 4000 BP 00.5 11.5 2 x 104 0 0.5 1 1.5 2 x 1040 0.5 1 1.5 2 x 104 0 0.5 1 1.5 2 x 1040 0.5 1 1.5 2 x 104 0 0.5 1 1.5 2 x 1040 0.5 1 1.5 2 x 104 0 200 400 600 800 1000 1200 1400 Gauge TG−2 Lisbon 4000 BP Lisbon 2200 BP CWF 4000 BP CWF 2200 BP GC06G 4000 BP 0 200 400 600 800 Gauge TG−3 Period (s) Power spectral density Lisbon 4000 BP Lisbon 2200 BP CWF 4000 BP CWF 2200 BP GC06G 4000 BP 0 50 100 150 Gauge TG−4 Lisbon 4000 BP Lisbon 2200 BP CWF 4000 BP CWF 2200 BP GC06G 4000 BP A) Group A: Lisbon and CWF sources for the two bathymetries 0 10 20 30 40 50 60 70 Gauge TG−1 Power spectral density GBF PBF HSF MPF 0 5 10 15 20 25 30 Gauge TG−2 GBF PBF HSF MPF 0 20 40 60 80 100 Gauge TG−3 Period (s) Power spectral density GBF PBF HSF MPF 0 1 2 3 4 5 6 7Gauge TG−4 Period (s) Period (s) GBF PBF HSF MPF B) Group B: rest of sources for 4000 years BP bathymetry. Fig. 12. Computed power spectral density (using FFT) from time series of water elevations (Δt= 0.6 s and lengths from 5 to 12 h) at selected locations (see Fig. 2).