A Numerical Modelling Study on the Potential Role of Tsunamis in the Biblical Exodus

J. Mar. Sci. Eng. 2015, 3, 745-771; doi:10.3390/jmse3030745 Journal of Marine Science and Engineering ISSN 2077-1312 www.mdpi.com/journal/jmse Article A Numerical Modelling Study on the Potential Role of Tsunamis in the Biblical Exodus José M. Abril † and Raúl Periáñez †,* Departamento de Física Aplicada I, ETSIA, Universidad de Sevilla, Carretera de Utrera km 1, D.P. 41013 Seville, Spain; E-Mail: [email protected] † These authors contributed equally to this work. * Author to whom correspondence should be addressed; E-Mail: [email protected]; Tel.: +34-954-486474; Fax: +34-954-486436. Academic Editors: Valentin Heller and Billy Edge Received: 1 June 2015 / Accepted: 13 July 2015 / Published: 28 July 2015 Abstract: The reliability of the narrative of the Biblical Exodus has been subject of heated debate for decades. Recent archaeological studies seem to provide new insight of the exodus path, and although with a still controversial chronology, the effects of the Minoan Santorini eruption have been proposed as a likely explanation of the biblical plagues. Particularly, it has been suggested that flooding by the associated tsunamis could explain the first plague and the sea parting. Recent modelling studies have shown that Santorini’s tsunami effects were negligible in the eastern Nile Delta, but the released tectonic stress could have triggered local tsunamigenic sources in sequence. This paper is aimed to a quantitative assessment of the potential role of tsunamis in the biblical parting of the sea. Several “best case” scenarios are tested through the application of a numerical model for tsunami propagation that has been previously validated. The former paleogeographic conditions of the eastern Nile Delta have been implemented based upon recent geological studies; and several feasible local sources for tsunamis are proposed. Tsunamis triggered by submarine landslides of 10–30 km3 could have severely impacted the northern Sinai and southern Levantine coasts but with weak effects in the eastern Nile Delta coastline. The lack of noticeable flooding in this area under the most favorable conditions for tsunamis, along with the time sequence of water elevations, OPEN ACCESS J. Mar. Sci. Eng. 2015, 3 746 make difficult to accept them as a plausible and literally explanation of the first plague and of the drowning of the Egyptian army in the surroundings of the former Shi-Hor Lagoon. Keywords: Biblical Exodus; tsunamis; numerical modelling; sea parting; Reed Sea 1. Introduction The Minoan eruption of Thera (Santorini), dated around 1613 BC, was one of the largest Plinian eruptions on earth in the past 10,000 years [1]. Stanley and Sheng (1986) [2] reported the evidence for the presence of ash ejected from this eruption in sediment cores from the eastern Nile Delta, and they suggested that it could be associated with the Exodus plague of darkness. The connection between this extreme volcanic event and the biblical plagues preceding the Exodus, has been largely argued in the scientific literature [3–7]. Recent archaeological findings [8,9] provide new insight on the Exodus route and on the location of the sea crossing. Thus, the site of Migdol (Exodus 14:2; see Figure 1) would equate one of the military fortresses in the Horus Way, placed in the shoreline of the Shi-Hor, a paleo-lagoon opened to the Mediterranean Sea. The water-body that the Israelites crossed when leaving Egypt, called yam suph, the Sea of Reeds, should have been a large body of water in the area of Migdol. Nevertheless, the above scenario would correspond to the context of two centuries after the Minoan Santorini eruption. Indeed, the chronology of the biblical Exodus has been a matter of controversial [10] (see for example the review by Sivertsen, 2009 [4]). This last author proposes that there are two distinct episodes that merged in a single narrative: the exodus-flight (a group of western Semites leaving the Wadi Tumilat area slightly before 1600 BC, during the Hyksos rule in Egypt), and the exodus-expulsion (of Israelite slaves during the co-rulership of Tuthmossis III and Amenophis II, about 1450 BC). Figure 1. Map of the eastern Nile Delta and the Suez Canal, including the main geographical references cited in the text. The zoon-box shows the palaeogeographical reconstruction of the former shoreline at the second millennium BC around the Migdol site and the Shi Hor lagoon (after Hoffmeier, 2005, [9]). J. Mar. Sci. Eng. 2015, 3 747 Sivertsen (2009) [4] had suggested that the first plague would have been the likely result of the Minoan eruption and its associated tsunamis. According to this author, in the exodus-flight episode persisting strong wind conditions would have dried any of the two submerged ridges lying south of Lake Timshah, allowing for the sea crossing, while in the exodus-expulsion, a volcanic eruption in the area of Yali and Nisyros (eastern Aegean Sea) would have produced the tsunami that drowned the Egyptian army at the Shi-Hor Lagoon. The modern reef at 29.88° N extends about 10 km under the Gulf of Suez, and it is not of uniform depth. Voltzinger and Androsov (2003) [11] published a numerical model for the drying of the reef by strong winds using a simplified bathymetry for the Gulf and a shallow reef with a uniform depth of 3 m, getting an exposure time of 4 h under winds of 33 m/s. More recently, Drews and Han (2010) [12] pointed out the strong effect of a non-uniform depth. Thus, they showed that a 3 m deep notch (1 km length) within a reef of 2 m depth would have never been exposed under a wind of 33 m/s. Furthermore, these modelling exercises look for steady-state solutions with typical simulation times of the order of one day; but they omit the tidal oscillations, which are up to 1.5 m range in this area [13,14] and a major factor in the involved physical problem. Thus, this scenario for the sea crossing in the exodus-flight seems quite unlikely. Drews and Han (2010) [12] studied an alternative scenario for a sea parting by strong winds: the Lake of Tanis, northeast of the Shi-Hor palaeo-lagoon. They used the paleogeographic reconstructions by Hoffmeier (2005) [9] in which the Lake of Tanis did not exist and the mouth of the Shi-Hor lagoon (the Kedua passage) was connected to the open sea. Instead, in their modified topography a sand bank delimited the northern extension of Lake Tanis, opened to the sea through the narrow Pelusium mouth. This geometry critically controls the water flux to allow a partial drying under persisting strong winds. Despite of these hypothetical and critical topographic settings, this scenario does not fit well with the exodus path and the location of the Egyptian military fortresses [9]. The idea of a tsunami producing the sea parting in the Exodus has been suggested by several authors [4,15,16], although not supported by quantitative analysis. It is known that the Minoan Santorini eruption generated tsunamis by the entry of pyroclastic flows and by caldera collapse. They left their fingerprints, identified through sedimentary deposits, in most of the coastlines of the Aegean Sea, SW Turkey and Crete [17], Cyprus [18], the coastal area of Israel and Gaza [19–21] and eastern Sicily [22,23]. While model results for these tsunamis are able to reasonably account the estimated runups in the Aegean Sea coasts, they fail to explain the observed effects in other areas (e.g., eastern Sicily and Israel and Gaza coasts); and, particularly, they predict negligible effects in the eastern Nile Delta [17,24,25]. In fact, the Aegean Sea is a semi-enclosed water body, and the islands of the Hellenic arc comprise a physical barrier that dissipates most of the tsunami energy. Moreover, the outer border of the Nile Delta, where a sharp drop in water depth occurs, produces partial reflection and energy dissipation in the tsunami waves (as it can be seen in the dispersion pattern of the Cyprus 1222 earthquake, modelled by Periáñez and Abril, 2014a, [25]). To explain the whole set of tsunamigenic deposits from the Santorini event, Periáñez and Abril (2014a) [25] suggested a scenario of sequential tsunamis linked to intense tectonic stress release. Particularly, a tsunami triggered by a submarine landslide in the eastern area of the Nile Delta would have been able to generate the isochronous tsunamigenic deposits found in Israel and Gaza. Nevertheless, this does not necessarily support the hypothesis of tsunamis flooding large areas in the J. Mar. Sci. Eng. 2015, 3 748 eastern Nile Delta, as required in the Sivertsen hypothesis for the first plague, or in the Jacobovici and Cameron interpretation of the sea crossing. A quantitative analysis is needed, taking into account the former paleogeographic conditions and the potential tsunamigenic sources compatible with the known geological settings of the eastern Nile Delta. This is the aim of the present paper, which tries to improve our insight on whether these extraordinary natural disasters could had been only a source of inspiration for the Exodus narrative or they could have conformed reliable scenarios for the sea parting and the drowning of the Egyptian army (e.g., the heated debate between minimalist and maximalist viewpoints). Going further beyond the previous work [25], the present study includes a more detailed paleogeographical reconstruction of the former coastline, including the target area of the Shi-Hor Lagoon [9]. As potential tsunamigenic sources (Section 2) this work considers submarine landslides in the eastern edge of the stable Nile Delta shelf, defined as large as possible within the known constrains stated by geological studies. A hypothetical and sudden colmation of a branch of the former Damietta canyon has been also considered. In the tree cases sensitivity tests have been conducted for the total volume and the prescribed motion of the landslide, up to conform a set of 14 simulations. Two tsunamigenic scenarios triggered by geological faults are also considered along the Pelussium Mega-Shear line, and the combination of two simultaneous sources. This variety of sources, selected as the “best cases” for the drowning of the Egyptian army in the Shi-Hor area, reveal that due to the physical settings of the former eastern Nile Delta, most of the tsunami energy dissipates within the stable shelf or it travels out. The inner shoreline remains only weakly affected by flooding (Section 3). 2. Materials and Methods 2.1. Model Description Numerical models for tsunami propagation are a relatively well established methodology which has been widely validated against recorded data from historical events over the world [26–29]. For the present study, the selected tsunami propagation model is based on the 2D depth-averaged nonlinear barotropic shallow water equations (see for instance Kowalik and Murty, 1993, [30]): ( ) ( ) 0Hu Hv t x y        (1) 22 22 u u u u u u u v g v A t x y x H x y                         (2) 22 22 v v v v v v u v g u A t x y y H x y                         (3) where u and v are the depth averaged water velocities along the x and y axes, h is the depth of water below the mean sea level, ζ is the displacement of the water surface above the mean sea level measured upwards, H = h + ζ is the total water depth, Ω is the Coriolis parameter (Ω = 2w·sinλ, where w is the Earth rotational angular velocity and λ is latitude), g is acceleration due to gravity, ρ is a mean value of water density and A is the horizontal eddy viscosity. τu and τv are friction stresses which have been written in terms of a quadratic law: J. Mar. Sci. Eng. 2015, 3 749 2 2 2 2 ; uv k u u v k v u v        (4) where k is 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. Horizontal viscosity has been set as A = 10 m2/s and the bed friction coefficient as k = 0.0025. Good model results when estimating wave amplitudes and runups have been obtained with these values. Actually, model results (amplitudes, runups and wave arrival times) have been compared with observations in previous works [28]. Consequently, these values have been retained in the present application. Moreover, the model sensitivity to the bed friction coefficient has been studied in Periáñez and Abril (2014b) [31]. Still waters are used as initial conditions. As boundary conditions, water flow towards a dry grid cell is not allowed, and a gravity wave radiation condition is imposed at the open boundaries [32], which is implemented in an implicit form. Due to the CFL stability condition [30] time step for model integration was fixed as Δt = 2 s. 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) [33] and in Periáñez and Abril (2013) [28] has been adopted. Wet grid cells are defined as those with a total water depth larger than a threshold value typically set as a few centimeters. Dry cells are defined as cells where depth is smaller than the threshold one. Flooding and drying is 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 the threshold depth, which has been set as 10 cm [28]. All the equations are solved using explicit finite difference schemes [30] with second order accuracy. In particular, the MSOU (Monotonic Second Order Upstream) is used for the advective non-linear terms in the momentum equations. For tsunamis produced by earthquakes in geological faults, the vertical sea-floor deformation is considered as the initial condition for the tsunami calculation, and it is computed using the classical Okada formulae [34]. Inputs for this equation are fault plane strike, rake, dip, slip, location, length and width, as well as seismic moment and rigidity. Validation and applications of this model can be found in Periáñez and Abril (2013 [28], 2014a [25], 2014b [31]), and Abril et al. (2013) [35]. The methodology of Harbitz (1992) [36] and Cecioni and Bellotti (2010) [37] has been adopted to simulate the generation of tsunamis by submarine landslides, along with the modification provided by Periáñez and Abril (2014a) [25] to describe varying slopes and asymmetric velocity profiles. In the first stage, the slide can be described as a solid body whose downslope movement locally modifies the bathymetry, resulting in an almost instantaneous and equal change in the level of the overlapping waters, which propagates as gravity waves. The adopted geometry is a box of length L, width B and maximum height Δh, and with an exponential smoothing over a distance S in the front and rear and B/2 on the flanks (see scheme in Figure 2). The resulting volume is V = 0.90B·Δh(L + 0.90S) [36]. As a practical approach, the former slide at its initial position will be superposed onto the present-day bathymetry of the source area. The greatest part of the energy transfer to the water column takes place during the first stage of the slide displacement (approximately until the time of its maximum velocity) [38]. Details of J. Mar. Sci. Eng. 2015, 3 750 the subsequent history (including deformation and breaking up of the slide, and the accurate update of the seafloor bathymetry) will contribute only as second order corrections. Thus, corrections of the present day bathymetry in the depositional areas have been omitted in this modelling approach. Away from the coastal area, changes in bathymetry during the last 4 ka are negligible in this context, and the former paleogeography conditions of the eastern coastline of the Nile Delta have been re-created for this work based upon available geological studies, as explained further below. Figure 2. Sketch with the geometrical parameters defining the submarine landslide (following the formalism proposed by Harbitz, 1992, [36]). The adopted geometry is a box of length L, width B and maximum height Δh, and with an exponential smoothing over a distance S in the front and rear, and B/2 on the flanks. The motion of the slide can be known after solving the governing dynamic equations [39,40], or it can be defined by imposing a prescribed motion based upon a maximum velocity, Umax, and the displacement, R [36]. Here we adopt the approach by Harbitz (1992) [36], in which Umax is estimated as a function of the slope angle, α, the average thickness of the slide, h , its density,  (~1.7 × 103 kg·m−3), the density of turbidity currents, t   (~1.1 × 103 kg·m−3), and the friction (μ) and drag ( u D C ) coefficients: max 1 2 ( ) (sin cos ) t u Dt gh UC      (5) with u D C , the drag coefficient along the upper surface of the slide, being estimated from the roughness length parameter, k (in the range of 0.01 to 0.1 m): 5 2 () 1.89 1.62log u D LS Ck       (6) The value of Umax strongly depends on the estimation of the Coulomb friction coefficient, μ, within an acceptable range (being its upper limit tan st    ), and Umax must remain within the range of the reported values in scientific literature [39,41–43]. In many cases a first and large slope angle is involved in the triggering mechanism. After a displacement R1, the slope angle decreases, but the moving masses still complete a second displacement R2. For each J. Mar. Sci. Eng. 2015, 3 751 slope angle the maximum velocity Umax,1 and Umax,2 are estimated as commented above, and the following function of time is imposed for the slide velocity, vs (Periáñez and Abril, 2014a, [25]): max,1 1 sin( ) 0 s trans t v U t t T     (7) 2 max,2 2 sin[ ( ) ] 22 s trans trans trans T v U t t t t t T        (8) with 1 1 max,1 2 R TU   12 2 max,2 [ ( )] trans R R S t TU     max.2 1 max,1 arcsin trans U T tU     and S(t) is the instantaneous position of the slide front at time t. The “two slope angle” kinematics is a model choice, being more general than the single sinus function used by Harbitz (1992) [36], but containing it as a particular case, and it allows generating asymmetric velocities profiles (as the ones used by Lastras et al., 2005, [39]). Applications of this model can be found in Periáñez and Abril (2014a) [25] and in Abril and Periáñez (2015) [38]. Friction stresses over the moving slide are formulated in terms of relative speed, as in Harbitz (1992) [36], in such a way that a slice moving faster than the water column can transfer energy to it. The slide model requires four input (but not free) parameters: maximum velocities (governed by μ) and displacements. The candidate source areas will be characterized by depth profiles along their respective transects, what allows the identification of slope angles. The displacements can be partially estimated from the previous profiles, or introduced as plausible values. These displacements have to be understood as effective run-out distances (displacement of the sliding block), over which the transfer of energy to the tsunami takes place. More details will be provided along with the source definitions. A further validation for the performance of the numerical model with its flood-drying algorithm can be found in Periáñez and Abril (2015) [44]. The submarine landslide submodel has been tested against independent modelling works and subject to extensive sensitivity tests [38,44]. The effects of model resolution and sensitivity to the friction coefficient have been studied in Periáñez and Abril (2014b) [31]. 2.2. Model Domain and Bathymetric Map of the Former Eastern Nile Delta The computational domain used for the numerical modelling of tsunami propagation comprises from 27.0° E to 36.0° E in longitude and from 30.5° N to 37.3° N in latitude (see partial views in Figure 3). Water depths have been obtained from the GEBCO08 digital atlas, available on-line, with a resolution of 30 s of arc both in longitude and latitude. The emerged lands appear with a minimum elevation of 1 m, the minimum depth for waters is 1 m, and both are provided with a resolution of 1 m. J. Mar. Sci. Eng. 2015, 3 752 Figure 3. Details of the domain used for the tsunami propagation model: General view of the former Nile Delta and the Levantine coastlines (up); detailed view of the eastern Nile Delta and the Shi-Hor lagoon (bottom). Water depths (m) with 30 s of arc resolution from GEBCO08 bathymetry and paleogeographic reconstructions (see text) are drawn. Black line is the present coastline, the red and green ones are, respectively, the former coastline and the marshland limit at 3500 year BP after Coutellier and Stanly (1987) [45]. The blue line delimits wetlands at northern Ballah Lake. Tjaru (Hebua I) and Migdol are Egyptian military fortresses around the former Shi-Hor paleolagoon [9]. Dots 1 to 6 are synthetic gauges. The locations of tsunamigenic sources (Tables 1 and 2) are also depicted: yellow boxes for submarine landslides and red lines for faults. During the last 5 ka, with a stable sea level, the major changes in the coastline of the Nile Delta have been driven by fluvial deposits, coastal erosion and land subsidence. The last is estimated as 2.0 mm/year in the area of Alexandria [46]. Coutellier and Stanley (1987) [45] estimated mean Holocene sedimentation rates of 5.0 m/ka in the eastern Nile Delta. As a result, the coastal margin in the studied area migrated northward some 50 km during the past 5 ka. As a practical and crude approach for the reconstruction of the former coastline in the entire Nile Delta at ca 3500 years BP, we followed the J. Mar. Sci. Eng. 2015, 3 753 methodology described by Periáñez and Abril (2014a) [25]. Briefly, we applied a correction to the present depths/elevations, linearly decreasing from the modern to the former coastline (this last as defined by Coutellier and Stanley, 1987, [45]). The small-scale details of the coastline will not significantly affect the general patterns of tsunami propagation, but they are of capital importance for assessing local run-ups. The detailed paelogeographical reconstruction published by Hoffmeier (2005) [9] for the area of the Shi-Hor Lagoon has been digitalized and merged with the previous corrected bathymetry, applying then a filter for smoothing the bathymetry among adjacent grid-cells. Results are shown in Figure 3. 2.3. Tsunamingenic Sources in the Eastern Nile Delta Works and reviews on the geodynamics of the eastern Mediterranean can be found, among others, in Barka and Reilinger (1997) [47], Mascle et al., (2000) [48], or Yolsal and Taymaz (2012) [29]. Its complex tectonic is responsible for intense earthquake and volcanism activity, triggering large tsunamis in the past [17,29,49,50]. The study by Papadopoulos et al., (2007) [51] on tsunami hazards in the Eastern Mediterranean concluded that the mean recurrence of strong tsunamis (most of them were caused by earthquakes) is likely equal to about 142 years. The assessment of the probability for landslide tsunami scenarios is a quite difficult task, due to the insufficient statistics from the recorded events, as discussed in detail by Harbitz et al., (2014) [52]. Such probabilities of occurrence must not be confused with the probability of that a given coastal area is affected by a tsunami of certain level (e.g., with a maximum water height). The first map of the Nile deep-sea fan has not been completed until recent years [53,54]. The Messinian salinity crisis led to the deposition of salt and anhydrite throughout the Mediterranean basin, while the proto-Nile river excavated deep canyons and transported offshore large volumes of terrigenous sediments. Thick Plio-Quaternary sediments covered then the ductile evaporitic layers, what triggered some giant gravity-driven salt tectonics [54–58]. A series of fault trends have been identified across the stable shelf of the Nile Delta, e.g., Rosetta, Baltim and Temsah fault trends, and the Pelusium Mega- Shear Fault system [59–61]. The Nile Delta is a major gas and condensate province with several mud volcanoes and geodynamics associated to fluid seepage [62,63]. A general overview on multi-scale slope instabilities along the Nile deep-sea fan can be found in Loncke et al., (2009) [57], and in Urgeles and Camerlenghi (2013) [64]. Garziglia et al., (2008) [60] identified and dated seven mass-transport deposits on the Western province, northern to the Rosetta Canyon, with volumes ranging from 3 to 500 km3, mean thickness from 11 to 77 m and run-out distances from 18 to 150 km, being the youngest deposit older than 8940 ± 30 cal. year BP. Ducassou et al., (2009), based upon 42 sediment cores collected across the entire Nile deep sea fan, identified several slump deposits and turbidities in the last 2000 ka BP [55], but none in the Mid and Late Holocene (ca 5 ka to present). Thus, any candidate source area for submarine landslides must be compatible with the “empty spaces” within this cloud of cores. Recently, Ducassou et al., (2013) [56] reported four highly mobile debris flows in the Nile deep-sea fan system, with a chronology confined between 5599 and 6915 cal. years BP for the most recent event, in the Rosetta province. Results from numerical models have proved that the effects of the Minoan Santorini tsunamis (triggered by entry of pyroclastic flows and caldera collapse) were negligible in north-eastern Egypt and J. Mar. Sci. Eng. 2015, 3 760 Figure 6. Computed time series of water elevations at tidal gauges TG-3, TG-4 and TG-5 for the submarine landslide SL-2, with runs R1 to R4; and CSL, with runs R1 to R3 (Table 1). See Figure 3 for the location of the synthetic gauges. Initial water depths were 9.0, 2.1 and 8.0 m for TG-3, TG-4 and TG-5, respectively. For tsunamis triggered by geological faults the maximum energy corresponds to the initial potential energy due to the deformation of the free water surface. It was 4.6 × 1013 J for F-1and 1.8 × 1014 J for F-2. The two studied sources are less energetic events than the submarine landslides. Figure 7 shows the maximum computed wave amplitude for tsunamis F-1, F-2, and for the simultaneous source F-2+SL-2 R2 (a landslide triggered in sequence with the earthquake). In this last case, the landslide component dominates the propagation pattern. In tsunamis F-1 and F-2, the highest amplitudes remains confided within the stable shelf of the Nile Delta, and they show strong directionality towards the coastline J. Mar. Sci. Eng. 2015, 3 761 delimiting the Bardawil Lake, but with weak effects around the Shi-Hor lagoon, as shown in the computed time series of Figure 8. Figure 7. Computed maximum amplitude for water elevations (m) due to geological faults F-1, F-2 and the double source F-2 + SL-2 R2. Source parameters are presented in Table 1. Simulation time is 5 h. J. Mar. Sci. Eng. 2015, 3 762 Figure 8. Computed time series of water elevations at tidal gauges TG-3, TG-4 and TG-5 for the geological faults F-1. F-2, and the double source F-2 + SL-2 R2 (Table 1). See Figure 3 for the location of the synthetic gauges. Initial water depths were 9.0, 2.1 and 8.0 m for TG-3, TG-4 and TG-5, respectively. Figure A2 shows the computed time series of water elevations at the synthetic gauges TG-1, TG-2 and TG-6 (see Figure 3) for a sub-set of tsunamis. Tsunamis triggered by geological faults (F-1 and F-2) have negligible effects in the Levantine coasts and in Alexandria (TG-6), in the western Nile Delta. The most energetic submarine landslides impact the Levantine coasts with waves above 15 m height, while their effects in the western Nile Delta are much weaker, with wave amplitudes of 0.5 m. For the studied landslide tsunamis, water currents exceed 5 to 10 m/s along the eastern delta shelf edge and the southern Levantine coastline (Figure 9 shows some examples). These are the areas where the major energy dissipation (Equations (2)–(4)) and, consequently, the greatest impacts occur. In all the cases, wave amplitudes only slightly exceed 1–1.5 m for the most energetic tsunamis at the Shi-Hor Lagoon (synthetic gauges 4 in Figures 5 and 6). The paleo-lagoon never dries out, nor the Kedua passage. Thus, around Shi-Hor Lagoon, only the coastal areas less than 1.0–1.5 m above sea level would have had some risk of being flooded. These dynamical tsunami scenarios (Figures 5 and 6), and particularly the time sequence of water elevations, perhaps could produce some scene of fear and chaos J. Mar. Sci. Eng. 2015, 3 763 in a military group camped at the Tjaru’s seaside, but in our opinion they hardly can explain a sudden drowning of the Egyptian army. Figure 9. Computed maximum water currents (m/s) for tsunamis SL-1 R6, SL-2 R1 and CSL R1. Source parameters are presented in Table 1. Simulation time is 5 h. J. Mar. Sci. Eng. 2015, 3 764 Figure 10. 3D-view of the free water surface deformation, computed after 10 min of the beginning of landslides SL-1 R1 (up) and SL-2 R2 (bottom). The above results are a consequence of the particular geological setting of the studied marine system. A large area of the Nile Delta shelf extends in front of the Shi Hor Lagoon. In these shallow waters, tsunami waves travels at low speed, while the bed friction is high (it is inversely proportional to water depth). The implemented formulation for friction (Equations (2)–(4)) has been widely used and validated in many oceanographic studies, and particularly in models for tidal and tsunami propagation (see for instance Periáñez and Abril, 2014b, [31]). Figure A3 shows a sensitivity test for the friction coefficient, applied to the most energetic tsunami from our set of simulations (SL-1 R6). When the nominal value of the friction coefficient (k = 0.0025) is increased / decreased by a 50%, water elevations at the gauges TG-3 to TG-5 slightly decrease/increase, without affecting in the essential to the main conclusions of this study. Scenarios for severe flooding can be found along the northern Sinai coasts, as shown in Figure 4. Nevertheless, they do not match with the plausible paths for the exodus. It is worth noting the great exposure to most of the tsunamis of the sand barrier delimiting the present Bardawil Lagoon. It is known that it was a military path in Persian times, and it has been proposed by some authors as the exodus path, J. Mar. Sci. Eng. 2015, 3 765 although the most recent studies concluded that this sand barrier was not a passable path at the second millennium BC [9]. Tsunamis triggered by submarine landslides are characterized by a depression of the free water surface over the line of rupture while a huge wave appears at the front of the moving slice. For some exceptional witnesses on the northern Sinai shoreline or in cabotage sailing in this area, a submarine landslide along the eastern Nile Delta shelf edge could have been seen as a parting of the sea (as shown in the 3D-view of Figure 10). The fingerprint of noticeable tsunami impacts in coastal areas can be potentially found in their characteristic sedimentary deposits, although they are not always well preserved; and huge submarine landslides can be detected and dated through their associated turbidite deposits. The selected tsunamigenic sources in this paper are compatible with the main geological features of the Nile Delta, and they are not in contradiction with our present knowledge from field studies. Obviously, this is not a proof for their feasibility, what can only be provided by empirical evidences. They represent extreme and near-field cases in an attempt to provide the best scenarios for the maximalist hypothesis of tsunamis behind the biblical Exodus. But, even for these favourable scenarios, the lack of noticeable flooding of the area, along with the time sequence of water elevations, make difficult to accept a tsunami as a reliable hypothesis neither for the first plague nor for the drowning of the Egyptian army in the surroundings of the Shi-Hor Lagoon. The tsunamigenic deposits from the Santorini event in the coastal area of Israel and Gaza [19] can be explained by the occurrence of tsunamis with source in the eastern Mediterranean, likely triggered by submarine landslides. In our opinion this could have been a source of inspiration for the Exodus narrative, which has to be understood as the likely result of the merging of several stories from different sources [4]. 4. Conclusions A numerical model for tsunami propagation, based on the 2D shallow water equations and previously applied and validated in a wide set of scenarios, has been adapted, from available geological and archaeological studies, to the former paleogeographic conditions in the eastern Nile Delta. Seventeen potential tsunamigenic sources in the eastern Nile Delta have been suggested as possible candidates within the Minoan tsunami sequence, comprising faults along the Pelusium Mega-shear system and several submarine landslides at the shelf edge and along the hypothetical remains of the Messinian Damieta canyon. They have been defined as large as possible within the known constrains stated by geological studies, and comprise sensitivity tests for the total volume and the prescribed motion of the landslides. Tsunamis triggered by submarine landslides could have severely impacted the northern Sinai and southern Levantine coasts with waves over 10–15 m height and water currents over 5–10 m/s, but with weak effects in the surroundings of the Shi-Hor Lagoon, the proposed scenario for the Biblical sea parting. The real occurrence of tsunamis can only be proved by empirical evidences; but as the proposed ones are representative of the most favorable scenarios for the hypothesis to be tested, the lack of any noticeable flooding, and the particular time sequence of water elevations, make difficult to accept them J. Mar. Sci. Eng. 2015, 3 766 as plausible and literally explanation for the first plague and for the drowning of the Egyptian army in the surroundings of the Shi-Hor Lagoon, although they might have been a source of inspiration for the Biblical narrative. Supplementary Materials Snapshots for the propagation of the tsunami SL-2 R2. Scale in the color bar is given in m. Source parameters are presented in Table 1. Supplementary materials can be accessed at: http://www.mdpi.com/2077-1312/3/3/745/s1. Author Contributions Conceived and designed the experiments: JMA, RP. Performed the experiments: JMA, RP. Analyzed the data: JMA, RP. Wrote the paper: JMA, RP. Conflicts of Interest The authors declare no conflict of interest. Appendix Electronic Supplementary Material Figure A1. Prescribed time series for the slice velocity (Equations (7)–(8)) for some of the submarine landslide scenarios defined in Table 1. J. Mar. Sci. Eng. 2015, 3 767 Figure A2. Computed time series of water elevations at tidal gauges TG-1, TG-2 and TG-6 for six of the tsunamigenic sources given in Tables 1 and 2. See Figure 3 for the location of the synthetic gauges. Initial water depths were 36, 33 and 36 m for TG-1, TG-2 and TG-6, respectively. Figure A3. Computed time series of water elevations at tidal gauges TG-3, TG-4 and TG-5 for tsunami SL-1 R6 by increasing/decreasing by a 50% the nominal value of the friction coefficient (Equations (2)–(4)). See Figure 3 for the location of the synthetic gauges. Initial water depths were 9.0, 2.1 and 8.0 m for TG-3, TG-4 and TG-5, respectively. J. Mar. Sci. Eng. 2015, 3 768 References 1. Friedrich, W.L. The Minoan Eruption of Santorini around 1613 B.C. and its consequences. Tagung en des Landesmuseums für Vorges chichte Halle. Band 2013, 9, 37–48. 2. Stanley, J.D.; Sheng, H. Volcanic shards from Santorini (Upper Minoan ash) in the Nile Delta, Egypt. Nature 1986, 320, 733–735. 3. Bruins, H.J.; van der Plicht, J. The Exodus enigma. Nature 1996, 382, 213–214. 4. Sivertsen, B.J. The Parting of the Sea. How volcanoes, Earthquakes, and Plagues Shaped the Story of Exodus; Princeton University Press: Princeton, NJ, USA, 2009. 5. Trevisanato, S.I. Treatments for burns in the London Medical Papyrus show the first seven biblical plagues of Egypt are coherent with Santorini’s volcanic fallout. Med. Hypotheses 2006, 66, 193–196. 6. Trevisanato, S.I. Six medical papyri describe the effects of Santorini’s volcanic ash, and provide Egyptian parallels to the so-called biblical plagues. Med. Hypotheses 2006, 67, 187–190. 7. Trevisanato, S.I. Medical papyri describe the effects of the Santorini eruption on human health, and date the eruption to August 1603–March 1601 BC. Med. Hypotheses 2007, 68, 446–449. 8. Byers, G.A. New evidence from Egypt on the location of the Exodus Sea crossing, part 1. Bible Spade 2006, 19, 14–22. 9. Hoffmeier, J.K. Ancient Israel in Sinai. The Evidence for the Authenticity of the Wilderness Tradition; Oxford University Press Inc.: New York, NY, USA, 2005. 10. Wood, B.G. Debunking “The Exodus Decoded”. ABR, 20 September 2006. Available online: Http://www.biblearchaeology.org/post/2006/09/20/Debunking-The-Exodus-Decoded.aspx (accessed on 15 May 2015). 11. Voltzinger, N.; Androsov, A. Modeling the Hydrodynamic Situation of the Exodus. Izv. Atmos. Ocean. Phys. 2003, 39, 482–496. 12. Drews, C.; Han, W. Dynamics of Wind Setdown at Suez and the Eastern Nile Delta. PLoS ONE 2010, 5, e12481, doi:10.1371/journal.pone.0012481. 13. Abril, J.M.; Abdel-Aal, M.M. Marine Radioactivity Studies in the Suez Canal, Part I: Hydrodynamics and Transit Times. Estuar. Coast. Shelf Sci. 2000, 50, 489–502. 14. Abril, J.M.; Abdel-Aal, M.M. Marine Radioactivity Studies in the Suez Canal, Part II: Field Experiments and a Modelling Study of Dispersion. Estuar. Coast. Shelf Sci. 2000, 50, 503–514. 15. Jacobovici, S.; Cameron, J. The Exodus Decoded; Documentary Film: Toronto, ON, Canada, 16 April 2006. 16. Salzman, R. Mega-Tsunami: The True Story of the Hebrew Exodus from Egypt; iUniverse: Lincoln, NE, USA, 2005. 17. Novikova, T.; Papadopoulos, G.A.; McCoy, F.W. Modelling of tsunami generated by the giant Late Bronze Age eruption of Thera, South Aegean Sea, Greece. Geophys. J. Int. 2011, 186, 665–680. 18. Mészáros, S. Some words on theMinoan tsunami of Santorini. In Thera and the Aegean World II; Doumas, C., Ed.; The Thera Foundation: London, UK, 1978; pp. 257–262. 19. Goodman-Tchernov, B.N.; Dey, H.W.; Reinhardt, E.G.; McCoy, F.; Mart, Y. Tsunami waves generated by the Santorini eruption reached Eastern Mediterranean shores. Geology 2009, 37, 943–946. J. Mar. Sci. Eng. 2015, 3 769 20. Neev, D.; Balker, N.; Emery, K.O. Mediterranean Coast of Israel and Sinai, Holocene Tectonics from Geology and Geophysics and Archaeology; Taylor and Francis: New York, NY, USA, 1987; p. 130. 21. Pfannenstiel, M. Erlauterungen zu den bathymetrischen Karten des ostlichen Mittelmeeres. Bull. Inst. Oceanogr. 1960, 57, 60. 22. De Martini, P.M.; Barbano, M.S.; Smedile, A.; Gerardi, F.; Pantosti, D.; del Carlo, P.; Pirrotta, C. A unique 4000 year long geological record of multiple tsunami inundations in the Augusta Bay (eastern Sicily, Italy). Mar. Geol. 2010, 276, 42–57. 23. De Martini, P.M.; Barbano, M.S.; Pantosti, D.; Smedile, A.; Pirrotta, C.; del Carlo, P.; Pinzi, S. Geological evidence for paleotsunamis along eastern Sicily (Italy): An overview. Nat. Hazards Earth Syst. Sci. 2012, 12, 2569–2580. 24. Ozel, N.M.; Ocal, N.; Cevdet, Y.A.; Dogan, K.; Mustafa, E. Tsunami hazard in the Eastern Mediterranean and its connected seas: Towards a tsunami warning center in Turkey. Soil Dyn. Earthq. Eng. 2011, 31, 598–610. 25. Periáñez, R.; Abril, J.M. Modelling tsunamis in the Eastern Mediterranean Sea. Application to the Minoan Santorini tsunami sequence as a potential scenario for the biblical Exodus. J. Mar. Syst. 2014, 139, 91–102. 26. Choi, B.H.; Pelinovsky, E.; Kim, K.O.; Lee, J.S. Simulation of the trans-oceanic tsunami propagation due to the 1883 Krakatau volcanic eruption. Nat. Hazards Earth Syst. Sci. 2003, 3, 321–332. 27. Kowalik, Z.; Knight, W.; Logan, T.; Whitmore, P. Numerical modelling of the Indian Ocean tsunami. In The Indian Ocean Tsunami; Murty, T.S., Aswathanarayana, U., Nipurama, N., Eds.; Taylor and Francis Group: London, UK, 2007; pp. 97–122. 28. Periáñez, R.; Abril, J.M. Modeling tsunami propagation in the Iberia–Africa plate boundary: Historical events, regional exposure and the case-study of the former Gulf of Tartessos. J. Mar. Syst. 2013, 111–112, 223–234. 29. Yolsal-Cevikbilen, S.; Taymaz, T. Earthquake source parameters along the Hellenic subduction zone and numerical simulations of historical tsunamis in the Eastern Mediterranean. Tectonophysics 2012, 536–537, 61–100. 30. Kowalik, Z.; Murty, T.S. Numerical Modelling of Ocean Dynamics; World Scientific: Singapore, 1993. 31. Periáñez, R.; Abril, J.M. A numerical modeling study on oceanographic conditions in the former Gulf of Tartessos (SW Iberia): Tides and tsunami propagation. J. Mar. Syst. 2014, 139, 68–78. 32. Herzfeld, M.; Schmidt, M.; Griffies, S.M.; Liang, Z. Realistic test cases for limited area ocean modelling. Ocean Model. 2011, 37, 1–34. 33. Kampf, J. Ocean Modelling for Beginners; Springer-Verlag: Heidelberg, Germany, 2009. 34. Okada, Y. Surface deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 1985, 75, 1135–1154. 35. Abril, J.M.; Periáñez, R.; Escacena, J.L. Modeling tides and tsunami propagation in the former Gulf of Tartessos, as a tool for Archaeological Science. J. Archaeol. Sci. 2013, 40, 4499–4508. 36. Harbitz, C.B. Model simulations of tsunamis generated by the Storegga slides. Mar. Geol. 1992, 105, 1–21. 37. Cecioni, C.; Bellotti, G. Modeling tsunamis generated by submerged landslides using depth integrated equations. Appl. Ocean Res. 2010, 32, 343–350.