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

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 km 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.


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][4][5][6][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). ).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][20][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 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).

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][27][28][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]): 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: 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 m 2 /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 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.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 × 10 3 kg•m −3 ), the density of turbidity currents, t   (~1.1 × 10 3 kg•m −3 ), and the friction (μ) and drag ( 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): 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][42][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 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 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].

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.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 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.

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][55][56][57][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][60][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 km 3 , 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 Levantine coasts [17,25], as already commented.Thus, for any less energetic event in the inner Aegean Sea (as the volcanic eruption in the area of Nisyros and Yale suggested by Sivertsen for the sea parting in the exodus-expulsion) a similar behaviour would be expected.Even if potential tsunami directionality is invoked, the Isle of Rhodes prevents any direct pathway for energy transfer from these sources towards the eastern Nile Delta.Moreover, for the Cyprus AD 1222 earthquake tsunami, with source area in southern Cyprus, model results (Periáñez and Abril, 2014a, [25]) show that the outer shelf of the Nile Delta acts as a natural barrier and slows waves at these shallow depths (<500 m), thus increasing wave amplitude in this area, but with negligible impacts in the inner shoreline.
Recent mass-wasting events can be recognized in bathymetric maps by the presence of head and footwall scars [57,64].Periáñez and Abril (2014a) [25] suggested several hypothetical candidate sites in the Nile Delta for submarine landslides of 9-10 km 3 , accomplishing for high slopes and being distant enough from the already studied areas, where any mass-wasting deposits in the last 5 ka can be discarded.They showed that tsunamis generated by sources in the western and northern Nile Delta did not significantly affect the coastal zones of Israel and Gaza, and their effects on the eastern area of the delta were equally negligible.Thus, the source area within the Nile Delta able to produce tsunamis potentially linked to the Exodus has to be confined to its eastern zone.
In this work the geometry of the moving boxes and their displacement will be defined with the criteria of being generously large, but compatible with the "empty areas" defined by the studied sediment cores.The value of μ will be subject to sensitivity tests; and 50 m/s has been adopted as the upper limit for Umax, according to Harbitz (1992) [36].The tsunamigenic sources selected for this study are summarized in Tables 1 and 2 and briefly discussed further below.
Landslide SL-1 is the one discussed and modelled by Periáñez and Abril (2014a) [25], but applied here along with the more detailed paleogeography.The slide front is placed at the eastern border of the stable Nile delta, facing a down-slope of 2.8 degrees.The second displacement prevents reaching the sites of cores studied by Ducassou et al. (2009) [55], although its numerical value has a minor effect in the tsunamigenic potential of the slide.The slide extends over a large area (B = 80 km; L = 20 km), but with a moderate value for its maximum height (6.0 m), accounting for a total volume of 9.80 km 3 .The value of μ has been fixed as 0.1, 0.3, 0.5, 0.8 and 0.9 of its maximum (static, μst) limit in model runs R1, R2, R3, R4 and R5, respectively.A second version of this slide uses a larger value for its maximum height (20.0 m), leading to a total volume of 32.7 km 3 .The larger height increases the slide speed, which reaches values up to 49.5 m/s for μ = 0.5μst (run R6), and of 31.3 m/s for μ = 0.8μst (run R7).The displacement of the slide over the second slope, with a smaller angle, has a minor contribution to its tsunamigenic potential; thus, and for the sake of simplicity, a value of μ = 0.75μst has been adopted for runs R1 to R7. 0.7 37.9 37.9 R3 0.9 21.9 21.9 # Defined as in Harbitz (1992) [36].Slide volume V = 0.9Bhm(L + 0.9S); B, width; L, length; S, smoothing distance; hm, average height; $ From the positive X direction (West to East); ¶ R1 and R2 are down-slope displacements, with slope angles α1 and α2 and maximum speeds Umax,1 and Umax,2, respectively.For SL-1 and SL-2 slides, µ/µst = 0.75 for the second slope; while a single slope scheme is adopted for CSL slides; * Implemented along with a change in bathymetry defining an hypothetical canyon lying E-W, 62 km length, 6.5 km wide and 35 m deep, excavated on the reconstructed 3500 year BP bathymetry, and centred at 31.27° E, 31.45°N. Slide SL-2 is also extracted from the present bathymetry.Its front (Table 1), at the south-eastern border of the stable shelf in the Nile Delta, is initially at 71 m water depth and displaces 4.3 km down a slope with angle 2.8°, and then 11.7 km along a slope of 0.6°.The displacement ends at a plateau.The area covered by this hypothetical mass-wasting does not disagree with the geological surveys above referred.With a maximum height of 26.0 m and a total volume of 10 km 3 , the moving slide reaches a maximum speed of 47.7 m/s for μ = 0.6μst (SL-2, run R1).Runs R2 to R4 for this slide use increasing values of μ in the first displacement, while the same criteria than for SL-1 has been adopted for the second displacement (Table 1).Loncke et al., (2006) [54] reported (their Figure 9) the paths of the Messinian canyons in the Nile Delta.Particularly, at the Damietta branch, a former canyon of some 60 km length and 6 km wide ran eastward.In the GEBCO08 bathymetry, the stable shelf area shows some holes of several hundred meters deep and few km wide (dark spots in Figure 3), likely remains of the former canyons.Due to its proximity to the former Shi-Hor Lagoon, and as a potential and extreme case, the source CSL (Table 1) simulates the effect of a large submarine landslide running into a still partially colmated canyon.The box geometry has been selected as large as possible, and its displacement has to be confined within the width of the former canyon.A maximum height of 20.4 m has been adopted to account for a total volume of 10.0 km 3 .For simplicity a uniform slope angle of 3.0 degrees has been selected, with three options for µ values, leading to maximum speeds of 48.9, 37.9 and 21.9 m/s in runs CSL R1, R2 and R3, respectively.
Gamal (2013) [59] provided a detailed study of the Pelusium Mega-Shear Fault system.The Pelusium line runs north-westwards crossing the south-eastern Nile Delta.Hypothetical earthquakes triggered by this fault with epicentres in this area of the Nile Delta could have been tsunamigenic sources whose potential effects on the nearby former shoreline of the Shi-Hor Lagoon will be studied here.Thus, sources F-1 and F-2 (Table 2) are defined with tentative values for length, width and slip, and angular parameters inspired in the studies of Badawy and Abdel Fattath (2001) [65], and Gamal (2013) [59].As the tsunamis triggered by geological faults are less energetic events when compared with the studied submarine landslides, the application of sensitivity tests has been discarded in this case.
Finally, the case of a fault-earthquake that triggers an almost simultaneous submarine landslide is also considered by combining the two most energetic sources, F-2 and SL-2 R2.

Results and Discussion
The numerical model has been run for all the tsunamigenic sources (Tables 1 and 2) with a typical simulation time of 5 h and a time step of 2 s. Figure A1 (in Electronic Supplementary Material) shows the time series of slide velocities for some of the studied submarine landslides (Table 1).
Figure 4 shows the computed maximum amplitudes for the submarine landslides SL-1 R6, SL-2 R1 and CSL R1 (Table 1), selected as the extreme cases of each tsunamigenic source.Sources SL-1 and SL-2 are both located at the outer edge of the stable Nile Delta, and they show strong directionality towards the coastline of Israel, were wave amplitudes exceeds 20 m in some areas.For source CSL, most of the tsunami energy propagates northwards.In the three cases, the former Egyptian coastline is much less impacted by tsunami waves.An example of tsunami propagation can be seen in Supplementary material, which shows an animation with some snapshots of the propagation of the tsunami generated by the submarine landslide SL-2 R2.   1. Simulation time is 5 h.
The energy of the tsunami can be evaluated at any time by integration over the whole domain of the potential and kinetic energies of the water column at each grid-cell.Submarine landslides transfer potential and kinetic energy to the water column, competing with energy dissipation by frictional stress at the seafloor (Equations ( 2)-( 4)).Thus, the peak energy usually appears around the time of the maximum velocity for the moving slide.In the scenario of landslide SL-1, with maximum thickness h = 6.0 m, the tsunami peak energy increases with Umax,1 from 5.3 × 10 14 J in R5 up to 4.2 × 10 15 J in R1 (see Table 1).When a maximum thickness h = 20.0 m is adopted, the tsunami peak energy reaches 8.4 × 10 15 J in R6 (the maximum computed wave amplitudes for this extreme event are shown in Figure 4).In the scenario of landslide SL-2, with maximum thickness h = 26.0m, the tsunami peak energy only slightly increases with Umax,1, from 2.6 × 10 15 J in R4 up to 3.7 × 10 15 J in R1 (see Table 1, and Figure 4 for the maximum wave amplitudes computed for R1).Similarly, for landslide CSL the tsunami peak energy increases from 2.0 × 10 15 J (R3) up to 3.9 × 10 15 J (R1).
Figures 5 and 6 show the computed time series of water elevations for all the landslide tsunamis (Table 1) at three of the synthetic tidal gauges TG-3, TG-4 and TG-5 (see Figure 3).TG-4, at the Shi-Hor Lagoon, is the candidate scenario for the drowning of the Egyptian army, as mentioned in the introduction.For all the runs with SL-1 landslide, a weak withdrawal of the sea of few tens cm takes place, and lasting about one hour.After that, water level suddenly rises about 0.5 m, and then continues increasing for 40-50 min, and reaching levels around or below 1 m for Runs R1 to R5, and below 1.5 m for R5 and R6.For runs R1 to R5, differences in the first maximum water height are of the order of 10-20 cm, what can be seen as a model sensitivity test for the Coulombian friction coefficient (and then for Umax,1).These differences are only of few cm for landslides SL-2 and CSL (see Figure 6).Concerning the evolution of water level at TG-4 for the series of SL-2 tsunamis (Figure 6), it is similar to the one commented for SL-1.For CSL landslides, the withdrawal of the sea lasts about two hours, but interfered by the arrival of a weak wave of few tens cm.At TG-5, in the inner shoreline (Figure 3), the tsunami signal arrives latter, and with a similar although smoother pattern than in TG-4; but the maximum water heights are amplified about 0.5 m for SL-1 and 0.25 m for SL-2 and CSL submarine landslides.At TG-3, landslides SL-1 closely follow the same pattern, with a withdrawal of the sea of about 4 m that lasts one hour, followed by a sudden rise of 2-3 m, and then a continuous increase over 50-60 min up to reach water levels of 3-5 m (this last for the 32.7 km 3 landslide).for the submarine landslide SL-1, runs R1 to R7 (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.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 × 10 13 J for F-1and 1.8 × 10 14 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 delimiting the Bardawil Lake, but with weak effects around the Shi-Hor lagoon, as shown in the computed time series of Figure 8.   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 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.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, 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].

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  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.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.

Figure 1 .
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]).

Figure 2 .
Figure 2. Sketch with the geometrical parameters defining the submarine landslide (following the formalism proposed byHarbitz, 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.

Figure 3 .
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 afterCoutellier 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 (Tables1 and 2) are also depicted: yellow boxes for submarine landslides and red lines for faults.

Figure 7 .
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 inTable1.Simulation time is 5 h.

Figure 10 .
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).

Figure A2 .
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 Tables1 and 2. See Figure3for the location of the synthetic gauges.Initial water depths were 36, 33 and 36 m for TG-1, TG-2 and TG-6, respectively.

Table 2 .
Fault parameters used in the simulations.Geographical coordinates correspond to the fault center.Rake is 90° in all cases.