An exploratory modelling study on sediment transport during the Zanclean flood of the Mediterranean

An exploratory modelling study onsediment transport duringtheZanclean flood oftheMediterranean R. Periáñez1 · J.M.Abril1· D.Garcia‑Castellanos2· F.Estrada3· G.Ercilla3 Abstract A nearly 400-km-long erosion channel through the Strait of Gibraltar has been interpreted as evidence for a catastrophic refill of the Mediterranean at the end of the Messinian salinity crisis, 5.33 million years ago. This channel extends from the Gulf of Cadiz to the Algerian Basin and implies the excavation of ca. 1000km3 of Miocene sediment from the Alboran Basin and bedrock from the Strait of Gibraltar. The fate of these eroded materials remains unknown. In a first attempt to predict the distribution of those flood deposits, we develop a numerical model to simulate the transport of material eroded from the Strait of Gibraltar. It is a Lagrangian model based upon standard sediment transport equations which is able to simulate suspended and bed load sediment transport. Water circulation during the flood was obtained from a hydrodynamic model of the whole Mediterranean Sea previously developed by the authors and applied to the Zanclean flood. Five particle sizes were considered for suspended load and three for bed load transport. Areas of sediment deposition in the Mediterranean Sea were determined. In the case of suspended load, these are related to hydrodynamic conditions: areas sheltered from the jet of incoming water by local topography and areas where water currents abruptly decrease due to a sudden increase in water depth. In the case of bed load transport, sediments follow water streamlines and deposits are much more localized than in the case of suspended load. Single channel seismic records were also analysed to identify and characterize flood-related deposits in the eastern Alboran Sea. Keywords Numerical model· Suspended load· Bed load· Deposition· Mediterranean Sea· Zanclean flood 1 Introduction The closure of the Guadalhorce and Rifian gateways (Fig.1a), which were the connections between the Atlan-tic Ocean and the Mediterranean Sea before the Messin-ian (7.2–5.3 Ma), limited the water exchange and led to the “Messinian Salinity Crisis” (MSC). During the MSC (5.97–5.33 Ma), the whole Mediterranean basin was at least partially isolated from the world ocean [16, 18, 44, 45], resulting in widespread salt precipitation and a decrease in the Mediterranean sea level of the order of thousand metres according to these authors. Following this extended interpretation, the Mediterranean Sea was later abruptly refilled during the so-called Zanclean flood. Discussions persists regarding the timing, triggering and modality mechanisms of this process (see the review by [44], and further discussion on several open questions in [1]). García-Castellanos etal. [17] reported strong evidence for a deep incision channel along the Gibraltar Strait from boreholes and seismic data generated in the frame of the Africa–Europe tunnel project. The erosion channel has a length of more than 400 km from the Gulf of Cadiz [13] to the Alboran Sea [14]—see Fig.1 for locations of geographic names mentioned in the text— with a varying width (2 to 8 km) and depth (200 to 600 m). García-Cas- tellanos etal. [17] postulated that the observed channel was excavated by the Zanclean flood (thus it is denoted as the Zanclean Channel) and applied a one-dimensional model which indicated that 90% of the water was transferred towards the Mediterranean in a short period, ranging from few months to two years. These results were later confirmed through numerical simulations carried out using a two-dimensional depth-averaged model of the whole Mediterranean Sea [39]. More recently, Abril and Periáñez [1] carried out simulations in which an erosion model was included within the fluid dynamics model, allowing to estimate how the erosion channel was excavated through time. This modelling study showed that interpreting the incised channel in the Camarinal area (connecting with the Zanclean Channel) in a scenario of a catastrophic flooding of the Mediterranean requires achieving peak water flows of at least several tens of Sv.1 Also, the main geological features of the Zanclean Channel, including a sill depth of a few hundred metres at Gibraltar, could be understood from a scenario of catastrophic flooding of the Mediterranean with initial conditions consisting of a wide sill surpassed by a thin water layer. In that work, the modelled scenario which better led to the mentioned facts implied a peak water flow of 70 Sv. This value was achieved when the water level at the Mediterranean was about 170 m below Fig. 1 a Western Mediterranean palaeogeography during the early Messinian [29]. b Map of the computational domain showing geographic names mentioned in the text and present-day (black) and Messinian (red) coastlines; the last from a 2400 m water level drop. The black square corresponds to the zoom area in Fig.3, and the blue one to the area in Fig.13. CamSill, Camarinal Sill; GSill, Gibraltar Sill; GS, Gibraltar Strait; AS, Alboran Sea; AB, Algerian Basin; Bal, Balearic Islands 1 1 Sv =106m3 /s. the Atlantic level, as will be discussed in Sect.2.4. At this stage, the giant jet of water crossing the Strait of Gibraltar produced bottom shear stresses of 1.8 ×104 Pa and incision rates of 1.4 m/day (see their Figs.5, 9 and 10) in this area. The authors estimated that the amount of removed material should have been of the order of 1 km 3 per day. At the same time, the height of the water column in the Alboran Sea was high enough to ensure small bottom shear stresses and negligible erosion rates in comparison with those of the Strait of Gibraltar. Thus, the remaining open question is: where the ca. 10 3 km3 of seafloor eroded by the flood was deposited? Answering this question may lead to an independent validation (or refutation) of the catastrophic flood hypothesis. Sediments were eroded due to the intense currents existing in the Strait during the flood and transported towards the Mediterranean, where they had to be deposited when currents were not strong enough to keep them in movement. Consequently, large deposits of sediments coming from the Strait of Gibraltar should be present somewhere in the Mediterranean Sea. Therefore, the purpose of this work is to investigate, using a sediment transport model, where sediments could have been deposited. Single channel seismic records (320 cubic inch) were analysed to identify and characterize flood-related deposits in the eastern Alboran Sea. Because we use a bathymetry reconstruction from the present-day seabed elevation as a proxy for the Miocene Mediterranean, model results must be interpreted with caution. The aim of the present paper is to show the general relationships between topography and sediment deposition. In any case, numerical models provide a reliable description of the real world as they conform virtual laboratories where different ideas can be tested. In relation to the Zanclean flood, the following conclusions were obtained: 1. The first models for a catastrophic Zanclean flood were independent of the initial conditions [17]. Abril and Periáñez [1] showed that when lateral erosion is described in terms of causal processes instead of a poorly justified empirical relation, initial conditions consisting of a narrow stream of water could not evolve towards catastrophic dimensions. 2. Gibraltar Sill was the only obstacle to water flow considered in the first (zero-dimensional) models for a catastrophic Zanclean flood. A 2D-hydrodynamic model revealed that Camarinal Sill, which is located about 10 km to the west from Gibraltar, exerted a significant limitation to the water flow [39]. 3. A channel excavated in the course of tens of ky by the continuous or pulsed inflow of Atlantic waters (required to explain the total amount of precipitated salts during the MSC), and by the runoff of surface waters (as it can be inferred from the work by [25]), would have existed in the desiccated Alboran Sea previous to the Zanclean flood. The Zanclean flood would have had superposed its erosional fingerprint over this palaeochannel, as demonstrated by the wetting–drying algorithm used in the hydrodynamic model by Periáñez and Abril [39], and Abril and Periáñez [1]. This algorithm was able to encounter the connected topographic-lows defining the path of the Zanclean Channel. 4. Despite some controversial on the continuity of the Zanclean Channel in the Camarinal Sill area, the accepted interpretation is that it consists of an erosive channel incised over the Messinian surface [17]. The western transect of the channel, deeply incised across the former isthmus of Gibraltar, can only be explained by numerical modelling in the scenario of a catastrophic Zanclean flood. This is, high water currents and bottom shear stresses in the Camarinal area appear only for Atlantic water flows of several tens of Sv entering into the Mediterranean [1]. 5. Controversy also exists around the triggering mechanism of the catastrophic flood. Numerical modelling [1] showed that the potential scenarios of initial conditions leading to a catastrophic flood are quite constricted. Thus, it is possible to discard those leading to slow and long-lasting erosive processes. The most likely scenario requires an initial stage consisting of a water layer a few metres thick flowing over the sill with a shear stress over the threshold value for erosion. The achievement of such scenario requires a catastrophic triggering mechanism such as tectonic activity through strike-slip faulting [5, 28]. 6. Once the water flow has been initiated, the only physical stop-condition is the vanishing of the hydraulic gradient between the Atlantic and Mediterranean. Thus, the transgression of Atlantic waters during the MSC should have been maintained in the range of reversibility, possible scenarios being the harmonic coupling between tectonics and erosion in Gibraltar [16] and/or eustatic sea level changes. Similarly, the suggested first step in the reflooding of the Mediterranean (preceding the second catastrophic step) should have involved only moderate water flows. Interpretation of the geophysical findings does not seem to be unambiguous, and thus, more refined modelling approaches are necessary to better illuminate the origin of the Strait of Gibraltar and the details of the Zanclean flood. The present work faces the fate of the eroded material at the former Gibraltar isthmus. It is based upon the “best hydrodynamic and incision scenario” for a catastrophic Zanclean flood, as constrained by previous studies [1]. Model outputs provide new insights which can be used to better support or for discarding the hypothesis of a catastrophic Zanclean flood of the Mediterranean. Thus, the model identifies a series of locations with high probability of containing materials deposited around the peak flow conditions and discards some others. The model is also more concerned in identifying target areas for future field research than for describing the structure of such sedimentary deposits. The model, which is based on standard formulations of sediment transport processes, is described in the Sect.2. First, the hydrodynamic model is briefly described for selfconsistency of the paper. Next, the sediment transport model, its numerical solution and parameter setup are described. Results are presented and discussed in Sect.3, including some sensitivity analysis results and experimental evidences of the presence of Zanclean deposits in the Mediterranean. 2 Model description A sediment transport model requires water depths and currents over the considered domain. These are generally provided by a hydrodynamic model. The hydrodynamic model is the one described in Periáñez and Abril [39], as applied to simulate the Zanclean flood of the Mediterranean. It is a two-dimensional depth-averaged model. The hydrodynamic model, the sediment transport model and the hydrodynamic setup for simulations are described in the following subsections. 2.1 Hydrodynamic model The 2D depth-averaged hydrodynamic equations are (see for instance [23]): where u and v are the depth-averaged water velocities along the x (west–east) and y (south–north) axis, 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 ( Ω=2𝜔sin 𝜆 , where 𝜔 is the Earth rotational angular velocity and 𝜆 is latitude), g is acceleration (1) 𝜕𝜁 𝜕t +𝜕 𝜕x (Hu)+ 𝜕 𝜕y (Hv)= 0; (2) 𝜕 u 𝜕 t+u𝜕u 𝜕x+v𝜕u 𝜕y+g𝜕𝜁 𝜕x−Ωv+𝜏u 𝜌H=A ( 𝜕2u 𝜕x2+𝜕2u 𝜕y2 ); (3) 𝜕 v 𝜕t +u𝜕v 𝜕 x +v𝜕v 𝜕 y +g𝜕𝜁 𝜕 y +Ωu+𝜏v 𝜌 H =A ( 𝜕 2 v 𝜕x 2+𝜕 2 v 𝜕y 2 ), 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. Parameter values were set to k=0.0025 (dimensionless) and A=10 m2 /s [39]. Equations are solved using explicit finite difference schemes [23] with second order accuracy. In particular, the Monotonic Second Order Upstream (MSOU) is used for the advective nonlinear terms in the momentum equations. The numerical scheme described in Kampf [22] was adopted for the description of wetting–drying processes. The solution of these equations provides the values of water depth H and components of the water velocity, u and v, over the model domain. Full details for the model application to the Zanclean flood may be seen in Periáñez and Abril [39]. This hydrodynamic model was also applied to simulate tsunami propagation in different areas [3, 35, 37, 40], as well as for tidal simulations [36, 38] and other massive flood events [2]. 2.2 Sediment transport The model is able to simulate the transport of particles in suspension (suspended load) and particles which are travelling immediately above the seabed (bed load), which occurs for the larger grain sizes. Equations for each transport mode are presented separately. The sediment transport model is formulated on a Lagrangian framework. Thus, the paths of particles are tracked along the simulation in both transport modes. The Lagrangian approach has been adopted to avoid the large numerical diffusion which would arise from the extremely high flow velocities if an Eulerian model were used. The process of mobilization of a large volume of sediment moving downslope as a mass flow driven by gravity, whose characteristics (from debris flow to a turbidity current) would depend upon a number of factors related to sediment nature and flow density, is not included in this work, as explained in the next section. 2.2.1 Suspended load Any substance present in seawater is transported due to water currents (advection) and mixing due to turbulence (diffusion). Advection is directly calculated from the currents provided by the hydrodynamic model. Horizontal and vertical diffusivities, which are also obtained from (4) 𝜏 u=k𝜌u √ u2+v2 ; 𝜏v=k𝜌v√u 2 +v 2 , the hydrodynamic model, are the required parameters to evaluate turbulent mixing. These formulations are explained below. In addition, particle settling and deposition must be considered. The hydrodynamic model is two-dimensional; thus, it does not calculate a vertical water velocity, u and v being depth-averaged as commented before. Nevertheless, both horizontal and vertical movements of particles are calculated. Water velocity can only produce a horizontal displacement of particles. However, turbulent mixing produces a displacement of particles in both horizontal and vertical directions. Thus, each particle is initially located at a position given by coordinates (x,y,z), where z is depth below the water surface. Horizontal coordinates of the particle change because of advection and horizontal diffusion, while the vertical coordinate of the particle may change due to vertical diffusion and particle settling. Sediment particles are released in the Strait of Gibraltar, just downstream the sill and homogeneously distributed over the transversal section of the Strait. The reason is that turbulence at peak flow conditions is so intense in the Strait (Sect.2.4) that particles would be instantaneously mixed over the water column. Then particles are transported by water currents and mixed by turbulence. Particles fall according to a settling velocity which depends on their size and are deposited on the seabed once they reach the bottom and if the bed stress is lower than a critical deposition stress. This critical stress depends on the particle size as well. Local bed stresses are provided by the hydrodynamic model as described in Sect.2.1. It must be noted that the erosion process itself is not modelled: only the paths of particles released in the Strait of Gibraltar are calculated and new particles are not incorporated to the water column from other regions. This was done since we are only interested in the fate of particles eroded from the Strait of Gibraltar. Advective horizontal transport is calculated from the following equation for each particle: where 𝐫=(x,y) is the horizontal position vector of the particle and 𝐪=(u,v) is the current vector at the particle position. The method used for its numerical solution is briefly described in Sect.2.3. Particle settling is evaluated according to the following equation: where ws is the settling velocity for the corresponding particle size (measured positive downwards) and z is the (5) 𝜕 𝐫 𝜕t =𝐪 , (6) dz dt =ws , vertical location of the particle (measured downwards from the local sea surface). When a particle falls on the seabed, it is deposited if the local bed stress is lower than a critical deposition stress, 𝜏cd , above which deposition does not occur. If deposition is not occurring, the particle is reflected back to the water column. Values for all these parameters are defined in Sect.2.5. The general equation describing turbulent mixing in a fluid is written as: where C is the concentration of the considered substance, H is total water depth and Kh and Kv are, respectively, horizontal and vertical diffusion coefficients. This way of expressing the diffusion equation implies that horizontal diffusion coefficients vary in the horizontal plane and the vertical diffusion coefficient is constant in depth (although it changes with location according to water turbulence). Their formulations are given in Sect.2.3. The standard method for solving this diffusion equation in a Lagrangian framework is briefly presented in Sect.2.3. 2.2.2 Bed load A number of equations to describe bed load transport exist in the literature (a brief review may be seen in [8]). However, they are based upon a bed load transport rate not suitable for a Lagrangian description. Consequently, the approach by Bilgili etal. [6] has been adopted, which can be directly used in a Lagrangian framework. In this approach, the critical flow velocity defining when the sediment movement starts is: where h is a characteristic water depth, g is acceleration due to gravity, d50 is the mean sediment diameter and dmax is the maximum one. Instead of using a characteristic depth, this has been replaced by the local water depth, thus h=H(x,y) , since water depths change in more than one order of magnitude over the model domain. Above the critical velocity, particles are assumed to travel at one-sixth of the depth-averaged current [6]. If the current decreases below Vcri the particle stops its movement. 2.3 Numerical solution A two-way linear interpolation method is used to evaluate water velocity, (u,v), at each particle position from (7) 𝜕 (CH) 𝜕t=𝜕 𝜕x ( HKh 𝜕C 𝜕x ) +𝜕 𝜕y ( HKh 𝜕C 𝜕y ) +HKv 𝜕2C 𝜕z2 , (8) V cri =1.4 √ gd50 ln � h 7d50 � dmax d50 � 1∕7 , the four nearest points to the particle where the hydrodynamic model provides values for the current [9]. Equation5, for horizontal advective transport, is expanded for numerical solution as follows: where Δt is time step used to integrate the model. An additional horizontal advective velocity vector (𝜕Kh∕𝜕x,𝜕Kh∕𝜕y) is included to avoid the accumulation of particles in regions of low horizontal diffusivity [42]. As mentioned, Kh and Kv are, respectively, the horizontal and vertical eddy diffusivities, which are deduced from water circulation. In particular, the Smagorinsky’s scheme [10] has been adopted to describe the horizontal diffusivity: where Δx and Δy are the grid cell sizes in the east–west and south–north directions respectively. Both values are 4 minutes of arc in the present application. The approach used by Lane [24] in a Lagrangian sediment transport model has been adopted for the vertical diffusion coefficient: where k (dimensionless) is the bed friction coefficient used in the hydrodynamic model and H is the local water depth. A standard value of k=0.0025 [43] has been used as already commented. Equation(6), for particle settling, is numerically solved as: where z is depth below the water surface, as mentioned before. A stochastic method is used to describe turbulent mixing. A given particle will jump due to turbulence in any direction of space. This movement is split in the horizontal plane and vertical direction since diffusion coefficients are different. Thus, it is considered that the maximum size of the horizontal step given by the particle, Dh , is [19, 34, 41]: in the direction 𝜃=2𝜋RAN , where RAN is a random number between 0 and 1. This equation gives the maximum size of the step. In practice, it is multiplied by RAN to obtain the real size at a given time and for a given particle. (9) x (t+Δt)=x(t)+uΔt+ ΔK h Δx Δt ; (10) y (t+Δt)=y(t)+vΔt+ ΔK h Δy Δt , (11) K h=ΔxΔy √( 𝜕u 𝜕x ) 2 + ( 𝜕v 𝜕y ) 2 +1 2 ( 𝜕u 𝜕y+𝜕v 𝜕x ) 2 , (12) Kv=k|𝐪|H, (13) z(t+Δt)=z(t)+wsΔt, (14) Dh = √ 12K h Δ t Similarly, the maximum size of the vertical step is [19, 34, 41]: which can be given towards the sea surface or bottom with the same probability. Time step used to solve the equations is fixed as Δt=60 s. This small value is required to satisfy stability conditions in the numerical solution [23, 33] due to the large water velocities occurring during the flood (Sect.2.4). Hydrodynamic and transport model are coupled offline. This means that the hydrodynamic model runs in advance. Then water current and diffusion coefficient fields are stored in files which are read by the transport model. An on-line coupling (simultaneous solution of both models) would be computationally too expensive. 2.4 Hydrodynamic conditions The hydrodynamic model provides the horizontal water currents (u,v) and water depths H over the domain, which are required to force the sediment transport model. The computational grid has been obtained from GEODAS database, available on-line, with a resolution of 4 minutes of arc, both in longitude and latitude. It extends from 29◦N to 46◦N and from 6◦W to 37◦E , thus covering the entire Mediterranean. It is worth noting that a higher spatial resolution also requires a smaller time step and thus a computational cost which can hardly be afforded to study the entire Mediterranean. Limiting the study area to the Alboran Sea or to the Western Mediterranean has the problem of providing reliable boundary conditions at the eastern open boundary, which affects the water circulation pattern. The goal of this 2D 4-arc-minutes model is to generate a reliable water circulation for the whole Mediterranean consistent with the water inflow at peak flow conditions predicted by the higher resolution model by Abril and Periáñez [1], which was developed for the Strait of Gibraltar and the western Mediterranean. To simulate the Messinian sea level, the base level of the present-day bathymetry was dropped to −2400 m. This value was selected since the equilibrium level of the isolated Mediterranean was between 1500 and 2700 m below present sea level, according to Blanc [7]. It was used in the hydrodynamic simulations presented in Periáñez and Abril [39]. The Messinian coastline obtained in this way is shown in Fig.1b (red line). It has to be mentioned that in the more recent calculations of Abril and Periáñez [1], which included an erosion model within the fluid dynamic model (covering only the western Mediterranean), the base level of the Mediterranean was dropped to −1500 m. This value was inferred from geophysical studies of the Messinian Erosional Surfaces [14, 47]. Such drop was also (15) Dv=√ 2 KvΔt, large enough to produce incision depths over 300 m (Zanclean channel depth) in the Gibraltar Sill with an energetically consistent erosion coefficient: in Abril and Periáñez [1] a novel approach for erosion in catastrophic floods was developed. The drop in the gravitational-potential energy of the water flow is partially devoted to erosion (through a yield factor). This way the friction coefficient is related with such yield factor and with the energy required to erode (through all the involved processes) a unit sediment volume. The relevant point here is that the same peak flow conditions were obtained with both values of the initial water level of the Mediterranean. Moreover, in the earlier stages, flow through the Strait of Gibraltar does not depend on the base of level of the dried Mediterranean. The reason is that the sill in the Strait acts as a broad crested weir and, consequently, flow is critical here [31, pp. 662–663]. This implies that information cannot travel upstream from the Mediterranean towards Gibraltar. It is worth noting that water-level and water current fields during peak flow conditions are the only required inputs for the sediment transport model used in this work, as will be detailed below. Model results attain peak flow conditions at the Strait of Gibraltar when water level at the Mediterranean was about 170 m below the Atlantic level, according to Abril and Periáñez [1]. This can be clearly seen in Fig.2, where time evolution of water flow through the Strait of Gibraltar, depth of the eroded sill and Mediterranean sea level are presented from the previous calculations. The shaded area indicates the maximum flow conditions. For these conditions, the accurate reconstruction of the Messinian bathymetry is expected to be less influencing. Instead of simulating sediment transport along the whole flood duration, currents obtained during the peak flow at Gibraltar were used. This is the moment when maximum erosion is produced since shear stress is also larger (equation4) and sediments are transported to longer distances. Moreover, it is not computationally feasible to simulate particle transport during the whole filling period. The peak flow is about 70 Sv (Fig.2) and corresponding currents are about 50 m/s in the Strait of Gibraltar [39]. A zoom of water depths and currents at this stage in the most western part of the Mediterranean may be seen in Fig.3. The hydrodynamic model calculates the water/ land boundaries through a wetting/drying algorithm, as explained. Thus, when the hydrodynamic simulation starts, all the Alboran Sea is dry (land). Water starts to flow from Gibraltar and circulates along the deepest narrows. The coloured areas in the top of Fig.3 are indicating the flooded area of the Strait at peak flow conditions and the limit between the colour and white regions indicates the coastline at that time, which is self-calculated by the model. The solid black line is present-day coastline as an indication. The general circulation pattern obtained in the whole Mediterranean Sea during peak flow conditions at Fig. 2 Computed [1] time evolution of water flow, depth of the eroded sill in Gibraltar and Mediterranean Sea level (measured downwards from the Atlantic Ocean level) along the flood process. The shaded box indicates maximum flow conditions 104 104.5 105 105.5 106 106.5 −600 −500 −400 −300 −200 −100 0 Sea level and sill depth (m) Mediterranean Sea level Sill depth 104 104.5 105 105.5 106 106.5 0 20 40 60 80 Time (y) Flow (Sv) Gibraltar can be seen in Fig.4. The horizontal and vertical diffusion coefficients (Eqs.11 and 12) resulting at this moment from the circulation in Fig.4 are presented in Fig.5. These diffusion coefficients are required to solve sediment transport. As an example, the time evolution of the computed bed stresses over the domain may be seen in Electronic Supplementary Material. Bed stress is an essential factor to define the regions where deposition may occur. Fig. 3 Zoom in the western Mediterranean at the considered stage of flooding. Top: water depths (m). Bottom: water current magnitude (m/s) in logarithmic scale. The area covered in these maps is shown in Fig.1 2.5 Model parameters driving sediment transport Sediment particles transported in suspension are released in the Strait of Gibraltar, just downstream the sill and homogeneously distributed over the transversal section of the strait (blue section in Fig.6), as explained before. Five characteristic sizes were simulated according to the scale in Open University Team [32], which is obtained from the Wentworth scale. These sediment classes are given in Table1. Twenty thousands particles are released for each class. Settling velocities for clay and silt are calculated from Stokes’s law. In the case of sands, experimental curves which give the settling velocity versus grain size were Fig. 4 Current field calculated in the Mediterranean at peak flow conditions. Only one of each 16 calculated vectors is drawn for more clarity. The red line indicates the Messinian coastline −5 0 5 10 15 20 25 30 35 30 32 34 36 38 40 42 44 46 20 m/s Longitude Latitude Fig. 5 Horizontal (top) and vertical (bottom) diffusion coefficients ( m2∕s ) resulting from water circulation during peak flow conditions (Fig.4) in logarithmic scale. The presentday coastline is shown 4 Conclusions A numerical model which simulates both suspended and bed load sediment transports during the Zanclean flood of the Mediterranean was developed. The hydrodynamics was obtained from a computational fluid dynamic model previously developed by the authors and applied to the Zanclean flood. The model was solved for the peak flow conditions under which erosion of the greater part of the former Gibraltar divide took place. The sediment transport model was implemented in a Lagrangian framework. Paths of sediment particles eroded from the Strait of Gibraltar during the flood were calculated until these particles are deposited on the seabed. Thus, the regions where Zanclean deposits could be present were determined. Because of the poorly constrained geography of the Mediterranean during the Late Miocene, the use of these results to predict locations with sediment accumulations must be done with caution. Nevertheless, results provide valuable clues with respect to the distance travelled by the sediments and the areas where deposits could be searched, as Zanclean chaotic deposits identified on seismic profiles in the Eastern Alboran Basin indicate, and/or areas which could be discarded in any attempt to find such Zanclean deposits. Results corresponded to a 20-day-long simulation, from the moment when particles are released in the Strait of Gibraltar. It was tested that results do not change for longer simulations, since most particles were deposited a few days after their release. The age of deposited particles was obtained from a “clock” attached to each particle. The post-flood distributions of sediment transported in suspension and bed load were obtained. In the case of suspended load, particles are deposited in regions of low water currents, i.e., low bed stress. These regions are related to zones sheltered from the water jet incoming the Mediterranean by local topography, areas where a sudden increase of water depth occurs and the centre of eddies. Thus, main deposits could be present at both sides of the Strait of Gibraltar (i.e., both sides of the erosion channel), the centre of the Alboran Sea and at both north and south sides of the Alboran Sea, at its connection with the Mediterranean. Of course, particles transported with the main jet fall according to their settling velocity and are finally deposited on the seabed. The distance of these deposits to the Strait of Gibraltar increases as particle size decreases, as could be expected. Thus, very extensive deposits of mainly silts and fine sands appear between the Balearic Islands and Sardinia. Long deposits are also apparent along Algeria coast. Sediment particles transported as bed load follow streamlines. Thus, deposits of Fig. 13 a Airgun seismic profile showing flood deposits (yellow areas) resting on the Zanclean erosive channel (red dashed line). Purple line represents the base of Pliocene. b Bathymetric map (see also Fig.1) showing seismic survey, red line, and patchy distribution of floodrelated deposits; c uninterpreted view of flood-related deposits. Legend: vertical scale in seconds (two way travel time); horizontal scale in metres very coarse sediment eroded from the Strait of Gibraltar and transported as bed load should be more localized in space than finer grain deposits. Sedimentary register on seismic records in the eastern Alboran Sea are consistent with the presence of such deposits that display a patchy distribution along the floor of the Zanclean Channel and on its flanks. Sediment in suspension is transported towards the east to distances reaching some 2000 km in the case of clay and silt and reaching some 1000 km in the case of coarse sand. In contrast, coarser sediment transported as bed load stays closer to the Strait of Gibraltar. Maximum travelled distances are of the order of some 500 km for granule and reduce to some 100 km for cobble. Acknowledgements Work partially supported by the Spanish Project FAUCES CMT2015-65461-C2-R (MINECO/FEDER). Compliance with ethical standards Conflict of interest The authors declare that they have no conflict of interest. References 1. Abril JM, Periáñez R (2016) Revisiting the time scale and size of the Zanclean flood of the Mediterranean (5.33 Ma) from CFD simulations. Mar Geol 382:242–256 2. Abril JM, Periáñez R, O’Connor JE, García-Castellanos D (2018) Computational fluid dynamics simulations of the Late Pleistocene Lake Bonneville Flood. J Hydrol 561:1–15 3. Abril JM, Periáñez R (2017) A modelling study on tsunami propagation in the Red Sea: historical events, potential hazards and spectral analysis. Ocean Eng 134:1–12 4. Bache F, Olivet JL, Gorini C, Rabineau M, Baztan J, Aslanian D, Suc JP (2009) Messinian erosional and salinity crises: view from the provence basin (Gulf of Lions, Western Mediterranean). Earth Planet Sci Lett 286:139–157 5. Balanyá JC, Crespo-Blanc A, Díaz-Azpiroz M, Expósito I, Luján M (2007) Structural trend line pattern and strain partitioning around the Gibraltar arc accretionary wedge: insights as to the mode of orogenic arc building. Tectonics 26:1–19 6. Bilgili A, Swift MR, Lynch DR, Ip JTC (2003) Modelling bed-load transport of coarse sediments in the Great Bay Estuary, New Hampshire. Estuar Coast Shelf Sci 58:937–950 7. Blanc PL (2006) Improved modelling of the Messinian Salinity Crisis and conceptual implications. Palaeogeogr Palaeoclimatol Palaeoecol 238:349–372 8. Camemen B, Larson M (2005) A general formula for noncohesive bed load sediment transport. Estuar Coast Shelf Sci 63:249–260 9. Clarke S (1995) Advective/diffusive processes in the Firth of Forth. Ph.D. thesis, University of Wales, Bangor, UK 10. Cushman-Roisin B, Beckers JM (2011) Introduction to geophysical fluid dynamics. Elsevier, Amsterdam 11. Do Couto D, Estrada F, Gorini C, Ercilla G, d’Acremont E, Ammar A, Jabour H (2018) 1. Alboran Sea. MSC surfaces and clastics. In: Lofi J (ed) Seismic Atlas of the Messinian Salinity Crisis markers in the Mediterranean Sea, vol 2. Mémoires de la Société Géologique de France. Commission for the Geological Map of the World, Paris, p 13 12. Eisma D (1993) Suspended matter in the aquatic environment. Springer, Berlin 13. Esteras M, Izquierdo J, Sandoval NG, Bahmad A (2000) Evolución Morfológica y Estratigráfica Pliocuaternaria del Umbral de Camarinal (Estrecho de Gibraltar) Basada en Sondeos Marinos. Rev Soc Geol España 13:539–550 14. Estrada F, Ercilla G, Gorini Chr, Alonso B, Vázquez JT, García- Castellanos D, Juan C, Maldonado A, Ammar A, Elabbassi M (2011) Impact of pulsed Atlantic water inflow into the Alboran Basin at the time of the Zanclean flooding. Geo-Mar Lett 31:361–376 15. Estrada F, Do Couto D, Ercilla G, Gorini C, Alonso B, d’Acremont E, Juan C (2018) 1. Alboran Sea. MSC surfaces. In: Lofi J (ed) Seismic Atlas of the Messinian Salinity Crisis markers in the Mediterranean Sea, vol 2. Mémoires de la Société Géologique de France. Commission for the Geological Map of the World, Paris, p 13 16. García-Castellanos D, Villaseñor A (2011) Messinian salinity crisis regulated by competing tectonics and erosion at the Gibraltar arc. Nat Lett 480:359–363 17. García-Castellanos D, Estrada F, Jiménez-Munt I, Gorini C, Fernández M, Vergés J, De Vicente R (2009) Catastrophic flood of the Mediterranean after the Messinian salinity crisis. Nature 462:778–781. https ://doi.org/10.1038/natur e0855 5 18. Hsü KJ, Cita MB, Ryan WBF (1973) The origin of the Mediterranean evaporites. In: Ryan WBF, Hsü KJ, Cita MB (eds) Initial Reports of the Deep Sea Drilling Project 13. US Government Printing Office, Washington, DC, pp 1203–1231 19. Hunter JR (1987) The application of Lagrangian particle tracking techniques to modelling of dispersion in the sea. In: Noye J (ed) Numerical modelling. Applications to marine systems. Elsevier, North-Holland, pp 257–269 20. Ji ZG (2008) Hydrodynamics and water quality. Modeling rivers, lakes and estuaries. Wiley, New Jersey 21. Juan C, Ercilla G, Hernández-Molina FJ, Estrada F, Alonso B, Casas D, García M, Farran M, Llave E, Palomino D, Vázquez JT, Medialdea T, Gorini C, D’Acremont E, El Moumni B, Ammar A (2016) Seismic evidence of current-controlled sedimentation in the Alboran Sea during the Pliocene and Quaternary: palaeoceanographic implications. Mar Geol 378:292–311 22. Kampf J (2009) Ocean modelling for beginners. Springer, Heidelberg 23. Kowalik Z, Murty TS (1993) Numerical modelling of ocean dynamics. World Scientific, Singapore 24. Lane A (2005) Development of a Lagrangian sediment model to reproduce the bathymetric evolution of the Mersey Estuary. Ocean Dyn 55:541–548 25. Loget N, Driessche JVD, Davy P (2005) How did the Messinian salinity crisis end? Terra Nova 17:414–419 26. Lofi J, Déverchère J, Gaullier V, Gillet H, Guennoc P, Gorini C, Loncke L, Maillard A, Sage F, Thinon I (2011) Seismic atlas of the Messinian salinity crisis markers in the offshore Mediterranean domain, CCGM, vol 179. Mémoires de la Société Géologique de France. Commission for the Geological Map of the World, Paris, p 72 27. Lofi J (ed) (2018) Seismic Atlas of the Messinian Salinity Crisis markers in the Mediterranean Sea, vol 2. Mémoires de la Société géologique de France. Commission for the Geological Map of the World, Paris, p 71 28. Luján M, Crespo-Blanc A, Comas M (2011) Morphology and structure of the Camarinal Sill from high-resolution bathymetry: evidence of fault zones in the Gibraltar Strait. Geo-Mar Lett 31:163–174 29. Martín JM, Braga JC, Betzler C (2001) The Messinian Guadalhorce corridor: the last northern, Atlantic–Mediterranean gateway. Terra Nova 13:418–424 30. Micallef A, Camerlenghi A, García-Castellanos D, Cunarro-Otero D, Gutsher MA, Barreca D, Spatola D, Facchin D, Geletti R, Krastel S, Gross F, Urlaub M (2018) Evidence of the Zanclean megaflood in the eastern Mediterranean Basin. Sci Rep 8(1):1078 31. Munson BR, Young DF, Okiishi TH (2002) Fundamentals of fluid mechanics. Wiley, New York 32. Open University Team (2005) Waves, tides and shallow water processes. Butterworth-Heinemann, Oxford 33. Periáñez R (2005a) Modelling the dispersion of radionuclides in the marine environment: an introduction. Springer, Berlin 34. Periáñez R, Elliott AJ (2002) A particle tracking method for simulating the dispersion of non conservative radionuclides in coastal waters. J Environ Radioact 58:13–33 35. Periáñez R, Abril JM (2013) Modelling 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 111–112:223–234 36. Periáñez R, Casas-Ruíz M, Bolívar JP (2013) Tidal circulation, sediment and pollutant transport in Cádiz Bay (SW Spain): a modelling study. Ocean Eng 69:60–69 37. Periáñez R, Abril JM (2014a) 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 139:91–102 38. Periáñez R, Abril JM (2014b) A numerical modelling study on oceanographic conditions in the former Gulf of Tartessos (SW Iberia): tides and tsunami propagation. J Mar Syst 139:68–78 39. Periáñez R, Abril JM (2015) Computational fluid dynamics simulations of the Zanclean catastrophic flood of the Mediterranean (5.33 Ma). Palaeogeogr Palaeoclimatol Palaeogeogr 424:49–60 40. Periáñez R, Cortés C (2018) A modelling study on tsunami propagation in the Caspian Sea. Pure Appl Geophys. https :// doi.org/10.1007/s0002 4-018-2057-9 41. Proctor R, Flather RA, Elliott AJ (1994) Modelling tides and surface drift in the Arabian Gulf: application to the Gulf oil spill. Cont Shelf Res 14:531–545 42. Proehl JA, Lynch DR, McGillicuddy DJ, Ledwell JR (2005) Modeling turbulent dispersion on the North Flank of Georges Bank using Lagrangian particle methods. Cont Shelf Res 25:875–900 43. Pugh DT (1987) Tides, surges and mean sea level. Wiley, London 44. Roveri M, Flecker R, Krijgsman W, Lofi J, Lugli S, Manzi V, Sierro FJ, Bertini A, Camerlenghi A, De Lange G, Govers R, Hilgen FJ, Hübscher Chr, Meijer PTh, Stoica M (2014) The Messinian Salinity Crisis: past and future of a great challenge for marine sciences. Mar Geol 352:25–58 45. Ryan WBF (2009) Decoding the Mediterranean salinity crisis. Sedimentology 56:95–136 46. Tattersall GR, Elliott AJ, Lynn NM (2003) Suspended sediment concentrations in the Tamar estuary. Estuar Coast Shelf Sci 57:679–688 47. Urgeles R, Camerlenghi A, García-Castellanos D, De Mol B, Garcés M, Vergés J, Haslam I, Hardman M (2011) New constraints on the Messinian sealevel drawdown from 3D seismic data of the Ebro Margin, western Mediterranean. Basin Res 23:123–145 48. Uusitalo L, Lehikoinen A, Helle I, Myrberg K (2015) An overview of methods to evaluate uncertainty of deterministic models in decision support. Environ Model Softw 63:24–31 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.