Bitácoras de Mares Numéricos: la Inundación Zancliense del Mediterráneo

183 Boletín de la Real Academia de Córdoba. BRAC, 165 (2016) 183-196 BITÁCORAS DE MARES NUMÉRICOS: LA INUNDACIÓN ZANCLIENSE DEL MEDITERRÁNEO José María Abril Hernández Académico Correspondiente Palabras clave Mar Mediterráneo. Inundación Zancliense. Crisis salina Messiniense. Simulación numérica. Dinámica de fluidos. RESUMEN El trabajo presenta un estudio de la inundación Zancliense del Mediterráneo (hace 5.3 Ma), aplicando técnicas de simulación numérica de la Mecánica de Fluidos. El modelo hidrodinámico resuelve situaciones dinámicas de inundación y calcula las tasas de erosión del lecho marino. La base de datos batimétrica incluye la configuración dinámica del estrecho de Gibraltar. Estudios previos habían concluido que el 90% final de su llenado acontece en unos dos años, con flujos de 90 ·106 m3/s y un ascenso del nivel del agua de hasta 7 m/día. La revisión de los modelos de incisión empleados junto con los resultados de los nuevos experimentos numéricos ponen en cuestión esta visión de la inundación Zancliense, a la vez que plantean nuevas líneas de investigación. Keywords Mediterranean Sea. Zanclean flood. Messinian salinity crisis. Numerical simulation. Fluid dynamics. ABSTRACT This work presents computational fluid dynamics simulations of the Zanclean catastrophic flood of the Mediterranean (5.33 Ma). The two-dimensional non-linear depth-averaged hydrodynamic model includes a flooding algorithm and the estimation of erosion of the seabed. The database for bathymetry includes the dynamic configuration of the Strait of Gibraltar. Previous studies concluded that the final filling of 90% of the Mediterranean volume occurs in two years with water inflows over 90 Sv, and the Mediterranean water level increasing at a rate of 7 m/day. We revisited the bedrock incision models; and new results from numerical experiments put into question the previous view of the Zanclean flooding and outline new research lines. INTRODUCCIÓN “Leguas de cielo y águilas: la Campiña desierta […], el seco mar geológico del Ligustinum”1. Cautivos de la mirada del insigne poeta Mario López, nos trasladamos al tiempo en el que las mareas salobres aún se mecían por estas béticas latitudes al paso de la Luna... Antes de la crisis salina del Messiniense (7.2-5.3 Ma) el valle del Guadalquivir era un brazo de mar que se abría paso hacia el Mediterráneo por el corredor del Guadal- 1 Mario López, “Torreparedones”, en Campo de Córdoba (1990). JOSÉ MARÍA ABRIL HERNÁNDEZ 184 horce2. El estrecho de Gibraltar aún no existía, y el Rif permanecía unido a la cordillera Penibética. Al sur, los corredores del Rif conformaban el segundo canal de comunicación entre el Atlántico y el Mediterráneo. El modelo termomecánico de Duggen et al.3 explica cómo una fractura en la corteza oceánica subducente en la región del estrecho habría originado la elevación tectónica en los márgenes continentales ibérico y africano, cerrando así los pasos marinos del Mioceno. El mar Mediterráneo quedó aislado del resto de los océanos del mundo. Con una tasa de evaporación mayor que el aporte de agua de lluvia y de ríos, el nivel del mar descendió cientos de metros, al tiempo que se precipitaban las sales. Los depósitos de evaporitas del Messiniense, con un volumen total de 106 km3 alcanzan espesores de hasta 1500 y 3500 m en las cuencas occidental y oriental del Mediterráneo, respectivamente4. Esta ingente cantidad de sales precipitadas habría requerido la evaporación de 50 veces el volumen del Mediterráneo, por lo que la evaporación debió coexistir con una etapa de aporte débil o intermitente5 desde el Atlántico. El nivel de equilibrio se alcanzó entre 1500 y 2700 m bajo el nivel actual del mar6, con la cuenca occidental prácticamente seca, a excepción del lago de agua hipersalina que cubría la Llanura Abisal Baleárica (Lago Mare). Se produjo la progresiva continentalización de los márgenes expuestos, y los ríos excavaron profundos cañones. Hace unos 5.3 millones de años, el Mediterráneo volvió a conectarse con el Atlántico a través del naciente estrecho de Gibraltar. El proceso es conocido como la inundación Zancliense del Mediterráneo. Se han sugerido diversas hipótesis para explicar sus causas, que van desde la elevación del nivel de las aguas de los océanos por deglaciación7 a la acción tectónica8; pero la hipótesis más aceptada atribuye su origen a la erosión causada por una corriente de agua que fluyó desde el Atlántico hacia el Mediterráneo9. 2 Martín, J.M., Braga, J.C., Betzler, Ch., 2001. “The Messinian Guadalhorce corridor: the last northern, Atlantic-Mediterranean gateway”. Terra Nova 13 (6), 418–424. 3 Duggen, S., Hoernie, K., van den Bogaard, P., Rüle, L., Morgan, J.Ph., 2003. “Deep roots of the Messinian salinity crisis”. Nature 422, 602–606. 4 Meijer, P.T., Krijgsman, W., 2005. “A quantitative analysis of the desiccation and re-fillingof the Mediterranean during the Messinian Salinity Crisis”. Earth Planet. Sci. Lett. 240, 510–520. 5 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. 6 Blanc, P.L., 2006. “Improved modelling of the Messinian Salinity Crisis and conceptual implications”. Palaeogeogr. Palaeoclimatol. Palaeoecol. 238, 349–372.Meijer y Krijgsman (2005), opus cit. 7 Pérez-Asensio, J.N., Aguirre, J., Jiménez-Moreno, G., Schmiedl, G., Civis, J., 2013. “Glacioeustatic control on the origin and cessation of the Messinian salinity crisis”. Glob. Planet. Chang. 111, 1–8. 8 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. 9 Blanc, P.L., 2002. “The opening of the Plio-Quaternary Gibraltar Strait: assessing the size of a cataclysm”. Geodin. Acta 15, 303–317. Loget, N., Den Driessche, J.V., 2006. “On the origin of the Strait of Gibraltar”. Sediment. Geol. 188–189, 341–356. Bitácoras de mares numéricos: la inundación zancliense del Mediterráneo 185 En la década de 1990, a raíz de los estudios geológicos realizados para el proyecto de túnel de conexión de España y Marruecos, aparecieron evidencias de canales erosivos excavados durante la inundación Zancliense10. García-Castella- nos et al. (2009) 11 apoyaron esta hipótesis a través de los resultados de un modelo de incisión validado en ríos de montaña. Sus resultados sugieren que el 90% del agua se transfirió hacia el Mediterráneo en un corto período de tiempo, que va desde unos meses a dos años. Este modelo, sin ningún nivel de detalle espacial (dimensión cero) y con hipótesis gruesas sobre la hidrodinámica en el canal y los procesos erosivos, aporta una visión básica de la inundación Zancliense. Los cálculos pueden mejorarse sustancialmente con las técnicas de la Mecánica Computacional de Fluidos que aplicamos en este trabajo, aportando una visión más amplia, detallada y crítica de este singular evento. En la sección de Metodología se presentan el modelo físico y las técnicas numéricas que nos permiten recrear el mar en los bits de una computadora. Se fundamentan e ilustran con referencias a diversos trabajos científicos y, en particular, con un amplio conjunto de aplicaciones participadas por el autor en las últimas décadas (nuestras bitácoras de mares numéricos). Los resultados del estudio de la inundación Zancliense del Mediterráneo se presentan y discuten en la Sección tercera. Por una parte, recogen algunos de los aspectos esenciales del trabajo realizado en colaboración con el profesor R. Periáñez, Computational fluid dynamics simulations of the Zanclean catastrophic flood of the Mediterranean (5.33 Ma)12. Se presentan, por otra parte, nuevos experimentos numéricos y nuevos análisis sobre la hidrodinámica y los protocolos para el cálculo de la erosión que nos harán revisar algunas de las ideas aceptadas hasta ahora. METODOLOGÍA Las ecuaciones de Navier-Stokes, que expresan los principios de conservación de la masa y del momento lineal para una partícula fluida, constituyen la piedra angular de la Mecánica de Fluidos. Para su resolución numérica en problemas de dinámica oceánica se aplican técnicas que involucran una discreti - zación espacial del dominio de interés, y también de la línea temporal (el tamaño de ambas está ligado por criterios de estabilidad numérica de las soluciones), sobre las que se promedian las ecuaciones fluido-dinámicas y se aplican hipótesis simplificadoras adaptadas al contexto del problema de interés. En problemas de transporte reactivo sobre escalas temporales de décadas en aguas someras no estratificadas resulta apropiada una descripción bidimensional Loget, N., Den Driessche, J.V., Davi, P., 2005. “How did the salinity crisis end?” Terra Nova 17, 414–419. 10 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. 11 Opus cit. 12 Periáñez, R., Abril, J.M., 2015. “Computational fluid dynamics simulations of the Zanclean catastrophic flood of the Mediterranean (5.33 Ma)”. Palaeogeography, Palaeoclimatology, Palaeoecology 424, 49–60. JOSÉ MARÍA ABRIL HERNÁNDEZ 186 de flujo residual, generado por la función de corriente (potencial escalar) 13 . Para la resolución de problemas que involucran dinámica mareal en aguas someras hemos de recurrir a un modelo bidimensional con ecuaciones integradas en la profundidad e imponer la oscilación mareal en los contornos abiertos del dominio14. El flujo puede estar también forzado por vientos o por diferencias en el nivel medio del mar en los contornos abiertos15. En sistemas marinos con fuerte estratificación, o siempre que se requiera resolver la estructura vertical de las corrientes, se hace necesaria la aplicación de modelos hidrodinámicos puramente tridimensionales16. Para abordar el problema de la inundación Zancliense del Mediterráneo, con una superficie de 2.5·106 km2 y con una escala temporal de varios años, utilizaremos las ecuaciones promediadas en profundidad: 0)()( = ∂ ∂ + ∂ ∂ + ∂ ∂Hv y Hu xt ζ (1)         ∂ ∂ + ∂ ∂ =+Ω− ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ 2 2 2 2 y u x u C H v x g y u v x u u t u A u ρ τ ζ (2)         ∂ ∂ + ∂ ∂ =+Ω+ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ 2 2 2 2 y v x v C H u y g y v v x v u t v A v ρ τ ζ (3) donde u y v son la componentes de la velocidad del agua, promediadas en profundidad, en la dirección de los ejes x (W-E) e y (S-N), h es la profundidad del agua bajo el nivel medio del mar, y ζ es el desplazamiento de la superficie libre sobre dicho nivel medio (medida positive hacia arriba). Así, H=h+ζ es la profundidad total del agua, Ω es el parámetro de Coriolis (Ω= 2w senλ, donde w es la velocidad angular de rotación de la Tierra, y λ es la latitud geográfica), g es la aceleración de la gravedad, ρ es la densidad del agua, y CA es el coeficiente horizontal de viscosidad turbillonaria. τu y τv son los esfuerzos de fricción con el lecho marino, que se escriben en términos de una ley cuadrática: 2222 ;vuvkvuuk fvfu +=+= ρτρτ (4) donde kf es el coeficiente de fricción con el fondo. Estas ecuaciones expresan la conservación de la masa (1) y del momento lineal (2-3), escritas en coordenadas Cartesianas. 13 Abril, J.M., García-León, M., 1993. “A 2D-4Phases Marine Dispersion Model for Non- Conservative Radionuclides. Part I: Conceptual and Computational Model”. J. Environ. Radioactivity 20, 71-88. Abril, J.M., García-León, M., 1993. “A 2D-4Phases Marine Dispersion Model for Non- Conservative Radionuclides. Part II: Two Applications”. J. Environ. Radioactivity 20, 89-115. 14 Periáñez, R., Abril, J.M., 1996. “Modelling the Dispersion of Non-Conservative Radionuclides in Tidal Waters - Part 2: Application to 226Ra Dispersion in an Estuarine System”. J. Environ. Radioactivity 31, 253-272. 15 Abril, J.M., Abdel-Aal, M.M., 2000. Marine Radioactivity Studies in the Suez Canal, Part I: Hydrodynamics and Transit Times. Estuar. Coast. Shelf S. 50, 489–502. 16 Toscano, M., García-Tenorio, R., Abril, J.M., 2005. “Making Predictions on the Evolution of Radioactive Spots in the Ocean. Validation in the Baltic Sea”. Pag. 41-48.En: Recent Advances in Multidisciplinary Applied Physics. Ed. 1. Oxford, Reino Unido. Elsevier Science. ISBN 0-08- 044648-5. Bitácoras de mares numéricos: la inundación zancliense del Mediterráneo 187 Las ecuaciones se resuelven utilizando un esquema explícito de diferencias finitas de segundo orden. En particular, el esquema MSOU (Monotonic Second Order Upstream) se usa para los términos advectivos no lineales en las ecuaciones del momento. Para CA y kf se han usado inicialmente valores de 10 m2/s y 0.0025, respectivamente, de validez probada en gran diversidad de aplicaciones. Como condiciones iniciales se establecen aguas en reposo con los niveles de equilibrio existentes en el Atlántico y en el Mediterráneo antes del comienzo de la inundación. El plano de referencia para ζ se establece en el nivel inicial del Mediterráneo. Como condiciones de contorno, se impide el flujo a través de contornos sólidos de tierra, y se prescribe el nivel de las aguas en el Atlántico. Conforme avanza la inundación, las aguas invaden zonas previamente desecas. Cuando el espesor de la lámina de agua que cubre una determinada celdilla de cálculo supera un valor umbral (típicamente 10 cm), ésta adquiere la etiqueta de celda de agua, aplicando así el esquena numérico de Kampf17. Este modelo ha sido utilizado con éxito para simular mareas y tsunamis generados por diferentes mecanismos (terremotos asociados a fallas geológicas, deslizamientos submarinos, la entrada de flujos piroclásticos en el mar y el colapso de calderas de volcán) en el Atlántico y el Mediterráneo18. Las curvas hipsométricas (superficie libre como función de la profundidad) del Mediterráneo en la actualidad y en el Mioceno tardío son similares, según muestran Meijer y Krijgsman (2005)19. Por esta razón la batimetría actual, que puede conocerse con gran nivel de detalle, será la utilizada en el presente estudio, con las modificaciones pertinentes en la zona del estrecho de Gibraltar. La malla de cálculo que cubre el Mediterráneo se extiende desde 29° N a 46° N y desde 6° W hasta 37° E. La batimetría y elevaciones se han obtenido de la base de datos GEODAS, con una resolución de 4 minutos de arco. Para el estudio detallado de las corrientes y la erosión en la zona Gibraltar-Alborán se ha definido un subdominio de mayor resolución (30” de arco, tomado de la batimetría GEBCO08). 17 Kampf, J., 2009. “Ocean Modelling For Beginners”. Springer-Verlag, Heidelberg. 18 Periáñez, R., Abril, J.M., 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. Abril, J.M., Periáñez, R., Escacena, J.L., 2013. “Modeling tides and tsunami propagation in the former Gulf of Tartessos, as a tool for Archaeological Science”. J. Archaeol. Sci. 40, 4499–4508. Periáñez, R., Abril, J.M., 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. Periáñez, R., Abril, J.M., 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. Abril, J.M., Periáñez, R., 2015. “A Numerical Modelling Study on the Potential Role of Tsunamis in the Biblical Exodus”. J. Mar. Sci. Eng. 3, 745 -771; doi:10.3390/jmse3030745. 19 Opus cit. JOSÉ MARÍA ABRIL HERNÁNDEZ 188 El modelo 0D de García-Castellanos et al. (2009)20 ha sido mejorado y corregido en varios aspectos, como la inclusión del efecto Bernouilli (que reduce el nivel libre del agua en un estrechamiento de la sección del flujo, hasta un centenar de metros en el caso que nos ocupa) y la disipación de energía y caída de nivel por fricción, pero conservando los protocolos de cálculo de la erosión. Este modelo permite estimar la evolución temporal de la sección transversal en el estrecho, e imponerla como condición de batimetría variable en el modelo hidrodinámico 2D antes descrito. La velocidad en el umbral de Gibraltar, vs, se estima a partir del radio hidráulico, Rh, y del gradiente hidráulico, S, por la fórmula de Manning: 2/13/2 1SR n vhs = ;τ b = ρgz s S; Q = zs Wv s ; a bb sk dt dz )( τ = ; w a w QkW = (5) donde n (~0.05) es el coeficiente de rugosidad. García-Castellanos et al. (2009) tomaron S=H/L, con L = 100 km y H = z0-z1 (la diferencia entre los niveles libres del agua en el Atlántico y el Mediterráneo occidental – véase esquena en Fig. 1), con un valor máximo de 1000 m, y Rh ≈zs-z0 , siendo zs la profundidad del canal. Nuestros cálculos distinguen la profundidad de la lámina de agua, zs*, lo que nos permite adaptar consistentemente las estimaciones de τb, Rh, y S. El modelo hidrodinámico corregido sigue siendo una simplificación bastante gruesa, pues el flujo discurre por un canal de sección mojada y pendientes variables, lo que plantea la necesidad de contrastar los cálculos con modelos hidrodinámicos más complejos, como el modelo 2D arriba presentado, capaces de resolver los gradientes hidráulicos a partir de los principios fundamentales de conservación de la masa y del momento. Independientemente del cálculo preciso del caudal que fluye bajo una configuración determinada de la sección en el estrecho y los niveles en el Atlántico y el Mediterráneo, ha de resolverse la tasa de erosión que provoca dicho flujo y el consiguiente incremento en la sección transversal del canal. En este punto mantenemos inicialmente las estimaciones de García-Castellanos et al. (2009), que hacen crecer irreversiblemente la anchura del canal con el caudal (W en Eq, 5) mientras la tasa de incremento de su profundidad (dzs/dt) es una función del esfuerzo cortante (Eq. 5). Para las constantes empíricas kw y aw se adoptan los valores 1.2 y 0.5, respectivamente, y para kb y a 1.30·10-4 m año-1 Pa-a y 1.5, respectivamente. Como condiciones iniciales se toman z1= 2500 m, z2=2700 m, que se corresponde con la situación al final de la crisis salina del Messiniense. Ambas cuencas están conectadas por el umbral de Sicilia, con una profundidad zSC= 430 m. Cuando se supera esta cota, el modelo de García-Castellanos et al (2009) permite que la totalidad del flujo que entra por Gibraltar se inyecte en el Mediterráneo oriental, omitiendo así los detalles de la topografía del estrecho de Sicilia. Para el cálculo del llenado del Mediterráneo usamos las curvas hipsométricas de Meijer y Krijgsman (2005)21. Los cálculos parten de una pequeña incisión con zs= 1 m (el valor inicial de W, dentro de unos márgenes realistas, no afecta al resultado). 20 Opus cit. 21 Opus Cit. Bitácoras de mares numéricos: la inundación zancliense del Mediterráneo 189 RESULTADOS Y DISCUSIÓN la energía de la inundaCión ZanClienSe A partir de las curvas de Meijer y Krijgsman (2005) de superficie libre en función del nivel de llenado, A(z), y conocidos los niveles iniciales, es posible estimar la energía potencial gravitatoria que se liberó en la inundación Zancliense (~9.2 millones de Megatones ): J108.3)( 22 0 max ⋅== ∫ z pdzzAzgE ρ SeCCión, flujo en el eStreCho de gibraltar y Curva de llenado CalCuladaS Con el Modelo 0d En Fig. 2 se muestran los resultados de la aplicación del modelo 0D para el cálculo del caudal y la sección transversal en el estrecho (curva M1). La escala de tiempos está centrada en el instante en el que el nivel de las aguas en la cuenca occidental supera el umbral de Sicilia (t = 26.11 años desde condiciones iniciales). Tras una etapa de bajo caudal, el efecto retroalimentado de erosión, aumento de sección e incremento del caudal, hace que éste crezca hasta alcanzar su valor pico de 180 Sv (1 Sv equivale a 106 m3/s). Cuando el nivel de las aguas en el Mediterráneo se aproxima a la sección mojada del estrecho, se amortigua el caudal de entrada hasta que se equilibran los niveles del Mediterráneo y el Atlántico (éste disminuye unos 9 m durante el proceso). La aplicación de modelos 2D de alta resolución, que se expondrá más abajo, revela que a medida que aumenta la incisión en la zona oriental del estrecho de Gibraltar, las profundidades más someras del umbral de Camarinal, al oeste del estrecho, se convierten en el factor limitante del flujo, estabilizándose éste alrededor de los 80-90 Sv hasta que se amortigua en el estadio final del llenado del Mediterráneo. Este valor umbral ha sido introducido en el modelo 0D, y los resultados aparecen en la curva M2 de la Fig. 2. De esta forma se incrementa el tiempo requerido para el llenado, y la sección transversal alcanza valores algo inferiores al caso de M1. En la actualidad el estrecho de Sicilia tiene una profundidad media de 269 m, pero su configuración pudo ser distinta en el momento de la inundación por los procesos geológicos que continúan moldeando esta zona22, así como por la propia acción del flujo Zancliense. La curva M3 en la Fig. 2 muestra los resultados del cálculo cuando se adopta una profundidad de 250 m para este estrecho. El mayor nivel que alcanza el Mediterráneo occidental limita el flujo en Gibraltar y su poder erosivo. Los tiempos de llenado del 90% final del volumen del Mediterráneo son 456, 580 y 740 días, para M1, M2 y M3, respectivamente. Por economía de espacio no se muestran las curvas de llenado. llenado del Mediterráneo (Modelo 2d Con reSoluCión de 4’ de arCo) Las ecuaciones 1-4 y sus desarrollos asociados constituyen una cadena de cálculo para resolver u, v y ζ en un instante de tiempo dado y en una celda de la malla. Con la resolución espacial de 4’ de arco el dominio comprende algo más 22 Corti, G., Cuffaro, M., Doglioni, C., Innocenti, F., Manetti, P., 2006. “Coexisting geodynamic processes in the Sicily Channel”. In: Y. Dilek & S. Pavlides (Eds), Postcollisional tectonics and magmatism in the Mediterranean region and Asia. Geol. Soc. Am. Spec. Paper, 409: 83-96. JOSÉ MARÍA ABRIL HERNÁNDEZ 190 de 165000 celdillas, y operamos con pasos de tiempo de entre 10 y 20 segundos. Así, para completar la simulación de un año de hidrodinámica es necesario resolver 3.5·1011 cadenas de cálculo. El umbral de Gibraltar se define por dos celdas centradas en 5.5°W. Para cada valor de la profundidad puede encontrarse el caudal cuasi-estacionario con simulaciones de varios días. La baja resolución espacial solo puede proporcionar una aproximación gruesa a la hidrodinámica en el estrecho de Gibraltar, pero el ajuste de las parejas de caudal y sección permite construir una función analítica de respuesta, que usaremos para imponer las condiciones de flujo que se deseen estudiar. Este ejercicio de modelización tiene como objetivo obtener una perspectiva espacial del proceso de llenado del Mediterráneo, como se discute en el trabajo de Periáñez y Abril (2015), al que nos remitimos en este punto. Las condiciones iniciales de la sección y el nivel en el Mediterráneo se toman del modelo 0D, lo que nos permite arrancar la simulación en cualquier instante de tiempo. A título de ejemplo, las figuras 3 y 4 recogen mapas de distribución espacial de niveles y corrientes en dos instantes del llenado. hidrodináMiCa en diStintoS eStadioS de evoluCión del eStreCho de gibraltar (Modelo 2d Con reSoluCión de 30’’ de arCo) El dominio se define entre las longitudes 6.000º W a 0.808º W y las latitudes 35.000º N y 37.608º N. Comprende 195000 celdas de 30” de arco de lado. En la frontera abierta del Mediterráneo se prescribe el nivel del agua como condición de contorno, que se actualiza en función del flujo que se calcula en el estrecho y usando la curva hipsométrica. El estrecho de Gibraltar se define aquí por el conjunto de celdas comprendidas entre 5.56º W (donde se encuentran las menores profundidades y que definirán la sección limitante) y los 5.39º W. Para recrear de una forma simple distintos estadios de evolución del estrecho, las profundidades actuales en la sección limitante se multiplican por un factor constante kS < 1, y para las colindantes se define un perfil de profundidades parabólico que converge a las profundidades actuales en la sección de 5.56º W. Fijado un valor de kS y un nivel inicial en el Mediterráneo occidental, se resuelven las ecuaciones hidrodinámicas hasta alcanzar un flujo estacionario (típicamente en menos de 24 horas). Aunque la tasa de erosión se calcula en cada paso de tiempo, los cambios resultan despreciables en el ámbito temporal de la simulación. El paso de tiempo se ajusta en función de los criterios de estabilidad, siendo el valor mínimo usado de 0.25 s. Con este paso de tiempo la simulación de 24 horas requiere la resolución de 6.7·1010 cadenas de cálculo, y estos ejercicios numéricos se repiten decenas de veces con distintas configuraciones del estrecho y del nivel del mar. A título de ejemplo, en Fig. 5 se muestran las corrientes cuasi-estacionarias calculadas para kS= 0.4 (la sección transversal limitante es de 2.6 km2). El flujo estacionario alcanza los 61 Sv, y las velocidades superan los 30-50 m/s a lo largo del canal de incisión identificado por García-Castellanos et al. (2009). En Fig. 6 comparamos los flujos calculados en función de la sección transversal en 5.56º W (esto es, para distintos valores de kS) usando resoluciones espaciales de 30” y 60” de arco con la fricción nominal kf =0.0025. El nivel del Mediterráneo se fija en -1500 m a fin de que no limite el flujo. Una baja resolución espacial sobreestima el caudal máximo, aunque produce resultados satisfactorios por debajo de los 70 Sv. Cuando la sección crece por encima de los 3.7 km2 las aguas más someras en el umbral de Camarinal (5.74º W) toman el Bitácoras de mares numéricos: la inundación zancliense del Mediterráneo 191 control en la limitación del flujo. Cuando la fricción se incrementa (k f =0.01) el flujo máximo queda por debajo de los 70 Sv. En la misma figura se muestran los caudales calculados con el modelo 0D-M2 en función de la sección del umbral después de las correcciones introducidas en la hidrodinámica23 (aquí el valor límite del flujo es impuesto). Las discrepancias son mayores en el cálculo del esfuerzo cortante τ b , que obviamente depende del coeficiente de fricción k f (Eq, 4). La estimación del modelo 0D (Eq. 5) produce valores mayores para caudales bajos, lo que incide en la rapidez de la incisión. Sobre esta cuestión volveremos para la estimación de las tasas de erosión. La atenuación del flujo al subir el nivel en el Mediterráneo sigue aproximadamente un perfil logarítmico para cada valor de k S , y puede obtenerse en secuencias de simulaciones con desniveles crecientes (no se muestran por economía de espacio). La resolución de 30” tiene también sus limitaciones para valores muy bajos de kS , pues una estrecha corriente de agua no puede resolverse de forma precisa con un número muy pequeño de celdas. Finalmente, el flujo alcanza situaciones supercríticas (número de Froude >1) que obligan a adaptar las propias ecuaciones hidrodinámicas (esta partes del trabajo se encuentran en progreso). No obstante, los resultados obtenidos son suficientes para acometer una revisión crítica del tamaño de la inundación catastrófica del Zancliense. Efectivamente, hemos de notar que la resolución de la hidrodinámica para distintas secciones del estrecho no supone resolver “la ecuación del tiempo” (los lapsos de tiempo involucrados entre una configuración y otra). Para abordar esta cuestión debemos volver sobre los protocolos de cálculo de la erosión. reviSión del Modelo de inCiSión y de la eSCala de la inundaCión ZanClienSe Los procesos implicados en la incisión de lechos rocosos en ríos de montaña, que constituyen el símil más próximo para la inundación Zancliense, son de una gran complejidad, e involucran la abrasión por el flujo de sedimentos que transporta la corriente, la cavitación, fractura y desplomes de ladera, entre otros24. En particular, la cuantificación de la erosión lateral sigue siendo un problema abierto. Las bases teóricas aún son débiles, y dominan los enfoques empíricos con validez restringida a las condiciones experimentales. Resulta conveniente un análisis energético considerando la energía por unidad de volumen requerida para erosionar la roca, εV, que según Lamb et al. (2015)25, para los procesos de abrasión es una función del módulo de Young, Y, y de la resistencia a la tracción del material, σT: Yk TVV / 2 σε = , con kV~106. El orden de magnitud de εV es ~109 Jm-3. Resulta así posible relacionar la potencia por unidad de área empleada en la erosión con la potencia disipada por la corriente de agua: vf dt dz s V τε = Aquí f <1 expresa la fracción de la potencia disipada que se consume en la erosión, y depende de forma compleja de las condiciones del flujo, de la carga de 23 Periáñez y Abril (2015). Opus cit. 24 Lague, D. 2015.The stream power river incision model: evidence, theory and beyond. Earth Surfcade Processes and Landforms. (In press). 25 Lamb, M.P., Finnegan, N.J., Scheingross, J.S., Sklar, L.S. 2015. New insights into the mechanics of fluvial bedrock erosion through flume experiments and theory. Geomorphology 244: 33-55.