Modelización tridimensional del transporte oceánico de 137cs forzado por vientos. Validación en el mar báltico tras el accidente nuclear de Chernobyl

UNIVERSIDAD DE SEVILLA Departamento de F´ısica At´omica, Molecular y Nuclear Modelizaci´on Tridimensional del Transporte Oce´anico de 137Cs Forzado por Vientos. Validaci´on en el Mar B´altico tras el Accidente Nuclear de Chernobyl TESIS DOCTORAL Autor: Manuel Toscano Jim´enez Director: Jos´eMar´ıa Abril Hern´andez Tutor: Manuel Garc´ıa Le´on Universidad de Sevilla, 2012 AGRADECIMIENTOS A Nani, mi mujer, por su apoyo incondicional a lo largo de esta larga y dif´ıcil carrera investigadora. Siempre est´aah´ı, luchadora las veinticuatro horas del d´ıa, positiva e inteligente. Sin su paciencia infinita, esto no hubiera sido posible. A mis padres, Manuel y Concha, personas de bien hechas a s´ı mismas que no conocen otro camino que el de la sencillez y la entrega a sus hijos en cuerpo y alma. De ellos sigo aprendido a disfrutar con alegr´ıa de las cosas peque˜nas y a mirar con serenidad, a´un en los peores momentos. A mis hermanos, por todo lo que he aprendido de ellos y con ellos. La primera mitad de mi carrera investigadora (1998 – 2005) se consolid´o con el reconocimiento de un sexenio de investigaci´on por la Comisi´on Nacional de Evaluaci´on de la Actividad Investigadora (CNEAI), que pude solicitar por mi condici´on de profesor titular de escuela universitaria. La segunda mitad (2005 – 2012) es la que trata de consolidarse con esta memoria de tesis doctoral. AJos´eMar´ıa Abril, que apareci´o como agua de mayo en la segunda mitad de mi carrera investigadora, por su intachable labor de revisi´on de esta memoria de tesis doctoral, y a Manuel Garc´ıa Le´on, por ser mi tutor desde la ´epoca de los cursos de doctorado en la Facultad de F´ısica, de tan tierno recuerdo. A todos mis compa˜neros del Departamento de F´ısica Aplicada de la Escuela T´ecnica Superior de Ingenier´ıa de la Universidad de Sevilla. Con menci´on especial al grupo joven, nacidos en los a˜nos sesenta y setenta, que ingresaron en el departamento en su ´epoca expansiva dirigida por Marcelo Rodr´ıguez Danta, y que son, sin lugar a dudas, una buena muestra de lo mejor que pas´o por las aulas universitarias espa˜nolas durante los a˜nos ochenta y noventa. De ellos es muy f´acil aprender, y son los testigos m´as fidedignos de lo que aparece en estos agradecimientos y de la evoluci´on de mi carrera cient´ıfica y docente. A mis amigos, especialmente a Pepe, Tete, Pablo y Santiago, por su apoyo y respeto permanentes. A mis hijos, Inmaculada y Fernando, por la alegr´ıa que nos dan cada d´ıa. Iniciar una carrera cient´ıfica cuando se est´a trabajando nunca ha sido f´acil, como sabemos los que hemos emprendido este camino cient´ıfico siendo ya profesores en una escuela t´ecnica espa˜nola, donde el trabajo docente es muy acaparador, y donde personas de gran val´ıa intelectual conocidas y apreciadas por nosotros han tenido que abandonar para no desatender su vocaci´on y obligaci´on docentes. Hasta ah´ı, la cosa puede ser llevadera. Lo malo es que si surgen obst´aculos como los que autor de esta memoria ha tenido que superar, el sendero se empina much´ısimo m´as, se embarra y se vuelve, como lleg´o a decir un amigo m´ıo “heroico”. Yo no dir´ıa tanto, pero es probable que a muchas personas que aparec´eis impl´ıcita o expl´ıcitamente en estos agradecimientos, y que me hab´eis ayudado y acompa˜nado estos largos a˜nos, muchas veces desde el m´as respetuoso silencio, os haya dado un pellizquito en el coraz´on como me ha ocurrido a m´ı al releer este p´arrafo. A todos vosotros, muchas gracias. Tesis Do ctoral ´ Indice de contenidos 1.Introducci´on General y Objetivos de la Tesis ....................................... 11 1. Normativa y formato de esta tesis doctoral ........................................... 11 2. Contexto acad´emico de la tesis ....................................................... 12 3. Elecci´on del tema y objetivo general de la tesis ....................................... 14 4. Objetivos espec´ıficos de la tesis ...................................................... 19 2.Hidrodin´amica General de un Fluido Oce´anico ..................................... 23 1. Introducci´on ......................................................................... 23 2. Ecuaci´on de continuidad y variaci´on de la cantidad de movimiento ................... 24 3. Influencia de la viscosidad ........................................................... 31 4. Turbulencia y ecuaciones de Navier-Stokes promediadas .............................. 37 5. Ecuaciones de movimiento del fluido oce´anico ........................................ 52 6. Caso particular: el problema de Ekman simple ....................................... 54 7. Ecuaci´on de advecci´on-difusi´on: transporte de materia ............................... 67 3.Descripci´on General del Mar B´altico y de sus Procesos F´ısicos ................... 70 1. Situaci´on geogr´afica, dimensiones, batimetr´ıa y origen glacial del mar B´altico ......... 70 2. El espectro de corrientes en el mar B´altico. Aplicaci´on de filtros y difusi´on ........... 74 3. Circulaci´on media en el mar B´altico inducida por vientos ............................. 77 4. Descripci´on de los vientos sobre el mar B´altico ....................................... 79 5. Distribuci´on de la salinidad en el mar B´altico ........................................ 80 6. La influencia del hielo en la circulaci´on del mar B´altico .............................. 83 7 Tesis Do ctoral ´ Indice de contenidos 4.La Circulaci´on en el Mar B´altico y su modelizaci´on ................................ 89 1. Introducci´on. Modelos hidrodin´amicos previos del mar B´altico ....................... 89 2. Modelo previo 3D del mar B´altico: modelo de circulaci´on media ...................... 93 3. Modelo previo de corrientes estacionarias del mar B´altico: Modelo de 8 vientos ...... 108 4. Adaptaci´on de los modelos previos a nuestra malla de c´alculo ....................... 135 5. Generalidades de los dos modelos de circulaci´on adaptados a esta tesis .............. 136 5.La Difusi´on en el Mar B´altico y su Modelizaci´on .................................. 153 1. La energ´ıa de las corrientes en el Mar B´altico ....................................... 153 2. Difusi´on diaria: factores influyentes, datos experimentales y estimaciones ............ 156 3. Sobre la relaci´on entre las diferentes escalas de turbulencia .......................... 169 4. La modelizaci´on de la difusi´on ...................................................... 171 5. Estabilidad y difusi´on vertical ...................................................... 174 6. Otras fuentes de difusi´on seg´un resoluci´on espacial o temporal del modelo ........... 182 6.M´etodos matem´aticos para las ecuaciones de transporte ......................... 191 1. Introducci´on ........................................................................ 191 2. M´etodo de diferencias finitas y algoritmos de limitaci´on de flujo en una, dos y tres dimensiones .............................................. 192 3. Introducci´on al m´etodo lagrangiano de Monte Carlo ................................. 205 4. Calibraci´on y aplicaciones del m´etodo lagrangiano de Monte Carlo .................. 207 7.Modelo Lagrangiano de Transporte Conservativo ................................. 226 8 Tesis Do ctoral ´ Indice de contenidos 1. Introducci´on general ................................................................ 226 2. Estad´ıstica de los vientos en el mar B´altico y su modelizaci´on ....................... 227 3. Aplicaci´on primera: la mancha inicial del Golfo de Finlandia ........................ 232 4. Aplicaci´on segunda: el Mar B´altico ................................................. 244 5. Aplicaci´on tercera: aplicaci´on al r´ıo Guadalquivir. Transporte de salinidad........... 254 8.Materia en Suspensi´on y su Modelizaci´on .......................................... 259 1. Aspectos generales de la materia en suspensi´on en el oc´eano ........................ 259 2. Din´amica general de part´ıculas suspendidas en el oc´eano ............................ 262 3. C´alculo de la tensi´on de fondo a partir del modelo hidrodin´amico ................... 264 4. Modelo de transporte de la materia en suspensi´on ................................... 274 9.Modelo del Transporte No Conservativo de Radion´uclidos ....................... 297 1. Introducci´on ........................................................................ 297 2. Fases: disuelta y s´olida en suspensi´on o sedimento .................................. 297 3. Equilibrio y c´alculo de las probabilidades de transici´on de fase ...................... 301 4. Modelo de transporte, resultados y validaci´on ....................................... 304 10.Sumario, Conclusiones y Trabajo Futuro ........................................... 316 1. Sumario y conclusiones ............................................................. 316 2. Trabajo futuro en relaci´on a esta tesis doctoral ...................................... 321 A. Anexos 1. Bibliograf´ıa ......................................................................... 326 9 16 CAP´ ITULO 1 cient´ıfico internacional RV GAUSS que tom´o medidas en 94 puntos del sistema a diferentes profundidades. km km SWEDEN FINLAND BALTIC REPUBLICS 50 100 200 300 400 500 600 800 800 500 30 0 500 700 200 300 200 50 50 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 0 100 200 300 400 500 600 700 800 900 1000 Figura 1.1: Actividad de 137Cs (Bq m–3) en la superficie del Mar B´altico durante Octubre de 1986, interpolada [189] a partir de medidas experimentales contenidas en la literatura [53, 54, 131]. Una cuesti´on destacable en estas medidas fue la influencia de la picnoclina (publica- Introducci´on General y Objetivos de la Tesis 17 ci´on de la revista Nature [88] de 1987) que actu´o como barrera limitando mucho el flujo de radion´uclidos hacia zonas por debajo de la misma, fen´omeno que tiene una mayor importancia en la mitad sur del Mar B´altico. Respecto a la actividad previa en agua y en sedimentos, el dato importante es que la actividad de 137Cs anterior a Chernobyl, seg´un indica la literatura [76, 66], aproximadamente uniforme en el B´altico y no muy variable con la profundidad. El valor de fondo de referencia se sit´ua en torno a 10–20 Bq m–3 lo que nos indica que los niveles de radiactividad en agua se multiplicaron en media por un factor 30 y convirtieron al Mar B´altico en el ecosistema fuera de la Uni´on Sovi´etica m´as afectado por la cat´astrofe de Chernobyl. El accidente de Chernobyl puso de manifiesto la debilidad [157] de los modelos existentes para hacer predicciones f´ısicas en estas situaciones de emergencia. Exist´ıan a mediados de los a˜nos ochenta modelos oce´anicos del Mar B´altico de diversa naturaleza alimentados por la tradici´on de la escuela escandinava de oceanograf´ıa y meteorolog´ıa (V.W. Ekman 1874–1954, C.G. Rossby 1898–1957, V. Bjerkness 1862–1951, H. Sverdrup 1888–1957, M.Knudsen 1871–1949) vinculada desde anta˜no a las grandes empresas exploradoras (F. Nansen 1861–1930, R. Amundsen 1872–1928) de principios del siglo XX, y por otro lado, se dispon´ıa de una gran informaci´on experimental sobre la radiactividad del sistema antes y despu´es del accidente, gracias en gran parte al largo recorrido de la seguridad nuclear en esta zona del mundo iniciado en la ´epocadeSievert,quefuequien dise˜no en los a˜nos cincuenta la pionera red de estaciones de radiaci´on gamma origen de toda la estructura actual y la ley de seguridad nuclear sueca que data del a˜no 1941. Lo parad´ojico es que cuando se produce el accidente de Chernobyl y se detectan las manchas de la figura 1.1 no hab´ıa un modelo matem´atico espec´ıfico que permitiera hacer predicciones sobre la evoluci´on de las mismas. Y m´as parad´ojico a´un es que m´as de una d´ecada despu´es, el autor de esta tesis hace una estancia en la universidad de Uppsala (Suecia) y este problema sigue sin haberse empezado a resolver. Ese es el motivo por el que eleg´ı este tema para iniciar mi carrera investigadora, la elaboraci´on de un primer modelo del transporte de 137Cs en el Mar B´altico, validable con las medidas disponibles en la literatura, que pudiera suponer un apoyo m´as en la seguridad nuclear de escenarios similares en el futuro; y en sistemas oce´anicos con analog´ıas suficientes en la hidrodin´amica, en los m´etodos matem´aticos empleados, y en los fen´omenos de transporte involucrados, que 18 CAP´ ITULO 1 pueden afectar por igual a la materia radiactiva y a la que no lo es, como la contaminaci´on qu´ımica, los nutrientes y las manchas de petr´oleo. La necesidad de aportaciones a este campo floreciente de la f´ısica aplicada ha sido uno de los motores que han inspirado y motivado el trabajo que se presenta en esta tesis doctoral, haciendo trabajar en equipo a dos ramas de la ciencia como la Oceanograf´ıa y la Seguridad Nuclear, de largo recorrido en este escenario, pero no remando, al menos de una forma tan intensa como aqu´ı, en la misma direcci´on, con el accidente de Chernobyl como catalizador. La modelizaci´on del transporte radion´uclidos tiene una primera etapa de modelos de cajas (o compartimentales) como el de Hallstadius [60] para el Mar del Norte y zonas adyacentes, o el de Nielsen [130] para todo el Atl´antico Noreste. A esta familia pertenecen tambi´en el modelo de Scott [164] para los mares ant´articos, el modelo de Sanchez–Cabeza [160] para el mar Mediterr´aneo, y el de Evans [33], espec´ıfico para el sistema que aqu´ınos ocupa, el mar B´altico. Estos modelos de cajas son una primera aproximaci´on al problema radioecolog´ıa en el oc´eano. A este primer grupo de modelos le sucede una familia que son los modelos din´amicos de radion´uclidos, que tienen en cuenta la posibilidad de variaci´on temporal de los flujos del sistema, y se caracterizan por tener una componente hidrodin´amica mayor. Hay que destacar el modelo de Prandle [149] para el Mar del Norte y los modelos de Abril [3, 4, 5, 6] sobre el mar de Irlanda, los dos ´ultimos incluyendo ya incluso el transporte no conservativo de radion´uclidos. La intensidad de modelos asociados al mar del Mar del Norte esta ligada a las cercanas centrales de reprocesamiento nuclear en Sellafield (Reino Unido) y La Hague (Francia). Hay que citar tambi´en el modelo de transporte para el Pac´ıfico Sur [63] ligado aladispersi´on de radion´uclidos provenientes de las pruebas nucleares que Francia efectu´o en su isla colonial de Mururoa desde los a˜nos sesenta hasta mediados de a˜nos noventa. En relaci´on con otros sistemas oce´anicos del planeta, tambi´en se publicaron, entre otros, un modelo sobre el canal de Suez [10, 11] y otro sobre un estuario atl´antico de la zona SO de Espa˜na [140]. Puede decirse que el modelo de Abril [3, 4] sobre el mar de Irlanda fue, tras su precedente de modelos de cajas de Hallstadius [60], el primer modelo radioecol´ogico din´amico sobre un mar de corrientes gobernadas por las mareas. De manera paralela, el modelo aplicado al mar B´altico que se presenta en esta tesis doctoral y en sus publicaciones pre- Introducci´on General y Objetivos de la Tesis 19 vias [187, 189], representar´ıa, tras sus precedentes de modelos de cajas [33, 130], el primer modelo radioecol´ogico din´amico tridimensional sobre un mar de corrientes gobernadas por el viento. 4.OBJETIVOS ESPEC´ IFICOS DE LA TESIS Para llegar al objetivo general marcado en la tesis hay que ir cubriendo una serie de etapas que de una manera casi cronol´ogica se traducen en los cap´ıtulos del 4 al 9. Los cap´ıtulos 2 y 3 son de car´acter b´asico para que un lector cient´ıfico medio, no necesariamente especialista en esta ´area concreta, pueda adentrarse respectivamente en los fundamentos generales de la mec´anica de fluidos de los sistemas oce´anicos y en las particularidades f´ısicas del sistema central bajo estudio, el Mar B´altico, que abrir´an las puertas de los cap´ıtulos del 4 al 9, necesariamente m´as t´ecnicos. El cap´ıtulo10sededicar´aalas conclusiones generales y al trabajo futuro. Los cap´ıtulos 4 y 5 establecer´an el modelo hidrodin´amico sobre el que se sostendr´ael modelo completo y que adaptar´a el movimiento del fluido a nuestra malla de c´alculo, de baja resoluci´on, con celdas de dimensiones horizontales 20 km×20 km. El modelo radioecol´ogico de referencia existente para el Mar B´altico era un modelo compartimental de 7 cajas [33]; por ello, la elecci´on estructural de un modelo con esta resoluci´on es quiz´as el paso m´as natural, siguiendo la l´ınea de otros investigadores, como Abril [3, 4] en sus modelos dispersi´on de radion´uclidos aplicados al Mar de Irlanda. Nuestro modelo, junto a todas sus simulaciones, va a ser ejecutado en una computadora est´andar. Aunque nunca se ha descartado el uso de grandes computadoras, siempre ha habido cierta vocaci´on de usar computadoras normales porque uno de los objetivos de este trabajo es que aplicaciones futuras puedan ser realizadas por tomadores de decisiones en situaciones de emergencia. Hay que decir que los c´alculos finales de la tesis han puesto al l´ımite la memoria RAM de la computadora. Superada esta primera fase m´as conceptual, ser´a necesario, en las siguientes fases de este proyecto, alternar el uso de la peque˜na y la gran computadora, conjugando el c´alculo conceptual con el c´alculo m´as detallado. El modelo hidrodin´amico tendr´a dos partes: el modelo de circulaci´on y el modelo de difusi´on. (a) El modelo de circulaci´on (cap´ıtulo 4) ser´a el de corrientes medias, llamadas tambi´en 20 CAP´ ITULO 1 corrientes de circulaci´on, de advecci´on, de convecci´on, residuales, o simplemente corrientes. Estas corrientes son inducidas fundamentalmente por los vientos, ya que las mareas en el Mar B´altico son pr´acticamente despreciables. Este modelo ser´a desarrollado con dos variantes que son la adaptaci´on a nuestra malla de c´alculo de sendos modelos previos de circulaci´on del Mar B´altico: un modelo 3D [43] de corrientes medias anuales dise˜nado por Instituto Sueco de Meteorolog´ıa e Hidrolog´ıa (SMHI), que ser´a adaptado a nuestra malla num´erica tal cual, y se convertir´aenel modelo m´as usado en esta tesis (cap´ıtulo7,8y9);yunmodelode8vientos[80] para corrientes estacionarias que ser´a extendido simb´olicamente a unas condiciones meteorol´ogicas algo m´as generales as´ı como implementado y validado num´ericamente para un problema concreto 3D (cap´ıtulo 7). Este ´ultimo ser´a de gran utilidad cuando los vientos dominantes del fen´omeno a estudiar sean distintos de los vientos t´ıpicos del modelo anual. (b) El modelo de difusi´on (cap´ıtulo 5) tendr´a en cuenta los fen´omenos de transporte asociados a la desviaci´on del vector velocidad respecto a su valor medio; desviaci´on que tiene una parte asociada al promedio temporal y otra debida al promedio espacial como consecuencia del tama˜no de las celdas num´ericas con las que vamos a trabajar. Una causa importante de la difusi´on horizontal ser´aladesviaci´on del vector de viento respecto del valor medio, que variar´aseg´unsealaescalatemporaldela ventana de filtrado: anual, mensual, semanal, diaria, horaria, etc. Los coeficientes de difusi´on van a ser grandes porque lo son tanto las escalas espaciales como las escalas temporales del problema. La difusi´on vertical estar´a gobernada por el grado de estratificaci´on del sistema. (c) El cap´ıtulo 6 se encargar´adeponersobrelamesalametodolog´ıa num´erica necesaria para mover en el tiempo a las part´ıculas de agua seg´un el ritmo marcado por el modelo hidrodin´amico, cuyo campo de velocidades se toma como conocido seg´un detalles que se explican en la memoria. Se comparar´an y calibrar´an en este cap´ıtulo el m´etodo de diferencias finitas (MDC) y el m´etodo de Monte Carlo (MMC), coment´andose sus ventajas e inconvenientes. (d) El cap´ıtulo 7 supone una validaci´on del modelo de transporte conservativo (la part´ıculas se mueven en disoluci´on con las del agua) en el Mar B´ altico. Se validar´aelmodelo total en dos aplicaciones, la primera, en relaci´on con el Golfo de Finlandia, con Introducci´on General y Objetivos de la Tesis 21 una escala temporal de unas diez semanas, sometida a la circulaci´on del modelo de 8 vientos; y la segunda, en relaci´on al Mar B´altico en su totalidad, con una escala temporal de unos diez meses, sometida a la circulaci´on del modelo de corrientes anuales. Realmente, se trata del mismo problema dividido en dos, ya que las condiciones finales del primer problema son pr´acticamente las condiciones iniciales del segundo. En cada problema, o simulaci´on, dominan unos determinados procesos f´ısicos, de modo que ambas simulaciones se complementan entre s´ı para validar las diferentes partes del modelo en su globalidad, principalmente cuatro: circulaci´on, difusi´on horizontal y movimientos verticales de difusi´on e intercambio de materia en la picnoclina. Se muestra a continuaci´on una tabla resumen de las simulaciones m´as importantes del cap´ıtulo 7, ya que en ´elseponenenpr´actica gran parte de la metodolog´ıa presentada en los cap´ıtulos anteriores dedicados a la circulaci´on (cap´ıtulo 4), difusi´on (cap´ıtulo 5) y m´etodos matem´aticos (cap´ıtulo 6). Tabla 1.1: Procesos f´ısicos dominantes y modelos de circulaci´on aplicados en cada problema del cap´ıtulo 7 Simulaciones del cap´ıtulo 7 Golfo de Finlandia Mar B´altico completo Procesos Circulaci´on, Difusi´on horizontal, dominantes difusi´on vertical intercambio en la picnoclina Modelo de Circulaci´on 8vientos Media anual No todas ellas han sido validadas del golpe sino que, jugando con la informaci´on experimental disponible, se han hecho distintas simulaciones no todas igualmente sensibles a las distintas partes del modelo en su conjunto. Ambas simulaciones se corresponden con dos escalas de tiempo distintas, unas diez semanas y unos diez meses respectivamente i. En la simulaci´on del Golfo de Finlandia (principio del cap´ıtulo 7) se trabaja en una escala temporal m´as corta y se valida de una manera aceptable la metodolog´ıa de la modelo de circulaci´on, con el centroide de la reciente (semanas posterior al accidente de Chernobyl) mancha nuclear b´altica avanzando unos 150 km hacia el SO durante unos dos meses. Tambi´en aqu´ısehaceunbuenajustey calibraci´on de la difusi´on vertical en las capas m´as superficiales que gracias a su acoplamiento con la estructura vertical de la circulaci´on consiguen reproducir la 22 CAP´ ITULO 1 pr´actica detenci´on de la mancha en el mes de agosto de 1986, que a partir de entonces sigue evolucionando gobernada fundamentalmente por la difusi´on. ii. La evoluci´on posterior de esta mancha tambi´en se utiliza en la simulaci´on del mar B´altico completo durante una escala m´as larga (periodo 1986–1987), en la parte central del cap´ıtulo 7, para calibrar otras partes del modelo como la difusi´on vertical en capas m´as profundas y el intercambio en la picnoclina. Tambi´en se har´a una aplicaci´on local de este modelo al transporte de sal en el estuario del R´ıo Guadalquivir. (e) Los radion´uclidos son part´ıculas que en el medio marino se pueden encontrar en la fase disuelta y el fase s´olida tanto en la materia suspendida como en el sedimento. Por ello, hace falta elaborar un modelo independiente del transporte de la materia en suspensi´on, disciplina a la que se consagra el cap´ıtulo 8. (f) Finalmente, en el cap´ıtulo 9, se estudiar´a el transporte de radion´uclidos en sus fases disuelta y s´olida, aprovechando parte de las metodolog´ıas de los cap´ıtulos anteriores, ydejar´aelcap´ıtulo 10 reservado a las conclusiones y al trabajo futuro en relaci´on con esta tesis doctoral. Cap´ıtulo 2 Hidrodin´amica General de un Fluido Oc´eanico 1.INTRODUCCI´ ON En lo que ahora sigue1se va a considerar la respuesta de un fluido sometido a distintas fuerzas externas e internas, lo que nos llevar´aalaobtenci´on de las ecuaciones b´asicas que describen la din´amica del oc´eano. Desarrollos similares y complementarios a ´este pueden encontrarse en multitud de textos [14, 21, 174]. Como se sabe, la Mec´anica de Fluidos aplicada a la Oceanograf´ıa es fundamentalmente mec´anica newtoniana ampliada con la influencia de los complejos fen´omenos turbulentos. Esta disciplina establece su base en las leyes de variaci´on de la masa, de la cantidad de movimiento, del momento angular y de la energ´ıa. Las fuerzas dominantes en la din´amica del oc´eano son: gravedad, fricci´on y de Coriolis. Recordemos algunas de sus propiedades m´as notables: (a) Gravedad:presi´on, mareas, diferencias de densidad. El peso del agua en el oc´eano produce presi´on. Los cambios en la gravedad debido a los movimientos del sol y la luna relativos a la tierra producen mareas, corrientes 1El desarrollo de este cap´ıtulo introductorio se hace partiendo de un nivel algo inferior a lo habitual en una tesis doctoral, para que puedan seguirlo tambi´en alumnos de ´ultimos cursos universitarios con inter´es en este campo para su proyecto fin de carrera, especialmentelosalumnosdelaE.T.S.deIngenier´ıa de la Universidad de Sevilla, que ya han iniciado una l´ıneadetrabajoenestaramabajolasupervisi´on del doctorando. 23 24 CAP´ ITULO 2 mareales y difusi´on mareal en el interior del oc´eano. El empuje act´ua sobre cada elemento del fluido en funci´on de su diferencia de densidad con el medio circundante. Si por ejemplo empieza a soplar aire fr´ıo sobre la superficie del mar, esto har´aque elaguadelacapasuperiorsevuelvam´as densa que la de la capa inmediatamente inferior, tendiendo la primera a caer y a ocupar un lugar inferior m´as estable. (b) Fricci´on: viento, viscosidad, fondo oce´anico. Es interesante distinguir entre fricci´on interna y externa. •Fricci´on interna: la viscosidad. Su papel es fundamental por ser la responsable de la transformaci´on de la energ´ıa en el seno del fluido. •Fricci´on externa: el viento y el fondo del mar. En el fluido destacan dos capas l´ımite en sendas fronteras del sistema con el aire y el fondo oce´anico. Es notable el papel de las fuerzas externas por su influencia mec´anica en el movimiento global del sistema. En este grupo sobresalen: la fuerza del viento, y la fuerza en el fondo del mar. Adem´as de producir olas, el viento cambia la distribuci´on del campo de presiones y aumenta su influencia al resto del oc´eano. (c) Coriolis. La fuerza de Coriolis es una fuerza de inercia debida a la rotaci´on de la tierra y su influencia se extiende tanto a la atm´osfera como al oc´eano. A gran escala, es la m´axima responsable de la direcci´on y el sentido en que se curvan las l´ıneas de corriente en el sistema mar-aire. 2.ECUACI´ ON DE CONTINUIDAD Y VARIACI ´ ON DE LA CANTIDAD DE MOVIMIENTO Para poder aplicar las leyes de conservaci´on hay que tener en cuenta que un elevado n´umero de las magnitudes escalares y vectoriales van viajando en el fluido sometidas a cambios a lo largo del tiempo, por lo que ser´a necesario recordar el concepto de Derivada Total, D/Dt.Consid´erese una magnitud escalar qque es transportada por un punto material en una dimensi´on x, y coloquemos dos puertas de control separadas una cantidad δx en las que qalcanza los valores {qin,q out}invirtiendo un tiempo δt, quedando el balance como sigue: Hidrodin´amica General de un Fluido Oc´eanico 25 qout =qin +∂q ∂tδt +∂q ∂xδx. (2.1) En la figura 2.1 se hace una representaci´on esquem´atica del balance, conjuntamente con un sistema de ejes cartesianos para la posici´on {x, y, z}ylavelocidad{u, v, w}de las part´ıculas del fluido. Figura 2.1: Esquema ilustrativo de la derivada total en un fluido. El cambio sufrido por qen ese segmento de control vendr´a dado por: Dq Dt = lim δt→0 qout −qin δt =∂q ∂t +∂q ∂x δx δt .(2.2) Si se tiene en cuenta que la velocidad del fluido en la direcci´on xpuede calcularse como u=δx/δt, se llega a Dq Dt =∂q ∂t +u∂q ∂x ,(2.3) y generalizando a las tres dimensiones Dq Dt =∂q ∂t +u∂q ∂x +v∂q ∂y +w∂q ∂z ,(2.4) e introduciendo el operador nabla D Dt =∂ ∂t +v · ∇(2.5) se obtiene una expresi´on bastante m´as compacta. 32 CAP´ ITULO 2 ∂u ∂t +u∂u ∂x +v∂u ∂y +w∂u ∂z =−1 ρ ∂p ∂x +2Ωvsen ϕ+ffric x,(2.31) donde destaca ffric x, fuerza de fricci´on por unidad de masa. Veamos su forma en el caso en que ffric xsea debido a la influencia de la viscosidad, ffric x=fvis x. Retomemos la caja elemental de fluido de dimensiones {δx,δy,δz}para considerar las tensiones (figura 2.4) en la direcci´on de movimiento x. Figura 2.4: Esquema de tensiones en un volumen elemental del fluido. Ayud´andonos del elemento cartesiano de volumen de la figura anterior tendremos las siguientes tensiones en la direcci´on x: •La tensi´on Txz,eslacomponentexaplicada en una cara del tipo z=cte.Habr´a que tener en cuenta las caras superior e inferior de la caja elemental. •La tensi´on Txy,eslacomponentexaplicada en una cara del tipo y=cte.Habr´a que tener en cuenta las caras anterior e posterior de la caja elemental. •La tensi´on Txx,eslacomponentexaplicada en una cara del tipo x=cte.Habr´a que tener en cuenta las caras izquierda y derecha de la caja elemental. Esta tensi´on, a diferencia de las otras dos que son tangenciales, es normal a la cara sobre la que act´ua y est´a relacionada con la presi´on del modo Txx =−p. Hidrodin´amica General de un Fluido Oc´eanico 33 De esta manera, podremos generalizar y obtener las tensiones en la direcci´on yque ser´an del tipo {Tyx,T yy,T yz}, mientras las paralelas a Oz tomar´an la forma {Tzx,T zy,T zz}. Es por ello que para completar el problema necesitamos 3×3 = 9 tensiones representadas por el tensor Tij ,i=1,2,3,j=1,2,3,(2.32) entre las que se destacan por un lado las componentes normales , (Tii =−p)representantes de la presi´on py por otro, las tensiones tangenciales -o de corte- (Tij ,i=j) situadas en los elementos no diagonales del tensor. Si extendemos a todo el espacio tendremos un campo tensorial de tensiones Tij(x, y, z). Volviendo a la fuerza total en direcci´on xsobre el elemento fluido se tiene que: δFx=Txx(x=δx)δSx+Txy(y=δy)δSy+Txz(z=δz)δSz −[Txx(x=0)δSx+Txy(y=0)δSy+Txz(z=0)δSz].(2.33) Las tensiones representadas en la figura son las sufridas por la celda problema debido a celdas situadas en posiciones superiores teniendo en cuenta el sistema coordenado Oxyz. Una celda situada en una posici´on superior sufrir´a debido a la inmediatamente inferior y seg´un el tercer principio de la Din´amica, una tensi´on opuesta a la dibujada en la figura; es por ello que para tener en cuenta la influencia de las tres celdas inferiores sobre la celda problema hayamos tenido en cuenta el sentido correspondiente y de ah´ı el signo menos que aparece en la f´ormula anterior. Teniendo en cuenta los elementos de superficie y de volumen δSx=δyδz , δSy=δxδz , δSz=δxδy , δV =δxδyδz , (2.34) δSx=δV/δx , δSy=δV/δy , δSz=δV/δz , (2.35) los cambios en las tensiones se expresan como sigue: Txx(x=δx)−Txx(x=0)≡δTxx ,T xy(y=δy)−Txy(y=0)≡δTxy ,(2.36) Txz(z=δz)−Txz(z=0)≡δTxz .(2.37) 34 CAP´ ITULO 2 Sustituyendo estos cambios: δFx=δTxxδSx+δTxyδSy+δTxzδSz,(2.38) y por otro lado δTxx =∂Txx ∂x δx , δTxy =∂Txy ∂y δy , δTxz =∂Txz ∂z δz , (2.39) as´ı tendremos que δFxqueda como sigue: δFx=∂Txx ∂x +∂Txy ∂y +∂Txz ∂z δV . (2.40) Introduciendo la fuerza total por unidad de masa ftot x≡δFx/δm,ladensidadρ= δm/δV , e igualando: ftot x=1 ρ∂Txx ∂x +∂Txy ∂y +∂Txz ∂z .(2.41) Para generalizar y calcular las tensiones en las otras dos direcciones y,z,basta permutar xcon yyxcon zrespectivamente, quedando: ftot y=1 ρ∂Tyy ∂y +∂Tyx ∂x +∂Tyz ∂z ,(2.42) ftot z=1 ρ∂Tzz ∂z +∂Tzy ∂y +∂Tzx ∂x .(2.43) En desarrollos posteriores utilizaremos a veces la notaci´on matricial para las componentes cartesianas de algunos campos vectoriales notables: r =[x1,x 2,x 3],v=[v1,v 2,v 3], ftot =ftot 1,ftot 2,ftot 3,{Tij}.(2.44) De este modo queda m´as compacta ftot x: ftot x≡ftot 1=1 ρ 3  j=1 ∂T1j ∂xj ,(2.45) Hidrodin´amica General de un Fluido Oc´eanico 35 y generalizando para ftot i,i=1,2,3: ftot i=1 ρ 3  j=1 ∂Tij ∂xj∀i=1,2,3.(2.46) Para empezar, veamos cu´anto valen la componentes de la fuerza por unidad de masa fν x,fν y,fν zcuando la fricci´on es realizada por influencia de la viscosidad molecular. Conocida es la ley de viscosidad para fluidos newtonianos –como es nuestro caso– para flujo en la direcci´on xy contenido en el plano xz: Txz =ρν ∂u ∂z ,(2.47) donde la viscosidad molecular ρν se define, para fluidos newtonianos como el cociente entre la tensi´on tangencial Txz y el perfil de velocidades ∂u/∂z.Elcoeficienteνes la viscosidad cinem´atica que toma valores de t´ıpicos de ν=10 –6 m2s–1 para el agua a temperatura ambiente. Se generaliza la ley anterior a las tres dimensiones, y se aprovecha la notaci´on matricial, Tij =ρν ∂vi ∂xj∀i=j, (2.48) expresi´on que sustituimos en ftot i, teniendo en cuenta los elementos diagonales de la matriz relacionados con la presi´on, Tii =−p, y llegando a: ftot i=1 ρ 3  j=1 ∂ ∂xj Tij =1 ρ 3  j=1 ∂ ∂xjρν ∂vi ∂xj   ≡fν i −1 ρ ∂p ∂xi   ≡fp i ,i=1,2,3 (2.49) donde tenemos separados los t´erminos de la viscosidad y de la presi´on. Teniendo en cuenta la incompresibilidad del sistema fν i= 3  j=1 ∂ ∂xjν∂vi ∂xj∀i=1,2,3,(2.50) e introduciendo el operador nabla 36 CAP´ ITULO 2 fν i= ∇·ν ∇vi= ∇·ν ∇vi∀i=1,2,3,(2.51) basta sumar las tres componentes para construir  fν:  fν= ∇·ν ∇v. (2.52) Ahora se podr´ıan particularizar las componentes dispuestas a ser aplicadas a los problemas espec´ıficos: fν x=∂ ∂x ν∂u ∂x+∂ ∂y ν∂u ∂y+∂ ∂z ν∂u ∂z,(2.53) fν y=∂ ∂x ν∂v ∂x+∂ ∂y ν∂v ∂y+∂ ∂z ν∂v ∂z,(2.54) fν z=∂ ∂x ν∂w ∂x +∂ ∂y ν∂w ∂y +∂ ∂z ν∂w ∂z .(2.55) Si sustitu´ımos fν xen la ecuaci´on de la cantidad de movimiento: ∂u ∂t +u∂u ∂x +v∂u ∂y +w∂u ∂z =−1 ρ ∂p ∂x +2Ωvsen ϕ+fν x,(2.56) llegamos a la componente xdel movimiento del fluido: ∂u ∂t +u∂u ∂x +v∂u ∂y +w∂u ∂z =−1 ρ ∂p ∂x+ 2Ωvsen ϕ+∂ ∂x ν∂u ∂x+∂ ∂y ν∂u ∂y+∂ ∂z ν∂u ∂z.(2.57) An´alogamente, para las otras componentes {y,z}: ∂v ∂t +u∂v ∂x +v∂v ∂y +w∂v ∂z =−1 ρ ∂p ∂y −2Ωusen ϕ+∂ ∂x ν∂v ∂x+∂ ∂y ν∂v ∂y+∂ ∂z ν∂v ∂z,(2.58) ∂w ∂t +u∂w ∂x +v∂w ∂y +w∂w ∂z =−1 ρ ∂p ∂z Hidrodin´amica General de un Fluido Oc´eanico 37 +2Ωucos ϕ−g+∂ ∂x ν∂w ∂x +∂ ∂y ν∂w ∂y +∂ ∂z ν∂w ∂z .(2.59) Con la notaci´on vectorial queda mucho m´as compacta la expresi´on de la variaci´on temporal del campo de velocidades: ∂v ∂t +v · ∇v =g −1 ρ ∇p−2 Ω×v + ∇·ν ∇v, (2.60) que son nuestra primera versi´on de las ecuaciones de Navier–Stokes, correspondientes a fluidos newtonianos incompresibles. 4.TURBULENCIA Y ECUACIONES DE NAVIER-STOKES PROMEDIADAS 4.1. Introducci´on. La turbulencia es una de las m´as importantes propiedades de los sistemas fluidos, debida a la complejidad natural de un medio continuo cuyas ecuaciones de gobierno para la cantidad de movimiento contienen t´erminos no lineales (de la forma u∂u/∂x ,v∂u/∂y ,w∂u/∂z ,etc)quesonlosm´aximos responsables de la alta dificultad matem´atica que alberga su soluci´on, muy sensible a inestabilidades ligadas a problemas mec´anicos de esta naturaleza, que requieren, como veremos m´as adelante, un tratamiento especial. La importancia de esos t´erminos no lineales viene dada por n´umero el adimensional de Reynolds Re, que se define como el cociente entre el t´ermino no lineal y el t´ermino viscoso. Re =T´ermino no lineal T´ermino viscoso =|(v · ∇)v | | ∇·(ν ∇)v |.(2.61) Si tuvi´eramos por ejemplo un flujo unidireccional (paralelo a Ox) cuyas variaciones se dieran en una direcci´on gen´erica Oy, se tendr´ıa que: Re =T´ermino no lineal T´ermino viscoso =|(v · ∇)v | | ∇·(ν ∇)v |=|u∂u ∂y | |ν∂2u ∂y2|.(2.62) Es muy importante, en todos los problemas fluidomec´anicos, hacer estimaciones del n´umero Re con objeto de poder hacer simplificaciones en las ecuaciones generales de Navier–Stokes. En muchos casos es suficiente el orden de magnitud de Re para conseguir 38 CAP´ ITULO 2 el objetivo deseado. En el problema, relativamente simple, que se presenta ahora como de ejemplo, se cumple que: u≈U|∂u ∂y |≈U L|∂2u ∂y2|≈U L2.(2.63) donde Ues una velocidad caracter´ıstica del campo v,yLes la longitud caracter´ıstica de las variaciones del campo en la direcci´on Oy , que se puede estimar en general sin demasiada dificultad. Sustituyendo de este modo en Re: Re =T´ermino no lineal T´ermino viscoso =|u∂u ∂y | |ν∂2u ∂y2|=UU L νU L2 =UL ν.(2.64) Este n´umero Re adimensional tan notable fue introducido por el cient´ıfico Osborne Reynolds (1842–1912) en su hist´orico experimento [155] cuyo objetivo fundamental era acercarse emp´ıricamente al concepto de la turbulencia. Estudi´o el movimiento de agua fluyendo a lo largo de un tubo y utiliz´o una tinta como trazador. Hizo m´ultiples medidas barriendo desde las velocidades peque˜nas hasta las grandes. Como es sabido, para velocidades peque˜nas el flujo es claramente laminar, mientras que para las grandes es turbulento. La transici´on ocurre cuando Re =UL/ν ≈2000, donde Ues la velocidad promedio y Les el di´ametro del tubo. Una propiedad a destacar del n´umero Re es que dos fluidos con la misma geometr´ıa y el mismo n´umero Re tienen el mismo patr´on de flujo. Por ello fueron notables los trabajos de Richardson [157] sobre el movimiento de un fluido con velocidad t´ıpica Ualrededor de un cilindro circular de di´ametro D: es interesante cuando muestra la comparaci´on entre dos configuraciones tan distintas aparentemente como {D=10cm,U=1 cm s –1} y{D=1 cm,U=10cms –1}, que presentan id´entico patr´on ya que en ambos casos se tiene el mismo Re = 1000. Si hacemos una primera estimaci´on para un sistema oce´anico como el que vamos a estudiar en esta tesis doctoral (Mar B´altico) suponiendo que se cumple la condici´on Re ≈UL/ν, que la velocidad U≈0.1 m s–1, que la longitud t´ıpica de las variaciones L≈ 10 km, y que a temperatura ambiente el agua tiene ν≈10–6 m2s–1,resultaqueRe = 1010Re (cr´ıtico) ∼2000, por lo que nuestro sistema en particular es muy turbulento. Hidrodin´amica General de un Fluido Oc´eanico 39 Esta conclusi´on puede extenderse de manera general a los oc´eanos de nuestro planeta. No quiere decir esto que la viscosidad molecular sea siempre despreciable; ser´a importante en determinadas zonas del oc´eano y en algunos problemas particulares. 4.2. Tensiones de Turbulencia y Ecuaciones de Navier-Stokes Promediadas Una de las grandes dificultades a la hora de resolver las ecuaciones de Navier-Stokes surge debido al car´acter no lineal de este sistema. Para avanzar en este problema un camino tradicional consiste en tomar valores medios en estas ecuaciones previa descomposici´on de las magnitudes f´ısicas en un t´ermino promedio temporal m´as otro de naturaleza turbulenta; as´ı fue como trabaj´o O. Reynolds, el primero en acercarse matem´aticamente al problema de la turbulencia. De este modo, la descomposici´on quedar´ıa [14] del siguiente modo para el vector velocidad y sus componentes: v = V+vu=U+uv=V+vw=W+w,(2.65) donde las velocidades promedio se definen como sigue:  V(t)=v(t)=1 Tt+T/2 t−T/2v(τ)dτ , (2.66) U(t)=u(t)=1 Tt+T/2 t−T/2u(τ)dτ V (t)=v(t)=1 Tt+T/2 t−T/2v(τ)dτ , (2.67) W(t)=w(t)=1 Tt+T/2 t−T/2w(τ)dτ . (2.68) Es muy razonable asumir que en las condiciones oce´anicas con las que vamos a trabajar: p=pρ=ρg=g. (2.69) Dado que la funci´on promedio es una funci´on lineal, se cumplir´a la siguiente propiedad de linealidad para el campo de velocidades: 40 CAP´ ITULO 2 L(v)=L(v).(2.70) donde L(v) es una funci´on lineal en la velocidad. Debe quedar claro que las derivadas temporales importantes que aparecer´an en las ecuaciones de la cantidad de movimiento son operadores que trabajan en una banda de tiempo (Δt) mucho mayor que el tiempo t´ıpico del promedio T: ∂ ∂t()≡lim Δt→0 Δ|r= cte Δt() , V≡lim Δt→0 Δr Δtcon ΔtT.(2.71) Con estas herramientas, podemos pasar a promediar las ecuaciones de Navier-Stokes para la cantidad de movimiento: ∂v ∂t+v · ∇v=g−1 ρ ∇p−2 Ω×v+ ∇·ν ∇v,(2.72) sustituyendo y aplicando la propiedad anterior de linealidad, ∂v ∂t +v · ∇v=g −1 ρ ∇p−2 Ω×v+ ∇·ν ∇v,(2.73) e introduciendo los valores medios, ∂ V ∂t +v · ∇v=g −1 ρ ∇p−2 Ω× V+ ∇·ν ∇ V. (2.74) Falta por desarrollar el t´ermino no lineal: v · ∇v. Para ello, es conveniente demostrar previamente que: v · ∇v = ∇·[vv],(2.75) donde [vv] es la matriz [vivj] de las componentes cartesianas del campo de velocidades, tal que v = 3  i=1 viei= 3  j=1 vjej.(2.76) Hidrodin´amica General de un Fluido Oc´eanico 41 Recordando definiciones:  ∇≡∂ ∂x1 ∂ ∂x2 ∂ ∂x3,[vv]≡⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ v1v1v1v2v1v3 v2v1v2v2v2v3 v3v1v3v2v3v3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .(2.77) Desarrollemos el producto de nabla por la matriz de las velocidades  ∇·[vv]:  ∇·[vv]=∂ ∂x1 ∂ ∂x2 ∂ ∂x3·⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ v1v1v1v2v1v3 v2v1v2v2v2v3 v3v1v3v2v3v3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ =3  i=1 ∂(viv1) ∂xi 3  i=1 ∂(viv2) ∂xi 3  i=1 ∂(viv3) ∂xi = 3  j=1 3  i=1 ∂(vivj) ∂xi ej=⎛ ⎝ 3  j=1 vjej⎞ ⎠ = ∇·v   3  i=1 ∂vi ∂xi + 3  j=1 3  i=1 vi ∂vj ∂xi ej.(2.78) En condiciones normales oce´anicas de fluido incompresible (se cumplen las condiciones de Boussinesq) la ecuaci´on de continuidad se expresa simplemente con  ∇·v =0,por lo se tiene:  ∇·[vv]= 3  j=1 3  i=1 vi ∂vj ∂xi ej=3  i=1 vi ∂ ∂xi3  j=1 vjej=v · ∇v, (2.79) de esta manera, v · ∇v = ∇·[vv],(2.80) que era lo que se quer´ıa demostrar. Hay que resaltar tambi´en que esta propiedad es v´alida para todo campo de velocidades incompresible. Volvamos al desarrollo del t´ermino no lineal v · ∇v. Vamos a demostrar en las siguientes l´ıneas que: v · ∇v= ∇·[vv]= ∇· V V+ ∇· v v.(2.81) 48 CAP´ ITULO 2 Componente i: 2 3  j=1 3  k=1 ΩkVjεijk donde Ω1=0 Ω 2=ΩcosϕΩ3=Ωsenϕ. (2.116) Del t´ermino convectivo de la aceleraci´on ten´ıamos:  ∇· V V= 3  j=1 3  i=1 ∂ ∂xi (ViVj)ej.(2.117) Permutando icon j, y trabajando con nabla, aparecen las componentes ide la parte convectiva de la aceleraci´on: Componentes ide ∇· V V: 3  j=1 ∇j(VjVi).(2.118) Volviendo a escribir las ecuaciones N–S en t´erminos conceptuales tenemos que la aceleraci´on de un punto del sistema a es igual a la resultante  fde las fuerzas actuantes por unidad de masa: a = f⇒aloc +acon   aceleraciones = fg+ fp+ fcor + fν+ fR   fuerzas / unid. masa .(2.119) Recapitulando las tres componentes i={1,2,3}de las ecuaciones N−Spara la cantidad de movimiento: ∂Vi ∂t  aloc i + 3  j=1 ∇j(VjVi)   acon i =−gδi3  fg i −1 ρ∇ip   fp i + 2 3  j=1 3  k=1 ΩkVjεijk   fcor i + 3  j=1 ∇jTν ij ρ   fν i + 3  j=1 ∇jTR ij ρ   fR i ,(2.120) siendo Tν ij el tensor de viscosidad newtoniana y TR ij el tensor de turbulencia proveniente de un E.V.M. (Eddy Viscosity Model): Tν ij =ρν ∂Vi ∂xj =ρν ∇jViTR ij ≡−ρv jv i=E.V.M. =ρAj ∂Vi ∂xj .(2.121) Hidrodin´amica General de un Fluido Oc´eanico 49 4.3. Casos particulares relevantes La discusi´on que viene a continuaci´on trata sobre las aproximaciones que se pueden hacer a los t´erminos din´amicos fν i,f R iy simplificaciones que esto conlleva sobre las ecuaciones de movimiento del sistema. En general tenemos las siguientes expresiones para las fuerzas fν i,f R i, admitiendo la aproximaci´on newtoniana para la primera y el modelo de viscosidad turbillonaria para la segunda: fν i= 3  j=1 ∇jTν ij ρ=Modelo de Newton = 3  j=1 ∇j(ν∇jVi),(2.122) fR i= 3  j=1 ∇jTR ij ρ=Modelo Eddy Viscosity = 3  j=1 ∇j(Aj∇jVi).(2.123) expresiones matem´aticamente equivalentes donde los coeficientes {ν;Aj}representan respectivamente a la viscosidad molecular y a la viscosidad turbulenta. (a) Caso (1). Viscosidad molecular constante, ν=cte: fν i= 3  j=1 ∇j(ν∇jVi)=ν 3  j=1 ∇j(∇jVi)=ν⎛ ⎝ 3  j=1 ∇j∇j⎞ ⎠Vi=ν⎛ ⎝ 3  j=1 ∇2 j⎞ ⎠Vi= =ν ∇2Vi=def. laplaciano : ∇2≡Δ=νΔVi⇒ fν=νΔ V, (2.124) la expresi´on de las fuerzas ser´ıa: """"""fν i=νΔVi,f R i= 3  j=1 ∇j(Aj∇jVi)"""""".(2.125) (b) Caso (2).- Viscosidad molecular constante, ν=cte, y viscosidad turbulenta constante, Aj=ctej(j=1,2,3) : fR i= 3  j=1 ∇j(Aj∇jVi)= 3  j=1 Aj∇j(∇jVi)= 3  j=1 Aj∇2 jVi,(2.126) 50 CAP´ ITULO 2 para compactar la expresi´on se define el vector viscosidad turbulenta  Ay el vector laplaciano que ser´a notado por  Δ:  A≡[A1A2A3], Δ≡∇2 1∇2 2∇2 3,(2.127) quedando como sigue fR i: fR i= A· ΔVi⇒ fR= A· Δ V, (2.128) siendo las expresiones de las fuerzas: """" fν=νΔ V, fR= A· Δ V""""(2.129) (c) CASO (3).- Viscosidad molecular constante, ν=cte, viscosidad turbulenta constante, Aj=ctej(j=1,2,3) , y turbulencia horizontal is´otropa (A1≃A2≡Ah;A3≡ Av):  fR= A· Δ V=(A1Δ1+A2Δ2+A3Δ3) V=(Ah(Δ1+Δ 2)+AvΔ3) V= =def. laplaciano horizontal y vertical: Δ1+Δ 2≡ΔhΔ3≡Δv= =(AhΔh+AvΔv) V, (2.130) siendo las expresiones de las fuerzas: """ fν=νΔ V, fR=(AhΔh+AvΔv) V"""(2.131) que expresada en componentes cartesianas: fν i=νΔVifR i=(AhΔh+AvΔv)Vi.(2.132) Utilizando una notaci´on m´as oceanogr´afica tendremos que las fuerzas de fricci´on se expresar´an de la manera siguiente: •Viscosidad newtoniana.  fν=νΔ V⇒fν x=νΔU=ν∂2U ∂x2+∂2U ∂y2+∂2U ∂z2; (2.133) fν y=νΔV=ν∂2V ∂x2+∂2V ∂y2+∂2V ∂z2(2.134) fν z=νΔW=ν∂W ∂2x2+∂2W ∂y2+∂2W ∂z2.(2.135) Hidrodin´amica General de un Fluido Oc´eanico 51 •Viscosidad turbillonaria.  fR=(AhΔh+AvΔv) V⇒ fR x=(AhΔh+AvΔv)U=Ah∂2U ∂x2+∂2U ∂y2+Av ∂2U ∂z2; (2.136) fR y=(AhΔh+AvΔv)V=Ah∂2V ∂x2+∂2V ∂y2+Av ∂2V ∂z2; (2.137) fR z=(AhΔh+AvΔv)W=Ah∂2W ∂x2+∂W ∂2y2+Av ∂2W ∂z2.(2.138) Sustituyendo en las componentes ide las ecuaciones de N−S: ∂U ∂t + 3  j=1 ∇j(VjU)=−1 ρ ∂p ∂x +2(Ω 3V−Ω2W)+fν x+fR x,(2.139) ∂V ∂t + 3  j=1 ∇j(VjV)=−1 ρ ∂p ∂y −2Ω3U+fν y+fR y,(2.140) ∂W ∂t + 3  j=1 ∇j(VjW)=−g−1 ρ ∂p ∂z +2Ω 2U+fν z+fR z.(2.141) En la inmensa mayor´ıa de los problemas oceanogr´aficos se desprecia la viscosidad molecular frente a la turbulenta (νAi;i=1,2,3). Tambi´en es bien sabido que las dimensiones verticales del oc´eano son mucho menores que las horizontales, propiedad que condiciona tambi´en a las velocidades en las direcciones correspondientes (W{U, V }) y que nos permite despreciar uno de los t´erminos (2Ω2W)delaprimera ecuaci´on. La componente de la fuerza de Coriolis en la direcci´on vertical es mucho menor que el peso (2Ω2Ug) para las velocidades t´ıpicas del sistema, por lo que no ser´a tenida en cuenta. Introduciendo el factor de Coriolis (f≡2Ω3=2Ωsenϕ) y teniendo en cuenta todo lo anterior se van a expresar las ecuaciones N−Spara la cantidad de movimiento en notaci´on oceanogr´afica (U, V, W): ∂U ∂t +∂(UU) ∂x +∂(VU) ∂y +∂(WU) ∂z =−1 ρ ∂p ∂x+ +fV +Ah∂2U ∂x2+∂2U ∂y2+Av ∂2U ∂z2,(2.142) 52 CAP´ ITULO 2 ∂V ∂t +∂(UV) ∂x +∂(VV) ∂y +∂(WV) ∂z =−1 ρ ∂p ∂y+ −fU +Ah∂2V ∂x2+∂2V ∂y2+Av ∂2V ∂z2,(2.143) ∂W ∂t +∂(UW) ∂x +∂(VW) ∂y +∂(WW) ∂z =−g−1 ρ ∂p ∂z+ +Ah∂2W ∂x2+∂2W ∂y2+Av ∂2W ∂z2.(2.144) Al tratarse de un fluido incompresible tenemos el problema mec´anico y el t´ermico desacoplados. 5.ECUACIONES DE MOVIMIENTO DEL FLUIDO OCE´ ANICO Supongamos en principio un medio de densidad ligeramente variable, siempre cumpliendo las hip´otesis de Boussinesq; de modo que la presi´on pueda aproximarse por el modelo hidrost´atico: ρ=ρo+ρ→p=pa+gζ zρdη =pa+gρo(ζ−z)+gζ zρdη (2.145) donde ζes la elevaci´on del mar por encima del nivel medio local z=0,yHes la profundidad de modo que el punto inferior de la columna de l´ıquido est´aenlacota m´ınima z=−H. Comolaanomal´ıa de densidad ρes mucho menor que el valor medio ρo(ρρo), podemos aproximar la inversa de la densidad de esta manera: 1 ρ=1 ρo+ρ≃1 ρo (2.146) Por otro lado, podemos intentar simplificar la siguiente integral: gζ zρdη =g0 zρdη +gζ 0ρdη (2.147) Como en los casos que vamos a estudiar ζes del orden de pocos cent´ımetros y zdel orden de pocos o muchos metros, ser´a razonable suponer: gζ zρdη ≃g0 zρdη (2.148) Hidrodin´amica General de un Fluido Oc´eanico 53 Las componentes horizontales de la velocidad, la ecuaci´on de continuidad, y el campo de presiones estar´a gobernado por las siguientes ecuaciones: ∂U ∂t −fV +∂(UU) ∂x +∂(VU) ∂y +∂(WU) ∂z =−1 ρ ∂p ∂x +Av ∂2U ∂z2+Ah∂2U ∂x2+∂2U ∂y2(2.149) ∂V ∂t +fU+∂(UV) ∂x +∂(VV) ∂y +∂(WV) ∂z =−1 ρ ∂p ∂y +Av ∂2V ∂z2+Ah∂2V ∂x2+∂2V ∂y2(2.150)  ∇· V=0 ⇒∂U ∂x +∂V ∂y +∂W ∂z = 0 (2.151) p=pa+gρo(ζ−z)+g0 zρdη (2.152) Tengamos en cuenta que: ρ=ρo+ρ⇒∂ρ ∂xi =∂ρo ∂xi +∂ρ ∂xi⇒ ⇒∂ρ ∂xi =∂ρ ∂xi⇒(2.153) ⇒∂p ∂xi =∂pa ∂xi +gρo ∂ζ ∂xi +g0 z ∂ρ ∂xi dη ∀i=1,2,3 (2.154) Sustituyendo (2.154) en las ecuaciones (2.149) y (2.150) e introduciendo para abreviar el operador laplaciano horizontal Δh: ∂U ∂t −fV +∂(UU) ∂x +∂(VU) ∂y +∂(WU) ∂z = =−1 ρo∂pa ∂x +gρo ∂ζ ∂x +g0 z ∂ρ ∂xdη+Av ∂2U ∂z2+AhΔhU(2.155) ∂V ∂t +fU +∂(UV) ∂x +∂(VV) ∂y +∂(WV) ∂z = =−1 ρo∂pa ∂y +gρo ∂ζ ∂y +g0 z ∂ρ ∂ydη+Av ∂2V ∂z2+AhΔhV(2.156) 54 CAP´ ITULO 2 ∂U ∂x +∂V ∂y +∂W ∂z = 0 (2.157) De esta manera, queda planteado un sistema de tres ecuaciones (2.155), (2.156), (2.157) con tres inc´ognitas {U, V, W},dadalahip´otesis de que se conocieran las funciones restantes {ρ, pa,ζ,A z,A h}involucradas en el el sistema. 6. CASO PARTICULAR: EL PROBLEMA DE EKMAN SIMPLE Antes de entrar en soluciones num´ericas del problema, parece interesante mostrar algunas soluciones particulares como son las ecuaciones libre de contornos laterales y de fondo. Es el caso particular m´as simple de todos los problemas de circulaci´on inducida por vientos y sirve para darnos una primera estimaci´on de funciones tan importantes como: •El di´ametro φpara las ´orbitas inerciales. Esta magnitud se define para un flujo gobernado ´unicamente por la gravedad y la fuerza de Coriolis. A este primer caso podemos denominarlo flujo inercial simple. •El perfil vertical de velocidades (h´elice de Ekman) para flujo horizontal, homog´eneo y estacionario. A este segundo caso podemos denominarlo flujo Ekman simple. Aqu´ı se estudiar´atambi´en espesor de la capa Ekman, zona donde se absorbe la inmensa mayor´ıa de la energ´ıa externa que proporciona el viento. Esta estimaci´on va a ayudarnos a entender e interpretar mejor las soluciones estacionarias de los campos {ψ(x, y),ζ(x, y),vM(x, y),vs(x, y)}correspondientes respectivamente a la funci´on de corrientes, la funci´on elevaci´on, la velocidad media en la columna de l´ıquido, y la velocidad en superficie, a las que llega Jankowski [80] en su modelo integrado de corrientes bidimensionales inducidas por los 8 vientos en las direcciones principales, con una intensidad de 10 m s–1 a una altura de diez metros sobre la superficie del mar B´altico. Estas soluciones pueden extenderse a las tres dimensiones [80], quedando este campo estacionario de corrientes evaluado en funci´on de la profundidad, y preparado para futuras aplicaciones en esta tesis doctoral, donde se tendr´a en cuenta tambi´en la capa Ekman de fondo (o bent´onica), de una naturaleza matem´atica similar a la de superficie, pero inducida l´ogicamente por la tensi´on en la base de la columna de l´ıquido. Esta soluci´on de Ekman m´as general, que veremos en el cap´ıtulo 4, aportar´at´erminos dependientes del contorno y de la batimetr´ıa en general, destacando la contribuci´on de las corrientes geostr´oficas. Hidrodin´amica General de un Fluido Oc´eanico 55 6.1. Flujo inercial simple. Estamos ante el caso m´as simple del transporte Ekman. Para estudiar las corrientes cerca de la superficie del oc´eano es interesante empezar por aqu´ı, considerando este caso tan particular [146] en que el sistema tras ser sometido a un fuerte viento durante unas horas deja de recibir este impulso y sigue movi´endose por inercia, sometido ´unicamente a la fuerza de la gravedad y a la de Coriolis; despreciando la influencia de la fricci´on. Las ecuaciones de movimiento, teniendo en cuenta (2.18) y (2.19), quedar´ıan como sigue: du dt =−1 ρ ∂p ∂x +2Ωvsin ϕ(2.158) dv dt =−1 ρ ∂p ∂y −2Ωusin ϕ(2.159) dw dt =−1 ρ ∂p ∂z +2Ωucos ϕ−g(2.160) donde: Ω es la velocidad angular de rotaci´on de la Tierra, con Ω = 2π/Tsid = 7.292 10–5s–1 ;Tsid = 86164 s es el periodo sid´ereo ; y f=2Ωsinϕel factor de Coriolis dependiente de la latitud ϕ. Buscando soluciones simples a estas ecuaciones podemos analizar lo que ocurre si despreciamos el gradiente horizontal de presiones frente al termino de movimiento y al de la fuerza de Coriolis: ∂p ∂x =∂p ∂y = 0 (2.161) Si no se tienen en cuenta la influencia de los contornos puede decirse que el flujo ser´a horizontal, quedando como sigue la ecuaciones de movimiento: du dt =2Ωvsin ϕ=fv (2.162) dv dt =−2Ωvsin ϕ=−fu (2.163) Se trata de un sencillo sistema de ecuaciones diferenciales lineales acopladas que puede ser resuelto mediante t´ecnicas tradicionales. Derivando la segunda ecuaci´on y sustituyendo en la primera: 56 CAP´ ITULO 2 du dt =−1 f d2v dt2=fv (2.164) Reordenado t´erminos, puede reconocerse la ecuaci´on diferencial del oscilador arm´onico: d2v dt2+f2v= 0 (2.165) Este sistema se integra de una manera a´un m´as simple utilizando notaci´on vectorial. Usando el vector velocidad v =[u, v], el sistema de ecuaciones se expresa de la forma: ˙ v =−fez∧v =f[−v,u] (2.166) Salta a la vista que el vector aceleraci´on es perpendicular al vector velocidad: a =˙ v =−fez∧v ⇒a ⊥v (2.167) Esto quiere decir que la aceleraci´on tangencial es nula y s´olo queda la normal: a =aNaT=0 ⇒|v|=cte ⇒|aN|=|a|=f|v|=CTE (2.168) Como el movimiento es uniforme (|v|=cte) y con aceleraci´on normal constante |aN|=CTE, tenemos que la soluci´on es el movimiento circular uniforme cuyo radio R ser´a: R=|v|2 aN =|v| f=constante (2.169) El llamado periodo inercial viene dado por Tine =2πR/|v|y define el llamado d´ıa pendular en alusi´on al tiempo empleado por su plano de oscilaci´on al realizar un ciclo completo. Teniendo en cuenta que f=2ΩsinϕyΩ=2π/Tsid, resultan las siguientes expresiones para el periodo y el di´ametro inercial {φ, Tine}: φ=2|v|/f =|v|Tsid 2πsin ϕTine =Tsid 2sinϕ(2.170) Hidrodin´amica General de un Fluido Oc´eanico 57 A continuaci´on se muestra la gr´afica de un interesante experimento [196] sobre corrientes inerciales medidas a una profundidad de 15 m en el Pac´ıfico Norte, en octubre de 1987, durante un periodo de 25 d´ıas. Las posiciones fueron observadas 10–12 veces al d´ıa mediante el sistema Argos de sat´elites meteorol´ogicos de ´orbita polar de la NOAA. Las corrientes mayores se generaron el tercer d´ıa debido a una tormenta que tuvo lugar en esazonadeloc´eano. Figura 2.5: Corrientes inerciales medidas en el Pac´ıfico Norte [196]. Cada bucle corresponde con una oscilaci´on inercial, que ocurre aproximadamente cada 17 horas. Para ver la variaci´on del periodo y del di´ametro inerciales en funci´on de la latitud se muestra una tabla ilustrativa para una velocidad t´ıpica de deriva |v|=20cms –1: TABLA 2.1: ´ Orbitas inerciales para |v|=20cms –1: Latitud ϕ(o)Periodo inercial Tine (h) Di´ametro φ(km) 10 68.93 15.8 40 18.62 4.3 60 13.82 3.2 90 11.97 2.7 6.2. Flujo Ekman simple. Las soluciones para este flujo particular se obtienen imponiendo la condici´on de car´acter estacionario en las ecuaciones generales (2.155) y (2.156), a la que hay que 64 CAP´ ITULO 2 Otra forma de llamar a DEva a ser la siguiente: DE≡Dπ. Integrando en la columna de agua, se ha demostrado que el flujo medio en la capa Ekman viaja en la direcci´on que forma 45ocon la velocidad superficial; ortogonalmente por tanto a la velocidad del viento U10. Vamos a definir esa profundidad como Dπ/4.Sevef´acilmente que Dπ/4=Dπ/4≡DE/4. Tambi´en resulta interesante definir la profundidad D⊥:v⊥ Vo.Secalcula f´acilmente que D⊥=DE/2. Las tres profundidades definidas anteriormente {DE,D π/4,D ⊥}est´an relacionadas con la direcci´on del campo de velocidades del fluido. Es por ello que a continuaci´on se va a definir una profundidad relacionada con el m´odulo de las velocidades. Se define profundidad cr´ıtica Dcpara el decaimiento exponencial 1/e: z=−Dc→v(z)=vc=Vo/e (2.205) Es f´acil obtener que: Dc=DE/π . A continuaci´on se muestra una tabla particularizada para el Mar B´altico con valores t´ıpicos de estas profundidades y de la velocidad superficial Voen funci´on de la velocidad del viento U10. TABLA 2.3: Algunas caracter´ıcticas de la capa Ekman en el Mar B´altico U10(m/s); Viento Vo(cm/s); Oc´eano DE(m) Dπ/4(m) Dc(m) D⊥(m) 57.5 27 7 9 14 10 15 45 11 15 23 20 30 100 25 32 50 Una caracter´ıstica de la capa Ekman es que en su interior la fuerza de fricci´on y la de Coriolis alcanzan valores del mismo orden. Es por ello que se define el cociente entre ambas fuerzas, magnitud adimensional llamada n´umero de Ekman E: E=Fuerza de Fricci´on Fuerza de Coriolis =Az∂2U ∂z2 fU (2.206) Si introducimos una velocidad t´ıpica Uy una profundidad representativa den laqueseproduzcanlasvariacionesdeU, podemos aproximar el numerador como sigue: Hidrodin´amica General de un Fluido Oc´eanico 65 Az ∂2U ∂z2≈Az U d2⇒E=AzU/d2 fU =Az fd2(2.207) Si evaluamos Eal final de la capa Ekman: d=DE=#2π2Az f⇒E=Az f f 2π2Az =1 2π2∼ =1 20 (2.208) Como cab´ıa esperar, a esa profundidad el n´umero de Ekman ya es muy peque˜no. Se comprueba que Edecrece con la profundidad, tendencia que se mantiene hasta que aparece la fricci´on debida al fondo del oc´eano, es decir hasta que nos aproximamos a la otra capa l´ımite. Por debajo de la capa Ekman, y lejos del fondo marino, el n´umero Ees muy peque˜no (E1), predominando la influencia de la fuerza de Coriolis. 6.3. La capa Ekman atmosf´erica. La capa Ekman de la que acabamos de hablar es una zona superficial en el oc´eano sensible doblemente a la influencia del viento y de la fuerza de Coriolis. Aunque es la m´as importante en lo que a los problemas de transporte oc´eanico se refiere, es obligatorio mencionar otra notable capa Ekman que acompa˜na a la anterior en los problemas tierraoc´eano: Capa Ekman atmosf´erica, que es la capa l´ımite de aire en contacto con el oc´eano. La soluci´on matem´atica a este problema indica que el flujo neto de viento cerca de la interfase con el oc´eano, caracterizado por  U10,formaun´angulo de 45ohacia la izquierda -en el hemisferio nortecon respecto al flujo atmosf´erico de altura. Seg´un la aproximaci´on geostr´ofica el flujo atmosf´erico de altura se mueve perpendicularmente a los gradientes de presi´on siguiendo una direcci´on pr´acticamente isob´arica. Si tenemos en cuenta que las corrientes oce´anicas superficiales forman un ´angulo de 45oaladerechade  U10, de lo cual resulta que tenemos una curiosa compensaci´on de la suma de ambos efectos ambos efectos cuyo resultado es que lascorrientessuperficiales oce´anicas inducidas por los vientos siguen pr´acticamente las l´ıneas isob´aricas. Para ilustrar esta importante propiedad se va a mostrar a continuaci´on una gr´afica representativa de un experimento [120] realizado en el Oc´eano Pac´ıfico. Se representan 66 CAP´ ITULO 2 las isobaras promediadas para el mes de abril de 1978, y las trayectorias seguidas por las boyas de prueba a la deriva. Se observa que las trayectorias transcurren sobre caminos pr´acticamente isob´aricos, a excepci´on de la zona de la corriente de Kuroshio, cerca del archipi´elago japon´es, donde no se re´unen las condiciones necesarias para esta clase de transporte. Figura 2.7: Trayectorias de boyas a la deriva, pr´oximas a las l´ıneas isob´aricas [120]. 6.4. Algunas revisiones, limitaciones y validaciones del transporte Ekman. Son innumerables los trabajos que se han realizado con posterioridad a estos pioneros estudios en Oceanograf´ıa. Las correcciones m´as importantes que cabe destacar son: (a) Correcciones por cercan´ıa a la costa, debido a que la soluci´on de Ekman est´a libre de contornos. (b) Intervalo transitorio. El transporte Ekman cl´asico se preocupa del estado estacionario, lo cual requiere que los vientos soplen durante periodos superiores a un d´ıa pendular. De todas maneras, el propio Ekman calcul´otambi´en una soluci´on dependiente del tiempo, del mismo tipo a la que obtuvo a˜nos m´as tarde [61]. (c) Correcciones por estabilidad vertical. Es relativamente simple suponer que Azdepende ´unicamente de U2 10. Es mucho m´as realista suponer que el coeficiente de viscosidad sea dependiente tambi´en de la profundidad. Hidrodin´amica General de un Fluido Oc´eanico 67 (d) Correcciones por la variaci´on de la densidad. Este cuesti´on tambi´en afecta al problema de la estabilidad. Las cuestiones sobre la estabilidad y la difusi´on en general, ser´an ampliadas en el cap´ıtulo 5 de este trabajo. Una gran variedad de cuidadosos experimentos y medidas de corrientes se ha realizado a lo largo del siglo XX [200, 24, 153] confirmando la bondad de la Teor´ıa de Ekman y el papel fundamental de las corrientes inerciales. Son mencionables los trabajos realizados [153] en relaci´on con la profundidad Ekman DEy la velocidad en superficie Vo: Vo=C3 $sin |ϕ|U10 ;C3= 0.0068 ; |ϕ|≥10o(2.209) DE=C4 $sin |ϕ|U10 ;C4= 7.12 s (2.210) La primera es la mitad de la propuesta por Ekman, debido en gran parte a las cuestiones de estabilidad aludidas anteriormente; mientras que en la segunda se produce un acierto digno de menci´on dada la complejidad de estas medidas. El transporte Ekman, con sus limitaciones sobre las condiciones de contorno del problema, ser´a de gran utilidad cualitativa y cuantitativa en la metodolog´ıa y anal´ısis de nuestro modelo. 7. ECUACI ´ ON DE ADVECCI ´ ON–DIFUSI´ ON: TRANSPORTE DE MATERIA Esta ecuaci´on describe de una manera general los fen´omenos f´ısicos de transferencia de materia en un sistema f´ısico. En otros campos de la f´ısica es denominada frecuentemente como ecuaci´on de convecci´on–difusi´on o ecuaci´on de trasporte. En la mec´anica de fluidos se ha abordado siempre el problema de transporte de una funci´on escalar A(r, t), que f´ısicamente se corresponde con alguna funci´on densidad (densidad de materia asociada a una sustancia concreta –concentraci´on–, densidad de carga, actividad –densidad de desintegraciones por segundo–, etc.), denoninaci´on que usaremos en algunas partes de este documento. Las partes advectiva y difusiva del transporte est´an asociadas respectivamente al t´ermino promedio  Vy al turbulento (v)delavelocidad(v) del fluido en cada punto P 68 CAP´ ITULO 2 del espacio. Representemos el flujo mediante el vector  Φ,que nos da cuenta de la cantidad de materia (expresada en unidades extensivas tales como gramos o bequerelios) que por unidad de tiempo atraviesa una unidad de superficie colocada en Pperpendicularmente al vector  Φ . Su parte advectiva ser´a igual a  Φadv =A V. (2.211) La parte difusiva –en el problema unidimensional–, siguiendo la primera ley de Fick, toma la expresi´on  Φdif =−D ∇A, (2.212) donde Des la difusividad o coeficiente de difusi´on. Teniendo en cuenta la ecuaci´on de continuidad para el escalar A(r,t), ∂A ∂t + ∇· Φ=S(2.213) donde Ses la funci´on fuente volum´etrica, y sustituyendo el flujo total como suma del advectivo m´as el difusivo  Φ= Φadv + Φdif =A V−D ∇A, (2.214) se llega a ∂A ∂t + ∇·A V+ ∇·−D ∇A=S,(2.215) resultado que incluye como un caso particular (despreciables la advecci´on y las fuentes :  V= 0,S= 0) a la segunda ley de Fick. En el caso del transporte de un radion´uclido, como el que aqu´ı nos ocupa, 137Cs, de actividad Ay constante de desintegraci´on λ, se cumple que su t´ermino fuente es S=−λA. (2.216) Hidrodin´amica General de un Fluido Oc´eanico 69 En el problema tridimensional se suele separar la difusi´on en una parte horizontal y otra vertical  Φdif =−Dh ∇hA+−Dv ∇vA, (2.217) teniendo de este modo la siguiente ecuaci´on para el transporte 3D de la funci´on escalar A(r,t) ∂A ∂t + ∇·A V+ ∇h·−Dh ∇hA+ ∇v·−Dv ∇vA=S.(2.218) El problema hidrodin´amico y el problema de transporte pueden estar acoplados. El ejemplo m´as sencillo es imaginar que la funci´on escalar estudiada en el problema de transporte fuera la funci´on salinidad, materia que participa activamente en las ecuaciones hidrodin´amicas del sistema por su importancia en la funci´on densidad del fluido. En esta tesis doctoral, las funciones escalares principales cuyo transporte vamos a estudiar son: la densidad de materia en suspensi´on y la actividad de 137Cs en el oc´eano, escalares pasivos ambos, que nos permiten desacoplar la hidrodin´amica del problema de transporte. As´ı, desde un punto de vista conceptual, podemos decir que si resolvemos las ecuaciones hidrodin´amicas del apartado anterior, (2.155), (2.156), (2.157), y obtenemos a partir de ellas el campo de velocidades  V(r,t), podremos sustituirlo en la ecuaci´on diferencial de transporte, (2.218), y de all´ı, tras la integraci´on num´erica correspondiente, deducir la soluci´on del problema que nos ocupa aqu´ı, que no es otro que la funci´on escalar A(r,t). Cap´ıtulo 3 Descripci´on General del Mar B´altico y de sus Procesos F´ısicos 1.SITUACI´ ON GEOGR ´ AFICA, DIMENSIONES, BATIMETR´ IA Y ORIGEN GLACIAL DEL MAR B´ ALTICO Para hacer un primer acercamiento al sistema se presenta a continuaci´on una figura esquem´atica 3.1 donde se aprecia la situaci´on geogr´afica del Mar B´altico y una divisi´on del mismo en compartimentos. El B´altico conecta por mar a una gran parte de los pa´ıses del norte de Europa: Alemania, Dinamarca, Suecia, Finlandia, Rusia, Estonia, Letonia, Lituania y Polonia. El zona central y principal del Mar B´altico es denominada B´altico Estricto y aglutina a cinco de las quince zonas: B´altico Estricto Norte (6), Gotland Cuenca (7 y 8), Bornholm (11) y Arkona (12). En esta zona se encuentran los puntos m´as profundos de nuestro sistema, a ambos lados de la isla de Gotland, principal isla del Mar B´altico. La distribuci´on de zonas del B´altico con el ´area y el volumen de cada una se encuentra en la tabla 3.1. En el mapa batim´etrico del Mar B´altico (figura 3.2) se han representado las isol´ıneas de profundidad correspondientes a {0, 10, 20, 40, 60, 80, 100, 150, 200}m , de gran importancia en las corrientes por constituir el conjunto de condiciones de contorno m´as determinante del problema matem´atico que se tiene que resolver. 70 Descripci´on General del Mar B´altico y de sus Procesos F´ısicos 71 Fig.3.1. Divisi´on del Mar B´altico en quince zonas y su situaci´on geogr´afica. L´ogicamente, su influencia ser´a mayor en los elementos del fluido cercanos a la costa y en las capas inferiores de la columna de l´ıquido. TABLA 3.1: Dimensiones principales de compartimentos del Mar B´altico Tabla de dimensiones ´ Area (×103km2)Volumen (km3)prof. media (m) Bah´ıa de Botnia 37.0 1495 40.4 Mar de Botnia 79.0 4885 61.8 Golfo de Finlandia 29.6 1100 37.2 Golfo de Riga 18.8 408 21.7 B´altico Estricto 209 13000 62.2 Belt Sea y The Sound 20.4 290 14.2 Kattegat 22.0 507 23.0 Skagerrak 32.3 6780 210 TOTAL 448.1 28465 63.5 No se representa el Golfo de Riga en el mapa batim´etrico, porque no va considerarse en este golfo b´altico debido a su car´acter pr´acticamente aislado en las aproximaciones que vamos a manejar. Tampoco se representan las regiones 13, 14 y 15, ya que trabajaremos con la hip´otesis de mar cerrado como se justificar´a en el apartado 5 de este cap´ıtulo. 72 CAP´ ITULO 3 El origen del Mar B´altico es relativamente reciente y su nacimiento como mar se debe a la din´amica geol´ogica de la zona durante la ´ultima glaciaci´on, lo que explica, como veremos m´as adelante, gran parte de sus caracter´ısticas. 150 80 60 60 100 80 80 20 60 80 40 40 100 150 80 150 40 40 40 150 100 100 200 200 80 60 100 150 80 60 60 100 80 80 20 60 80 40 40 100 150 80 150 40 40 40 150 100 100 200 200 80 60 100 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 Bathimetric map of the Baltic Sea km km SWEDEN FINLAND BALTIC REPUBLICS 100 60 80 200 200 100 100 150 40 40 40 150 80 150 100 40 40 80 60 20 80 80 100 60 60 80 150 0 50 100 150 200 250 metres Fig.3.2. Mapa batim´etrico del Mar B´altico, con una resoluci´on horizontal de 20 km, que va utilizarse en las simulaciones de nuestro modelo. La ´ultima glaciaci´on afect´o severamente al B´altico Norte, y el sistema no termin´ode formarse como tal hasta hace unos siete mil a˜nos. Antes de adquirir su actual condici´on de mar salobre (7000 B.P.) pas´o incluso por una fase de lago (Ancylus Lake; [41]) durante Descripci´on General del Mar B´altico y de sus Procesos F´ısicos 73 unos dos mil a˜nos, hasta que se abrieron los estrechos de Dinamarca, en un proceso complejo de restauraci´on de la forma del fondo asociado a la fusi´on de las enormes monta˜nas de hielo durante el final de la glaciaci´on. Fig.3.3. Figuras esquem´aticas de la evoluci´on geol´ogica del Mar B´altico durante su formaci´on en los milenios finales de la ´ultima glaciaci´on. Una consecuencia actual de este pasado glacial es que el fondo del Mar B´altico a´un recupera su forma y sube a un ritmo de varios mil´ımetros al a˜no, mayor cuanto m´as al norte; la horquilla va de –1 mm a–1 en el sur a +9 mm a–1 en el norte. Por todo ello, y resumiendo, la batimetr´ıa es muy suave al sur, con presencia de arenas de todo tipo, e irregular y escarpado a medida que avanzamos hacia el norte como recuerdo de su pasado glacial. 80 CAP´ ITULO 3 En un a˜no meteorol´ogico medio, seg´un estas distribuciones, las frecuencias para vientos intensos, muy intensos, y extremos quedan encuadradas en la siguiente tabla orientativa, la tabla 3.2. TABLA 3.2: Tabla orientativa de frecuencias de vientos de intensos, muy intensos y extremos. U10 (m s–1)10 12–15 15–20 20–25 Frecuencia ∼1/semana ∼1/mes ∼1/a˜no ∼1/(5 a˜nos) Ser´a interesante tambi´en tener en cuenta estad´ısticas sobre los vientos dominantes durante cada mes. En el trabajo de Radziejewska [152] se recoge informaci´on a este respecto de a˜nos los 1984–1989, muy ´util para nuestras simulaciones, interesadas en el entorno temporal del accidente nuclear de Chernobyl, ocurrido en el a˜no 1986. 5.DISTRIBUCI´ ON DE LA SALINIDAD EN EL MAR B´ ALTICO El Mar B´altico es un sistema salobre con cierta estructura de gran estuario y una ca´ıda progresiva en la salinidad a medida que nos alejamos del oc´eano y nos adentramos hacia aguas m´as septentrionales. Fig.3.5. Funci´on salinidad del Mar B´altico en funci´on de la distancia al Oc´eano Atl´antico. Curva de sistema decrecientemente salobre [180]. Descripci´on General del Mar B´altico y de sus Procesos F´ısicos 81 La figura 3.5 muestra valores t´ıpicos y promediados en la columna de l´ıquido para la salinidad. El mar B´altico es un mar semicerrado con poco intercambio con el oc´eano. Las zonas denominadas Skagerrak y Kattegat (frente a las costas de Dinamarca) hacen alusi´on a los estrechos primeros que sirven de transici´on entre el Oc´eano Atl´antico y el Mar B´altico. Estos estrechos limitan bastante el intercambio de agua entre el mar y el oc´eano, siendo el tiempo t´ıpico de residencia –cuya definici´on es similar al tiempo de semivida de un is´otopo radiactivo– de una part´ıcula en el Mar B´altico del orden de treinta a˜nos [178, 180] por lo que en esta tesis doctoral se trabajar´abajolahip´otesis de sistema cerrado respecto al Oc´eano Atl´antico. Las simulaciones que van a realizarse en este trabajo estar´an en el orden de uno a diez meses. Describamos brevemente la variaci´on de la salinidad con la profundidad, de la que se muestra un esquema ilustrativo en la figura 3.6. La picnoclina –franja de la columna de l´ıquido con una variaci´on importante en la densidad– y la haloclina –franja de la columna de l´ıquido con una variaci´on importante en la salinidad– permanentes son pr´acticamente coincidentes en el mar B´altico. Fig.3.6. Perfil vertical t´ıpico de salinidad en el Mar B´altico [133]. La profundidad de la permanente no ocupa la misma posici´on en las diferentes zonas del B´altico. Otro hecho destacable es su car´acter despreciable en el norte del B´altico 82 CAP´ ITULO 3 (Mar de Botnia y Bah´ıa de Botnia), por ser la zona m´as alejada a la entrada de agua salada proveniente del Oc´eano Atl´antico. En la tabla 3.3 se presentar valores num´ericos del perfil de salinidad. TABLA 3.3: Posici´on, espesor y salinidad de la picnoclina en diferentes zonas del B´altico Zona D1(m) D2(m) Ssur (PSU) Sbot (PSU) Arkona 40 50 8.5 17 Bornholm 40 70 816 B´altico Centro 60 90 7.5 11.5 Golfo de Finlandia 40 70 6.5 10 Mar de Botnia 40 70 6 7 Bah´ıa de Botnia 40 70 3.5 4 En los c´alculos del modelos usaremos la tabla 3.3, sobre la profundidad (m) de la capa picnoclina en diferentes zonas del Mar B´altico elaborada en esta tesis a partir de la literatura [133]. TABLA 3.4: Profundidades t´ıpicas de la picnoclina permanentes seg´un la zona. Zona profundidad (m) de la picnoclina Arkona 40 Bornholm 40 Gotland Cuenca 60 Golfo de Finlandia 40 En la tabla 3.4 se han eliminado los compartimentos del norte (Mar de Botnia, Bah´ıa de Bothnia) ya que la debilidad de su picnoclina hace que podamos despreciarla en la modelizaci´on, como se ver´aenelcap´ıtulo 5, a tenor de la buena mezcla que se produce en las columnas de l´ıquido de estas zonas del B´altico. Descripci´on General del Mar B´altico y de sus Procesos F´ısicos 83 6.LA INFLUENCIA DEL HIELO EN LA CIRCULACI´ ON DEL MAR B ´ ALTICO Habr´an de tenerse en cuenta las condiciones especiales del sistema durante el invierno, pues de una manera parcial se produce la congelaci´on de los primeros cent´ımetros de la capa superficial del Mar B´altico, especialmente en las latitudes m´as septentrionales (figura 3.7). Como descripci´on general [103] del comportamiento de las capas de hielo, cabe destacar que salvo una delgada l´ınea de hielo que suele adherirse firmemente a la costa, en condiciones normales, el hielo interior se encuentra sometido a la acci´on del viento y de las corrientes sobre las que descansa. Tenemos incluso una clasificaci´on del hielo seg´un sea la ratio entre la superficie helada y la superficie total: hielo muy cerrado (9/10), cerrado (7/10, 8/10), abierto (4/10, 5/10, 6/10) y muy abierto (1/10, 2/10, 3/10). El hielo, en sus diferentes niveles de ratio se encuentra flotando y a la deriva, existiendo incluso modelos de deriva de hielo [195, 197], validados por observaciones de sat´elite. Unos trabajos de los profesores Omstedt [135, 134], Lepp¨aranta [113] y Zhan [209] se˜nalan que la influencia de la capa fina de hielo (10–50 cm) sobre la circulaci´on es menor de lo que esperamos a priori. Seg´un estos c´alculos, basados en una aproximaci´on de circulaci´on inducida por vientos, validada experimentalmente mediante series temporales de la elevaci´on del mar, la atenuaci´on de los flujos respecto a las condiciones sin hielo est´an claramente por debajo del cincuenta por ciento. Al contrario, en el trabajo de Lau [108] se concluye que los flujos de corrientes bajo el hielo est´an por encima (entre un quince y un treinta por ciento) de los flujos con superficie libre. Hay dos factores que favorecen este fen´omeno: (a) P´erdida nula de energ´ıa en la parte r´ıgida del sistema, salvo en los casos de fracturas por choques entre hielos. (b) Aumento de la superficie eficaz absorbente de energ´ıa atmosf´erica cuando la rugosidad del hielo superficial es alta. Otra referencia interesante, seg´un Kowalik [100], sobre circulaci´on bajo hielo en el Mar de Beufort (Oc´eano ´ Artico), tambi´en muestra t´ımidas diferencias entre las frecuencias de las ondas seiche con o sin hielo, comprendidas entre el 0 y el 10 por ciento. 84 CAP´ ITULO 3 Fig.3.7: Distribuci´on de hielo en el Mar B´altico en condiciones de invierno, suave y severo [135]. NAO (Norte Atlantic Oscillation) es un fen´omeno clim´atico que mide las fluctuaciones de la diferencia de presi´on atmosf´erica a nivel del mar entre dos puntos representativos de sendos centros semi-permanentes de presiones como son el Anticicl´on de las Azores y la Borrasca de Islandia; el ´ındice NAO nos indica la diferencia de presi´on atmosf´erica entre dos puntos de referencia, en unas unidades normalizadas. En esta tesis, se considerar´a la posibilidad de inviernos suaves, normales y severos, y se utilizar´a la funci´on NAO (Oscilaci´on del Atl´antico Norte) como variable meteorol´ogica Descripci´on General del Mar B´altico y de sus Procesos F´ısicos 85 auxiliar cuyos valores t´ıpicos [110] para esos tres tipos de inviernos son respectivamente {+2,+1,−1}. NAO es una funci´on del tiempo, manej´andose en Meteorolog´ıa una amplia variedad de promedios con filtros de diferentes escalas: semanal, mensual, anual y niveles superiores, seg´un exija el problema bajo estudio. La diferencia de presiones entre los dos polos del oscilador controla en buena parte la intensidad, direcci´on de los vientos y la din´amica atmosf´erica general en la zona atl´antica norte. Su variaci´on no obedece a una periodicidad particular ni especial: NAO positivo favorece los vientos del oeste y el bloqueo de aire fr´ıo en las zonas ´articas; NAO negativo lleva asociada mayores presiones en las zonas polares y aire de componente norte en las latitudes intermedias. Las condiciones m´as probables, o normales, nos dan altas presiones en la zona subtropical (Azores) y bajas en la zona subpolar (Islandia, Mar B´altico), lo que indica NAO positivo. Fig.3.8: Evoluci´on de la funci´on NAO (Oscilaci´on del Atl´antico Norte) entre 1980 y 1996 [110]. El fen´omeno NAO fue descubierto por G.T. Walker (1868-1958), que tras finalizar sus estudios en Cambridge se emple´o en el Departamento de Meteorolog´ıa de La India, investigando otros fen´omenos clim´aticos importantes como los Monzones del Oc´eano ´ Indico, la Oscilaci´on del Sur (”El Ni˜no”), el posible acoplamiento entre ambos, y el papel de este ´ultimo como indicador de algunos aspectos de la meteorolog´ıa global del planeta. 86 CAP´ ITULO 3 La l´ınea verde (l´ınea m´as suave) de la figura 3.8 representa la evoluci´on de la funci´on NAO invernal (JFM index) entre los a˜nos 1980 y 1996. Cabe destacar los inviernos severos {1985, 1987, 1996}ylossuaves{1989, 1990}; el resto pueden catalogarse como normales. Podemos utilizar como primera aproximaci´on para nuestro modelo, y bas´andonos en el trabajo de Omstedt [135], la siguiente tabla que relaciona la duraci´on en d´ıas de la capa de hielo en los diferentes compartimentos en funci´on de la severidad del invierno. TABLA 3.4: Tiempos de permanencia (d´ıas) del hielo en el Mar B´altico Regi´on b´altica invierno suave invierno normal invierno severo Bah´ıa de Botnia 90 120 150 Mar de Botnia 090 120 Golfo de Finlandia 90 120 150 B´altico Estricto 0 0 90 Otra cifra estimada a partir de este trabajo [135] es la fecha central de los intervalos helados, que podemos situar alrededor del 15 de febrero. Evidentemente, la anterior tabla es una aproximaci´on que admite muchos matices que han de incorporarse en futuras aproximaciones. Se deben se˜nalar trabajos muy detallados [73, 95, 96] que procuran esclarecer y desarrollar las correlaciones entre la funci´on NAO y la distribuci´on de hielo en el oc´eano durante el invierno. En la misma l´ınea, modelos acoplados atm´osfera-hielo-oc´eano [162, 122, 123] analizan situaciones de fuertes vientos yr´apida reducci´on de la concentraci´on de hielo; y situaciones de vientos d´ebiles donde usualmente se produce formaci´on de hielo. Para consultar temperaturas t´ıpicas mensuales en el Mar B´altico se pueden usar los anuarios de temperaturas de los institutos meteorol´ogicos de Suecia y Finlandia o el trabajo recopilatorio de Landsberg [107]. En la figura 3.9 se pueden observar las gr´aficas sobre las temperaturas medias mensuales medidas en una estaci´on meteorol´ogica representativa de cada zona del B´altico: B´altico Estricto Sur, B´altico Estricto Norte, Golfo de Finlandia, Mar de Botnia, Bah´ıa de Botnia. Descripci´on General del Mar B´altico y de sus Procesos F´ısicos 87 0 2 4 6 8 10 12 −10 −5 0 5 10 15 20 month of the year Monthly mean temperature ( o C ) Baltic Proper south Baltic Proper north Gulf of Finland Bothnian Sea Bothnian Bay Fig.3.9. Curvas de temperatura media mensual [107] en diferentes zonas del Mar B´altico. En condiciones de mar helado, y a la luz de la literatura cient´ıfica, nuestras simulaciones podr´ıan trabajar con el siguiente coeficiente de correcci´on λice  Vcor =λice  V,0.7<λ ice <1.2,(3.5) para el campo de velocidades, que nos da la velocidad en condiciones heladas  Vcor en relaci´on a la velocidad en condiciones normales  V. Es destacable el hecho de que la incertidumbre asociada a λice es menor que el asociado al par´ametro de arrastre λD atm´osfera-oc´eano que seg´un el cap´ıtulo 2 proporciona valores t´ıpicos para la velocidad en superficie que oscilan entre el 1.5% y el 3% de la del viento, incertidumbre, que junto a lo costoso que ser´ıa el ajuste de las velocidades por el efecto del hielo, nos invita a trabajar por defecto sin corregir el campo de velocidades, considerando por tanto λice =1. 88 CAP´ ITULO 3 El hielo totalmente cerrado (10/10) suele darse en las zonas costeras y forma l´ıneas paralelas a la costa de una anchura muy delgada [114] y generalmente despreciable frente alaresoluci´on horizontal (20 km) del modelo presentado en esta tesis doctoral. Es altamente improbable que la zona de hielo completamente cerrado complete todo un compartimento: s´olo ocurri´o una vez durante toda la segunda mitad del siglo XX, y fue durante unas semanas en el compartimento m´as septentrional (Bah´ıa de Botnia) del sistema, en el severo invierno de 1985. A´un en compartimentos en esas condiciones sever´ısimas, el grado de transferencia energ´etica entre la atm´osfera y el sistema hielo-oc´eano es relativamente alto seg´un la literatura [114] encontr´andose por encima del cincuenta por ciento. Una l´ınea interesante de futuro de esta tesis doctoral consiste en un acoplamiento robusto de nuestro modelo de dispersi´on con los modelos meteorol´ogicos a corto y medio plazo en el tiempo. Cap´ıtulo 4 La Circulaci´on en el Mar B´altico y su Modelizaci´on 1.INTRODUCCI´ ON. MODELOS HIDRODIN´ AMICOS PREVIOS DEL MAR B´ ALTICO La modelizaci´on hidrodin´amica de esta tesis doctoral tiene como objetivo alimentar a un modelo de transporte posterior cuya finalidad principal es cuantificar la evoluci´on del radion´uclido 137Cs en el Mar B´altico. Estamos interesados aqu´ı en modelar la dispersi´on de 137Cs en dos escalas de tiempos diferentes: (a) ESCALA 1: Escala mensual y corrientes estacionarias. Este problema estudia la evoluci´on de una funci´on escalar en el oc´eano en una escala de tiempo que va desde varias semanas a pocos meses. Tendremos en cuenta el viento dominante, y resolveremos problemas transitorios como el proceso de mezcla en la capa superficial y el arrastre del centroide de una distribuci´on de contaminante en la superficie del mar. Aqu´ı encajan bien un tipo de estudios hidrodin´amicos previos sensibles a la influencia del viento dominante. Dentro de esta rama, destacan los estudios sobre corrientes estacionarias inducidas por el viento, disciplina cl´asica en oceanograf´ıa desde los estudios pioneros de Ekman de 1905 [30]. Centr´andonos en el B´altico, destaca el primer estudio de esta naturaleza publicado por Kowalik en 1972 [99], autor del cl´asico libro mundialmente citado sobre modelizaci´on oce´anica [101], queensuprimeraetapacient´ıfica investig´o las corrientes estacionarias del mar B´altico, labor que continu´odespu´es [79] con otros colaboradores como Jankowski, y 89 96 CAP´ ITULO 4 zk−1 zk∂U ∂x +∂V ∂y +∂W ∂z dz =0.(4.8) Va a darse a continuaci´on el primer paso para integrar las ecuaciones anteriores en el espesor de cada capa (hk=zk−1−zk,k=1, ..., N), seg´un los procedimientos seguidos en la literatura [168, 169, 170] sobre los que nos apoyamos aqu´ı. ·Notaci´on de los t´erminos integrados en la columna de l´ıquido zk−1 zk aloc idz ≡azloc i,k ,zk−1 zk acon idz ≡azcon i,k ,(4.9) zk−1 zk fp idz ≡fzp i,k ,zk−1 zk fcor idz ≡fzcor i,k ,zk−1 zk fR idz ≡fzR i,k (4.10) ·Integraci´on de la componente Udel t´ermino de movimiento: azloc 1,k +azcon 1,k =hk⎡ ⎣∂U ∂t + 2  j=1 ∇j(VjU)⎤ ⎦k−1 2 +(UW)k−1−(UW)k(4.11) donde, los sub´ındices en los corchetes o par´entesis hacen referencia a la posici´on zken la columna de l´ıquido donde se est´a evaluando para la integraci´on. Si adem´as de sub´ındice aparece super´ındice, por ejemplo (UW)k−1 k,esta notaci´on querr´a indicar, (UW)k−1 k= (UW)k−1−(UW)k.En los extremos de la columna de l´ıquido, aparecen las velocidades verticales en la superficie (k−1 = 0) relacionadas con la funci´on elevaci´on, o su valor en el fondo oce´anico (k=C), necesariamente nulo: (W)k−1=0 =˙ ζ,(W)k=C=0.(4.12) N´otese tambi´en que el sumatorio en jvaahora1a2ynode1a3,comoocurr´ıa antes de integrar; esto se debe a que el t´ermino j= 3 ha sido ya integrado. ·De la presi´on y la densidad obtuvimos (2.145, 2.146, 2.154) lo siguiente: ρ=ρ0+ρ,1 ρ≈1 ρ0 ,p=pa+gρ0(ζ−z)+g0 zρdη , (4.13) La Circulaci´on en el Mar B´altico y su Modelizaci´on 97 ∇ip=∂p ∂xi =∂pa ∂xi +gρ0 ∂ζ ∂xi +g0 z ∂ρ ∂xi dη (4.14) ·Integraci´on de la componente Udel t´ermino de fuerzas: fzp 1,k =zk−1 zk−1 ρ ∂p ∂xdz =1 ρ0Pk x+Sk x,(4.15) donde se definen a continuaci´on Pk xySk x. Como ζypano son funciones de z Pk x≡−zk−1 zk (∇xpa+gρ0∇xζ)dz =−hk(∇xpa+gρ0∇xζ) (4.16) PH x≡−0 −H(∇xpa+gρ0∇xζ)dz =−H(∇xpa+gρ0∇xζ) (4.17) Por otro lado Sk x≡−gzk−1 zk dz 0 z∇xρdη , SH x≡−g0 −Hdz 0 z∇xρdη (4.18) Las dos primeras Pk x,PH xson funciones de la presi´on atmosf´erica y de la elevaci´on, mientras que las dos ´ultimas Sk x,SH xdependen de la estratificaci´on horizontal del fluido. Estas dos ´ultimas son las funciones dependientes del gradiente horizontal de densidad (salinidad) calculadas a priori en este modelo por suponerse estacionarias las distribuciones salinas y de temperatura. 2.4. Condiciones de contorno Teniendo en cuenta que zk−1 zk Av ∂2U ∂z2dz =Av ∂U ∂z k−1−Av ∂U ∂z k≡notaci´on ≡Av ∂U ∂z k−1 k (4.19) los t´erminos de Coriolis y Reynolds ser´an: fzcor 1,k +fzR 1,k =hk(fV)k−1 2+hkAh∂2U ∂x2+∂2U ∂y2k−1 2 +Av ∂U ∂z k−1 k (4.20) 98 CAP´ ITULO 4 donde, al evaluar en los extremos, aparecen las tensiones sobre la superficie (k−1=0) o sobre el fondo oce´anico (k=C) en forma de condiciones matem´aticas de contorno ρAv ∂U ∂z k−1=0 =Ts xz ,ρAv ∂U ∂z k=C =Tb xz ,(4.21) sin olvidar tambi´en la necesidad de imponer la nulidad de la velocidad en la direcci´on normal al contorno del sistema. Las tensiones Ts xz yTb xz ser´an modeladas en apartado 2.7 de este cap´ıtulo, apareciendo la influencia del viento y del contacto fluido-fondo. ·Reuniendo todos los t´erminos para la componente U: azloc 1,k +azcon 1,k =hk⎡ ⎣∂U ∂t + 2  j=1 ∇j(VjU)⎤ ⎦k−1 2 +(UW)k−1 k(4.22) fzp 1,k =1 ρ0Pk x(pa,ζ)+Sk x(ρ),(4.23) fzcor 1,k +fzR 1,k =hk(fV)k−1 2+hkAh(ΔhU)k−1 2+Av ∂U ∂z k−1 k .(4.24) 2.5. Esquema conceptual de la integraci´on en el tiempo Los esquemas de Simons [170] y Funquist [43] est´an adaptados a las dificultades num´ericas de estas ecuaciones diferenciales de Navier–Stokes que obligan a trabajar con diferentes pasos de tiempo, mallas alternadas y t´ecnicas de suavizaci´on, entre otras herramientas, para vencer las m´ultiples dificultades que pueden presentar estas ecuaciones. En el apartado 2.9 daremos m´as detalles de esta metodolog´ıa. Lo que toca a continuaci´on, es un esquema conceptual, una especie de gu´ıa simplificada esta metodolog´ıa [170, 43], desnuda de much´ısimos matices num´ericos, y con una intenci´on clarificadora en esta parte de la exposici´on. La integraci´on fundamental de las ecuaciones diferenciales de movimiento es una integraci´on en el tiempo. Conocidas las funciones {U(t),V(t),W(t),ζ(t)}en el instante t, podemos plantear un esquema muy simple para el c´alculo de su valor en el instante t+Δtque ser´ıa el que sigue: U(t+Δt)=U(t)+Δt hk−azloc 1,k +fzp 1,k +fzcor 1,k +fzR 1,k≡FU(U, V, W, ζ;t) (4.25) La Circulaci´on en el Mar B´altico y su Modelizaci´on 99 donde se ha definido FU, una funci´on conocida en el instante t, por depender del conjunto de funciones {U(t),V(t),W(t),ζ(t)}y de derivadas espaciales (∂/∂x, ∂2/∂x2, ∂/∂y , ...) que operan sobre dicho grupo. De una manera m´as abreviada, se tendr´ıa: U(t+Δt)=FU( V,ζ;t) donde  V=[U, V, W].(4.26) An´alogamente se definen FVyFWen (4.30) y (4.31), para las otras componentes del campo  V(r, t). ·Integraci´on de la componente Vde las ecuaciones N–S. Las ecuaciones resultantes son textualmente las mismas que las anteriores, cambiando Upor V,xpor y,adem´as del t´ermino de Coriolis fV que pasa a ser −fU . Reuni´endolas todas se llega a: azloc 2,k +azcon 2,k =hk⎡ ⎣∂V ∂t + 2  j=1 ∇j(VjV)⎤ ⎦k−1 2 +(VW)k−1 k(4.27) fzp 2,k =1 ρ0Pk y(pa,ζ)+Sk y(ρ),(4.28) fzcor 2,k +fzR 2,k =hk(−fU)k−1 2+hkAh(ΔhV)k−1 2+Av ∂V ∂z k−1 k .(4.29) ·Integraci´on en el tiempo de las ecuaci´on de movimiento para V. An´alogamente: V(t+Δt)=V(t)+Δt hk−azloc 2,k +fzp 2,k +fzcor 2,k +fzR 2,k≡FV( V,ζ;t).(4.30) 2.6. Integraci´on de la ecuaci´on de continuidad. Corrientes verticales zk−1 zk∂U ∂x +∂V ∂y +∂W ∂z dz =0 =⇒hk∂U ∂x +∂V ∂y k−1 2 +(W)k−1−(W)k=0, (W)k−1=(W)k−hk∂U ∂x +∂V ∂y k−1 2 ≡FW(4.31) 100 CAP´ ITULO 4 donde (W)k−1=0 =˙ ζ,(W)k=C=0.(4.32) Veremos en el apartado 2.8 que el esquema operar´a siempre de modo que FWsea una funci´on conocida, cuesti´on que se consigue empezando los c´alculos por la capa de fondo, aprovechando que es nula la velocidad vertical en su cara inferior, coincidente con la superficie de contacto fluido–fondo. 2.7. Modelo de las condiciones de contorno y de la viscosidad turbulenta ·Tensiones superficiales y de fondo En el trabajo de referencia que aqu´ı nos ocupa [43] y en su antecesor [170] las expresiones usadas para las tensiones superficial y de fondo son respectivamente las siguientes: −→ Ts=CDρaUa Ua,−→ Tb=BρH V2D,(4.33) donde CDes el coeficiente de arrastre, ρaes la densidad del aire,  Uaes la velocidad del aire a una altura de 10 metros sobre la superficie del mar, Bes un coeficiente dependiente de la profundidad H,y V2Des el promedio de  Ven la columna de l´ıquido, cuyos detalles num´ericos son: CD=1.6×10−3,ρ a=1.3kgm –3 ,B=b/H2,b=5×10−3m2s–1 .(4.34) Una parametrizaci´on de CDmuy usada [101, 143, 46], como funci´on de la velocidad del viento Uaes CD(Ua)=(γ1+γ2Ua)×10−3,γ 1=0.75 ,γ 2=0.067 .(4.35) Los modelos de referencia [170, 43] han trabajado con valores de Avque en la superficie han estado en el rango Av∼50 −100 cm2s–1 , mientras que en la zona de altos gradientes de densidad se han usado valores muy peque˜nos, Av<1cm 2s–1 .En lo que se refiera a la viscosidad turbulenta horizontal, los autores [170, 43] han considerado en sus modelos valores de Ah∼100 m2s–1 ,habiendo hecho c´alculos considerando Ah>100 m2s–1 sin que se noten cambios dignos de menci´on en los resultados. La Circulaci´on en el Mar B´altico y su Modelizaci´on 101 2.8. Algoritmo general ascendente por capas Una vez hecha una integraci´on en el tiempo con unas condiciones iniciales, habr´ıa que iterarla adecuadamente, aprovechando las condiciones de contorno. Vuelve a exponerse a continuaci´on un esquema sencillo, pero que b´asicamente reproduce la esencia del m´etodo usado en los referidos trabajos [170, 43]. Est´a dividido en diferentes pasos, a modo de bucle. (a) PASO 0: Condiciones iniciales. Inicialmente (t=t0) conocemos las funciones {U, V, W, ζ}quedependendelaposici´on, (x, y, z;t=t0) ; de una manera m´as abreviada,  V= V(r;t=t0) : Condiciones iniciales . (4.36) Seguidamente (Paso 1), ha de notarse que las variables conocidas {U, V, W, ζ}(t0) aparecen dentro de FUyFV, funciones diferenciales en el espacio R3, resultantes de separar UyVpara integrar en el tiempo. (b) PASO 1: t=t0+Δt, capa de fondo k=N. Del apartado anterior tenemos las ecuaciones que permiten obtener expl´ıcitamente UN−1 2(t0+Δt)=FU( V,ζ;t0),k=N(4.37) VN−1 2(t0+Δt)=FV( V,ζ;t0),k=N(4.38) que sustituidas a continuaci´on en la versi´on finita de la ecuaci´on de continuidad, nos dan la velocidad vertical en la cara superior de la caja WN−1(t0+Δt)=FW(U, V, WN;t0+Δt),k=N−1 (4.39) gracias a que permanentemente Wk=N(t)=0,∀ten la cara inferior de la misma. (c) PASO 2: t=t0+Δt, capa superior a la de fondo k=N−1. An´alogamente se pueden obtener expl´ıcitamente UN−3 2(t0+Δt)=FU( V,ζ;t0),k=N−1 (4.40) 102 CAP´ ITULO 4 VN−3 2(t0+Δt)=FV( V,ζ;t0),k=N−1 (4.41) que sustituidas a continuaci´on en la versi´on finita de la ecuaci´on de continuidad, nos dan la velocidad vertical en la cara superior de la caja WN−2(t0+Δt)=FW(U, V, WN−1;t0+Δt),k=N−2 (4.42) gracias a que conocemos WN−1(t0+Δt) en la cara inferior de la misma, como salida del paso anterior. Este paso se repite en la capa superior (k=N−2), y se reitera hasta llegar a la pen´ultima (k= 2), inclusive. (d) PASO 3: t=t0+Δt, capa superficial k=1. Del mismo modo se pueden obtener expl´ıcitamente U1 2(t0+Δt)=FU( V,ζ;t0),k= 1 (4.43) V1 2(t0+Δt)=FV( V,ζ;t0),k= 1 (4.44) que sustituidas a continuaci´on en la versi´on finita de la ecuaci´on de continuidad, nos dan la velocidad vertical en la superficie del fluido W0(t0+Δt)=FW(U, V, W1)|t0+Δt=˙ ζ(t0+Δt),k=0,(4.45) ζ(t0+Δt)=desarrollo en serie ≡Fζ˙ ζ,ζ;t0.(4.46) gracias a que conocemos W1(t0+Δt) enlacarainferiordelacapasuperficial,como salida del paso anterior. De este modo cerramos el ciclo, puesto que ya conocemos los campos  V(r;t=t0+Δt)yζ(x, y;t=t0+Δt) (4.47) llegando as´ı a las condiciones de partida, Paso 0, pero en el instante num´erico posterior (t=t0+Δt): Fin del ciclo. La Circulaci´on en el Mar B´altico y su Modelizaci´on 103 Aqu´ıseiniciar´ıa otro ciclo, exactamente igual que el anterior, que transcurrir´ıa entre t=t0+Δtyt=t0+2Δt;yas´ı sucesivamente hasta completar la integraci´on num´erica del problema. 2.9. M´as detalles del esquema num´erico de Simons [170] y del SMHI [43] para las ecuaciones N–S (a) Se trata de un esquema de diferencias finitas centradas. (b) La computaci´on empieza por la capa de fondo y va sucesivamente subiendo hasta la capa superficial. (c) Para tener en cuenta las soluciones oscilatorias, los pasos de tiempo tienen que ser suficientemente peque˜nos y adecuadamente combinados con los pasos espaciales, seg´un un criterio de Courant–Friedrichs–Lewy para ecuaciones diferenciales de tipo hiperb´olico, como es el caso: Δt Δs<1/√2 cmax (4.48) donde cmax es la velocidad de onda m´as r´apida en el modelo, que este caso el la onda gravitatoria externa cuyo valor es cmax ∼√gH . (d) Por cuestiones de estabilidad num´erica, se utilizaron pasos de tiempo distintos seg´un los procesos f´ısicos involucrados: el t´ermino del gradiente de presi´on, el de Coriolis y el no lineal han sido tratados con pasos de tiempo menores que los t´erminos asociados a la viscosidad, mediante un esquema de tipo Lax-Wendroff. (e) El t´ermino advectivo no lineal y el de Coriolis requieren un delicado tratamiento tambi´en del mallado espacial con objeto de alcanzar buenas condiciones de estabilidad, lo que se consigue combinando adecuadamente dos reticulados distintos buscando la m´axima estabilidad de la soluci´on. (f) Los t´erminos no lineales producen altos arm´onicos que se han minimizado tanto jugando con los coeficientes asociados a la turbulencia como realizando suavizaciones de soluciones conocidas del problema. (g) Se han hecho diversas simulaciones y se ha comprobado que los niveles medios de elevaci´on se alcanzan al cabo de algunas horas (5–20 horas) en funci´on de la velocidad del viento. 104 CAP´ ITULO 4 (h) La estabilidad num´erica de las soluciones se ha chequeado num´ericamente mediante la evaluaci´on de la energ´ıa total en ausencia tanto de t´erminos forzados como de disipativos. (i) Esta metodolog´ıa num´erica de diferencias finitas centradas en el espacio y en el tiempo fue validada y puesta a punto con varias geometr´ıas sencillas, concluy´endose que los errores por truncamiento con respecto a la soluci´on exacta son del orden de un 1 % . 2.10. Aspectos meteorol´ogicos del modelo del SMHI [43] Se ha dise˜nado un a˜no meteorol´ogico t´ıpico que contiene, con todas las dificultades de elaboraci´on, limitaciones y ventajas que esto conlleva, los eventos meteorol´ogicos m´as probables. 0 5 10 15 20 25 30 35 −10 −5 0 5 10 Wind spectrum for a typical winter month, from the SMHI day of the month Velocity scale (m/s) = 10 kt 5 m/s 10 m/s = 20 kt Figura 4.2: Espectro de vientos para un mes t´ıpico de invierno, dise˜nado por el SMHI [43]. En el eje vertical se representa una escala de la velocidad y en el horizontal el d´ıa del mes. El a˜no t´ıpico incluye las variaciones t´ıpicas de este intervalo de tiempo, de modo que las diferencias entre este a˜no t´ıpico y un a˜no cualquiera son de un orden menor respecto a las primeras.A la hora de hacer las simulaciones va a hacerse uso de la informaci´on La Circulaci´on en el Mar B´altico y su Modelizaci´on 105 meteorol´ogica registrada en el Instituto Meteorol´ogico de Suecia donde se han seleccionado datos acumulados durante cinco d´ecadas que van a ser de gran utilidad. 0 5 10 15 20 25 30 35 −8 −6 −4 −2 0 2 4 6 8 Wind spectrum for a typical spring month, from the SMHI day of the month Velocity scale (m/s) Figura 4.3: Espectro de vientos para un mes t´ıpico de primavera, dise˜nado por el SMHI [43]. En el eje vertical se representa una escala de la velocidad y en el horizontal el d´ıa del mes. 0 5 10 15 20 25 30 35 −6 −4 −2 0 2 4 6 Wind spectrum for a typical summer month, from the SMHI day of the month Velocity scale (m/s) Figura 4.4: Espectro de vientos para un mes t´ıpico de verano, dise˜nado por el SMHI [43]. En el eje vertical se representa una escala de la velocidad y en el horizontal el d´ıa del mes. Para cada una de las cuatro estaciones meteorol´ogicas se ha escogido un mes representativo formado por 25–30 d´ıas, constituido por 4–6 eventos que cubren las caracter´ısticas concretas del invierno, primavera, verano y oto˜no en esta zona de la superficie 112 CAP´ ITULO 4 ζ −H ∂V ∂y dz =∂ ∂y ζ −HVdz   My/ρ0 −Vζ ∂ζ ∂y −VH ∂H ∂y (4.69) donde {Uζ,U H,V ζ,V H}son los valores que toman las componentes horizontales de la velocidad en los contornos superior e inferior del problema. Acerca de la velocidad vertical Wen las partes superior e inferior de la columna de l´ıquido, se tiene lo siguiente: [W]ζ −H=[W]ζ−[W]−H(4.70) Hagamos las siguientes definiciones de funciones: [W]ζ≡Wζ≡Dζ Dt [W]−H≡WH≡D(−H) Dt =−DH Dt (4.71) Recordemos que se trata de funciones de dos variables {x, y}donde Wζ=Wζ(x, y) yWH=WH(x, y), por lo que el c´alculo de su diferencial absoluta puede obtenerse como la derivada local m´as el t´ermino convectivo: Wζ=Dζ Dt =∂ζ ∂t + Vζ· ∇ζW H=−DH Dt =−∂H ∂t  =0 − VH· ∇H(4.72) donde  Vζ; VHrepresentan al campo de velocidades en la parte superior e inferior de la columna de l´ıquido respectivamente. Sustituyendo en la velocidad vertical queda: [W]ζ −H=Wζ−WH=∂ζ ∂t + Vζ· ∇ζ+ VH· ∇H(4.73) Sustituyendo (36),(37) y (38) en la ecuaci´on de continuidad integrada (32): ∂Mx ∂x −ρ0Uζ ∂ζ ∂x +UH ∂H ∂x +∂My ∂y −ρ0Vζ ∂ζ ∂y +VH ∂H ∂y  +ρ0 ∂ζ ∂t +ρ0 Vζ· ∇ζ+ VH· ∇H= 0 (4.74) La Circulaci´on en el Mar B´altico y su Modelizaci´on 113 Reordenando sumandos: ∂Mx ∂x −ρ0Uζ ∂ζ ∂x +Vζ ∂ζ ∂y   ( Vζ· ∇)ζ +∂My ∂y −ρ0UH ∂H ∂x +VH ∂H ∂y    ( VH· ∇)H + +ρ0 ∂ζ ∂t +ρ0 Vζ· ∇ζ+ VH· ∇H= 0 (4.75) Tras anularse cuatro de los sumandos dos a dos, se llega a la expresi´on final para la ecuaci´on de continuidad del fluido oce´anico integrada en la dimensi´on z: ∂Mx ∂x +∂My ∂y +ρ0 ∂ζ ∂t = 0 (4.76) Yas´ı se llega al sistema final de tres ecuaciones diferenciales para las tres inc´ognitas {Mx,M y,ζ}: ∂ ∂tMx−fMy=−(ζ+H)∂pa ∂x +gρ0 ∂ζ ∂x+Sx+Ts xz −Tb xz +AhΔhMx(4.77) ∂ ∂tMy+fMx=−(ζ+H)∂pa ∂y +gρ0 ∂ζ ∂y+Sy+Ts yz −Tb yz +AhΔhMy(4.78) ∂Mx ∂x +∂My ∂y +ρ0 ∂ζ ∂t = 0 (4.79) La soluci´on estacionaria de este sistema 2D puede acoplarse despu´es (apartado 3.6 de este cap´ıtulo) al sistema que plantea las velocidades horizontales {U(x, y, z),V(x, y, z)} en cada capa del mar B´altico. Posteriormente, en el apartado 5.2 de este cap´ıtulo, implementaremos este planteamiento num´ericamente a nuestra malla de c´alculo 3D. 3.3. Parametrizaci´on de las TENSIONES Ts,TbydelaTURBULENCIA{Av,A h} Las tensiones en la superficie del oc´eano Ts xz,Ts yzhan sido estudiadas y parametrizadas en diferentes modelos num´ericos [101, 155, 168]: se va a suponer que el vector de tensiones debido al viento −→ Tslleva la direcci´on del viento −→ Uayquesum´odulo −→ Ts es proporcional al cuadrado de la velocidad −→ Ua.Aqu´ı se muestra su dependencia de la velocidad del viento −→ Ua, de sus componentes horizontales −→ Ua=[Ua,x Ua,y], de su m´odulo Ua, del coeficiente de arrastre CDy de la densidad del aire ρa: 114 CAP´ ITULO 4 −→ Tsua≡−→ Ua −→ Ua≡−→ Ua Ua−→ Ts∝U2 a⇒−→ Ts=ρCDU2 a −→ Ua Ua ⇒Ts xz =ρaCDUaUa,x Ts yz =ρaCDUaUa,y (4.80) donde uaes el unitario en la direcci´on del viento cuya velocidad tiene por m´odulo Ua= −→ Ua=U2 a,x +U2 a,y.Enlosc´alculos se tomar´aρa=1.3kgm –3 yparaelcoeficientede arrastre CD=2.6×10–3 [155, 168]. Existen diversas parametrizaciones en la literatura [101, 143] para CDcomo funci´on de la velocidad del viento. Por otro lado, las tensiones en el fondo del mar Tb xz,Tb yzvan a ser modeladas [101, 155, 168] para esta primera aproximaci´on de manera an´aloga a las de la superficie: el vector de tensiones debido al fondo oce´anico −→ Tblleva la direcci´on del vector transporte de masa −→ Mysum´odulo −→ Tbes proporcional al cuadrado del vector −→ My a la inversa del cuadrado de la profundidad total H+ζ.Aqu´ı se muestra su dependencia respecto a −→ M, su m´odulo y sus componentes, la profundidad y el coeficiente de fricci´on r: −→ Tbub≡−→ M −→ M−→ Tb∝M21 (H+ζ)2⇒−→ Tb=rM21 (H+ζ)2 −→ M −→ M ⇒Tb xz =rM (H+ζ)2MxTb yz =rM (H+ζ)2My(4.81) donde ubes el unitario en la direcci´on del vector −→ Mcuyo m´odulo viene dado por M= −→ M=M2 x+M2 y. El coeficiente de fricci´on rdepende [95] de la rugosidad del fondo y de las propiedades de la capa l´ımite de fondo, alcanzando valores t´ıpicamente en el intervalo 2−4×10–3 .Enlosc´alculos se tomar´ar=2.5×10–3 [155, 168] para el coeficiente de fricci´on. El coeficiente de viscosidad vertical Avva a ser estimado mediante la teor´ıa de Felsenbaum para corrientes inducidas por vientos, en que Avdepende [155, 179] de la velocidad del viento Ua,delpar´ametro de Coriolis fy de la profundidad H: Av=α1UaHsi H<H cr ´o Av=α2 U2 a fsi H ≥Hcr (4.82) La Circulaci´on en el Mar B´altico y su Modelizaci´on 115 donde Hcr =α3Ua/f es una profundidad cr´ıtica, y los coeficientes {α1,α 2,α 3}son iguales a {0.54×10–4,4.7×10–8,8.7×10–4}respectivamente. El par´ametro Hcr para el B´altico es del orden de la profundidad media y los valores que se obtienen para Aven ambas partes de la funci´on a trozos son similares, y ayudar´an en el dise˜no del modelo de dispersi´on presentado en esta tesis doctoral. La viscosidad turbulenta horizontal puede ser parametrizada dependiendo de la escala espacial de los procesos din´amicos intervinientes, teniendo en cuenta la conocida Ley 4/3 de Richardson Ah=coL4/3(4.83) donde Lrepresenta la escala horizontal de los torbellinos presentes en el problema, y co es una constante emp´ırica [26] que tiene un orden de magnitud de co≈510 –4 m2/3s–1. En los c´alculos bidimensionales del modelo del IOPAN [80] se trabaja tomando Ligual a varias veces el tama˜no de la resoluci´on horizontal (hx=hy=10km),tales que Ah> 1000 m2s–1 ,paraas´ı favorecer la obtenci´on num´erica del estado estacionario. 3.4. Detalles sobre la integraci´on num´erica de las ecuaciones 2D Una forma sencilla de integrar este sistema de ecuaciones (4.77)–(4.79) es mediante un esquema de M´etodo de Diferencias Finitas que han seguido diferentes autores [101, 155]: Mt+τ xm+1,n =1−2τRm+1,nMt−τ xm+1,n +2fτMy+2τGxm+1,n + −ρogHm+1,n τ hxζt m+2,n −ζt m,n+2τAhΔhMt+τ xm+1,n (4.84) Mt+τ ym,n+1 =1−2τRm,n+1Mt−τ ym,n+1 −2fτMx+2τGym,n+1 + −ρogHm,n+1 τ hyζt m,n+2 −ζt m,n+2τAhΔhMt+τ ym,n+1 (4.85) ζt+2τ m,n =ζt m,n −τ ρo1 hxMt+τ xm+1,n −Mt+τ xm−1,n +1 hyMt+τ ym,n+1 −Mt+τ ym,n−1(4.86) donde: 116 CAP´ ITULO 4 Gxm+1,n =pam+2,n −pam,n 2hx Hm+1,n +Ts xm+1,n +Sxm+1,n (4.87) Gym,n+1 =pam,n+2 −pam,n 2hy Hm,n+1 +Ts ym,n+1 +Sym,n+1 (4.88) My=Mt−τ ym+2,n+1 +Mt−τ ym+2,n−1+Mt−τ ym,n+1 +Mt−τ ym,n−1/4 (4.89) Mx=Mt−τ ym+1,n+2 +Mt−τ ym−1,n+1 +Mt−τ ym+1,n +Mt−τ ym−1,n /4 (4.90) ΔhMt+τ xm+1,n =Mt+τ xm+1,n+2 +Mt+τ xm+1,n−2+Mt+τ xm+3,n +Mt+τ xm−1,n −4Mt+τ xm+1,n /4h2 x(4.91) ΔhMt+τ ym,n+1 =Mt+τ ym,n+1 +Mt+τ ym,n+3 +Mt+τ ym,n−1+Mt+τ ym+2,n+1 −4Mt+τ ym−2,n+1 /4h2 y(4.92) y teniendo en cuenta las expresiones (23) para la tensi´on en el fondo de la columna de l´ıquido: Rm+1,n =r H+ζt m+2,n +ζt m,n/22My 2+Mt−τ xm+1,n 2(4.93) Rm,n+1 =r H+ζt m,n+2 +ζt m,n/22Mx 2+Mt−τ ym,n+1 2(4.94) El mallado de diferencias finitas utilizado ha sido de tipo Arakawa C. Ya se ha hablado de la importancia de la b´usqueda de soluciones estacionarias de las ecuaciones diferenciales que tenemos que resolver. Para ello, es aconsejable analizar la soluci´on transitoria del sistema (4.77)–(4.79) que se har´a con la ayuda del esquema num´erico mostrado en las expresiones (4.84)–(4.94), y que conducir´a, tras un tiempo determinado, a valores pr´acticamente invariables del transporte de masa y de la elevaci´on en cada punto, es decir de los campos Me x(x, y),M e y(x, y),ζ e(x, y), donde el super´ındice e hace alusi´on su car´acter estacionario. La Circulaci´on en el Mar B´altico y su Modelizaci´on 117 Si tenemos en cuenta que en estado estacionario se cumple: ∂ ∂tMe x=0 ∂ ∂tMe y=0 ∂ ∂tζe= 0 (4.95) y sustituimos en (4.77)–(4.79): −fMe y=−(ζe+H)∂pa ∂x +gρo ∂ζe ∂x +Sx+Ts xz −Tb xz +AhΔhMe x(4.96) fMe x=−(ζe+H)∂pa ∂y +gρo ∂ζe ∂y +Sy+Ts yz −Tb yz +AhΔhMe y(4.97) ∂Me x ∂x +∂Me y ∂y = 0 (4.98) •Funci´on de corrientes Ψ Dado el car´acter solenoidal ( ∇·  Me= 0) del campo estacionario de transporte de masa, podemos definir una funci´on escalar de corrientes Ψ a partir de la cual se pueden calcular las componentes del campo  Me, con las ventajas matem´aticas y gr´aficas –l´ıneas de nivel, grandes circu´ıtos de corriente, etc– que proporciona una funci´on del tipo Ψ(x, y) quecontengatodalainformaci´on anal´ıtica de las corrientes fundamentales que tienen lugar en el sistema. La funci´on de corrientes Ψ(x, y) se define de la siguiente manera: Me x=−∂Ψ ∂y Me y=∂Ψ ∂x (4.99) y si se expresa en forma vectorial:  Me=ez∧ ∇Ψ (4.100) donde ezes el vector unitario en el sentido creciente del eje OZ . 3.5. Descripci´on de algunas soluciones del problema 2D Ya se han introducido los fundamentos de los modelos hidrodin´amico bidimensional (2D). A continuaci´on se van a mostrar una parte de sus resultados para ilustrar el problema y acercarnos al comportamiento f´ısico del sistema. En primer lugar haremos una 118 CAP´ ITULO 4 descripci´on de las corrientes inducidas bajo vientos del Oeste-Este, para terminar con los resultados asociados a los vientos en sentido Sur-Norte. 3.5.1. Descripci´on de corrientes inducidas por vientos del Oeste-Este. Pasamos a describir soluciones del sistema de ecuaciones para corrientes 2D tales que  Vwex;Ua=10 ms –1. Antes de hacer un an´alisis por compartimentos, veremos las gr´aficas globales para el Mar B´altico, empezando por la soluci´on estacionaria de las funciones cuando tenemos las condiciones meteorol´ogicas citadas anteriormente. La gr´afica 4.6 presenta la funci´on de corrientes Ψ(x, y) con isol´ıneas de la forma Ψ=cte expresadas en la unidad Sverdrup (1 Sv = 106m3s–1 =1Hm 3s–1) en honor al insigne ocean´ografo escandinavo. Las l´ıneas de corrientes plumeadas indican regiones de circulaci´on cicl´onica, o lo que es lo mismo en sentido contrario a las agujas del reloj. Puede hablarse de cierto predominio de la circulaci´on cicl´onica, aunque como se ver´am´as detalladamente en las gr´aficas compartimentales, hay regiones importantes con circulaci´on anticicl´onica (l´ıneas lisas sin plumear). El campo de corrientes superficiales  VS(x, y) mostrado en la figura 4.7 es pr´acticamente uniforme en la direcci´on Oeste que es la direcci´on seguida por los vientos de altura  Vwex. Hay una ligera desviaci´on de esta tendencia sobre todo en los puntos de la costa este debido a la influencia de las condiciones de frontera. El mapa ilustrativo (figura 4.8) de las corrientes medias  V2D(x, y) promediadas verticalmente guarda por su propia naturaleza una estrecha relaci´on con el de las l´ıneas de corriente Ψ como se puede ver comparando ambos gr´aficos. Con un poco de intuici´on se vislumbran en la figura que viene a continuaci´on las l´ıneas de corrientes de la figura anterior. Aparecen interesantes detalles que se comentar´an m´as detenidamente en el an´alisis compartimental como es el protagonizado por altas velocidades medias en las zonas costeras del este del Mar B´altico. Resulta interesante tambi´en la representaci´on (figura 4.9) de la funci´on elevaci´on ζ(x, y)enelr´egimen estacionario y condiciones meteorol´ogicas que ya conocemos. Salta la vista la estructura del Mar B´altico en forma de ”Y” griega, con dos brazos, uno principal y otro secundario que da lugar al Golfo de Finlandia, configuraci´on que es la m´axima responsable del car´acter creciente de la funci´on ζ(x, y) hacia al norte en la rama principal y hacia el este en la secundaria. La Circulaci´on en el Mar B´altico y su Modelizaci´on 119 Aeste´ultimo factor hay que sumar tambi´en la influencia de la batimetr´ıa del fondo marino, por la conocida tendencia que tienen a alinearse las curvas de profundidad H(x, y)=cte y las de elevaci´on ζ(x, y)=cte, o mejor dicho la tendencia de las l´ıneas de elevaci´on a alinearse con las l´ıneas de profundidad. Figura 4.6: L´ıneas de la funci´on de corrientes Ψ(Sv) inducidas por vientos del oeste tales que Ua=10ms –1 [80]. 120 CAP´ ITULO 4 Figura 4.7: Corrientes superficiales  VS(x, y) inducidas por vientos del oeste tales que Ua=10ms –1 [80]. La Circulaci´on en el Mar B´altico y su Modelizaci´on 121 Figura 4.8: Corrientes medias  V2D(cm s–1) promediadas verticalmente e inducidas por vientos del Oeste tales que Ua=10ms –1 [80]. 128 CAP´ ITULO 4 Figura 4.12: Corrientes superficiales  VS(x, y) inducidas por vientos del sur tales que Ua=10ms –1 [80]. La Circulaci´on en el Mar B´altico y su Modelizaci´on 129 Figura 4.13: L´ıneas de elevaci´on ζ(cm) inducidas por vientos del sur tales que Ua=10ms –1 [80]. 130 CAP´ ITULO 4 A continuaci´on se van exponer detalles de las soluciones estacionarias de las ecuaciones en cada unos de los cinco compartimentos que forman el sistema. (a) Soluciones en el Golfo de Botnia En este compartimento se tiene una mayor presencia de la circulaci´on cicl´onica, sobre todo en el norte donde se aprecian dos dominios con l´ıneas antihorarias. En el sur el fluido circula m´as d´ebilmente y con sentido anticicl´onico como se muestra en la figura de corrientes medias. La elevaci´on del mar ζ(x, y) es claramente creciente hacia el norte, alcanz´andose en este compartimento los m´aximos del Mar B´altico para estas condiciones meteorol´ogicas, con valores que superan la mitad de un metro. Las corrientes superficiales van n´ıtidamente hacia el norte con ligeras desviaciones en las costas septentrionales en las que han de satisfacerse l´ogicamente las condiciones de contorno del problema. (b) Soluciones en el Mar de Botnia La caracter´ıstica m´as destacable del movimiento del fluido en esta zona es que las corrientes anticicl´onicas cubren la mayor parte del compartimento. Se tiene un gran torbellino horario en el sur y otro de tama˜no intermedio y antihorario en el norte. Hay que mencionar tambi´en la preponderancia de la circulaci´on anticl´onica en la parte central, acompa˜nada de peque˜nos remolinos cicl´onicos en la costa este, apreciables m´as claramente el mapa de corrientes medias. En las l´ıneas de elevaci´on ζ(x, y)=cte se nota la influencia de la topograf´ıa junto la tendencia creciente del campo en el sentido sur-norte, gracias a lo cual la funci´on tiene una forma pr´acticamente plana en la zona central del compartimento. Las variaciones de las corrientes superficiales  VS(x, y) respecto de la direcci´on predominante marcada por el meridiano local se deben a las corrientes de retorno que tienen lugar en las franjas norte y sur del Mar de Bothnia. (c) Soluciones en el Golfo de Finlandia Dadas las condiciones meteorol´ogicas con vientos del sur y un compartimento como ´este con una geometr´ıa tan claramente orientada en la direcci´on oeste-este se obtiene como soluci´on un campo de velocidades no muy intensas que hacen circular al sistema en sentido contrario a las agujas del reloj. Las corrientes superficiales  VS(x, y)van claramente en sentido norte, con las l´ogicas desviaciones en las cercan´ıas del contorno, La Circulaci´on en el Mar B´altico y su Modelizaci´on 131 ylasl´ıneas de elevaci´on presentan esta misma tendencia creciente hacia el norte como consecuencia de la acci´on de los vientos. (d) Soluciones en el B´altico Centro Cuando las corrientes atmosf´ericas soplan hacia el norte y el sistema afectado est´a orientado y canalizado fundamentalmente en esta direcci´on existen un grupo de soluciones matem´aticas a este problema [207, 23, 205, 22] que apuntan en la direcci´on de la formaci´on de dominios de circulaci´on cicl´onica en zona este y anticicl´onica en la oeste, con el acompa˜namiento de corrientes de retorno en las zonas norte y sur. Si a˜nadimos a esta propiedad la presencia de un importante isla en la zona (isla de Gotland) y las irregularidades de la costa de Suecia cercana al Archipi´elago, podremos interpetar bien la soluci´on del problema fluidomec´anico que nos cita aqu´ı. Se observa as´ı circulaci´on anticicl´onica en la costa oeste a la altura de la isla, en la zona oeste da la franja delimitada entre la gran isla y las costas b´alticas, y en el litoral noroccidental no invadido por el gran cicl´on generado frente las costas de Estocolmo que se expande debido a la forma del relieve costero. Hay una fuerte orientaci´on de las corrientes superficiales  VS(x, y) hacia el norte con un cambio de orientaci´on claro en norte (Archipi´elago) donde se producen la corrientes de retorno y los flujos m´as intensos. Las l´ıneas de nivel presentan la tendencia general creciente hacia el norte, mezclada con la influencia de la topograf´ıa del fondo del mar. Otro hecho destacable es relativamente peque˜no intercambio hidr´aulico con el Golfo de Finlandia, mucho mayor cuando los vientos soplan en la direcci´on perpendicular a la que nos ocupa ahora mismo. (e) Soluciones en el B´altico Sur En la parte suroeste (Isla de Bornholm) se tienen dos remolinos en sentido contrario al norte y sur de la peque˜na isla; un cicl´on frente a las costas suecas y anticicl´on frente a las alemanas. El resto del sistema est´a recorrido por corrientes medias  V2Dhacia el norte en la costa, cre´andose circulaci´on horaria en la mitad oeste y antihoraria en la mitad oriental. Son muy notables las corrientes fuertes de retorno en el sur (Golfo de Gdansk) con valores extremos de flujo. Es precisamente en esa zona donde las corrientes superficiales  VSse salen m´as de tendencia natural hacia el norte para adaptarse a la frontera y cumplir as´ı las condiciones de contorno. Esto se da a lo largo de toda la costa sur b´altica. Las isol´ıneas de elevaci´on est´an en la 132 CAP´ ITULO 4 t´onica general de crecimiento modelada por la topograf´ıa de fondo que aparece en el resto del Mar B´altico. 3.6. Segunda fase de integraci´on: Del modelo 2D a las corrientes 3D Se retoman a continuaci´on las ecuaciones generales del cap´ıtulo 2, (2.155) y (2.156). En condiciones bastante frecuentes de no proximidad a la costa, donde se da una relativa suavidad en campo horizontal de corrientes, es razonable considerar la aproximaci´on (AhΔhV1,2AvΔvV1,2), como veremos en la siguiente estimaci´on con valores t´ıpicos para el esquema num´erico del mar B´altico: ΔhU∼∂2U ∂x2∼ΔU h2 x∼U h2 x∼10−1ms –1 108m2∼10−9m–1 s–1 (4.101) ΔvU∼∂2U ∂z2∼ΔU h2 z∼U h2 2∼10−1ms –1 102m2∼10−3m–1 s–1 (4.102) Ah∼102m2s–1 ,A v∼10−2m2s–1 =⇒AvΔvU∼100 AhΔhU: (4.103) AhΔhV1,2AvΔvV1,2(4.104) Tras un tiempo suficientemente grande como para que el viento haya conseguido un r´egimen de circulaci´on suficientemente estable (∂V1,2/∂t ≃0) puede buscarse una soluci´on estacionaria al problema de las corrientes horizontales: Av ∂2U ∂z2+fV =1 ρo∂pa ∂x +gρo ∂ζ ∂x +g0 z ∂ρ ∂x dη(4.105) Av ∂2V ∂z2−fU =1 ρo∂pa ∂y +gρo ∂ζ ∂y +g0 z ∂ρ ∂y dη(4.106) similar al caso simple de Ekman estudiado en el cap´ıtulo 2 (2.171, 2.172), ecuaciones homog´eneas que generaban all´ı la espiral del mismo nombre. Las que nos ocupan aqu´ı (4.105, 4.106) tienen, a diferencia de las primeras, un t´ermino no homog´eneo a la derecha, dependiente sobre todo de la funci´on elevaci´on y del gradiente horizontal de salinidad. La soluci´on del sistema formado por (4.105) y (4.106) puede obtenerse anal´ıticamente conocidas las funciones {Az,p a,ζ,ρ}, por un lado; y la funciones {Mx,M y}por otro, La Circulaci´on en el Mar B´altico y su Modelizaci´on 133 necesarias para el modelo de las tensiones que aparecen en las condiciones de contorno: todas ellas son la soluci´on del problema estacionario bidimensional representado en el sistema de ecuaciones anterior (4.77, 4.78, 4.79) cuya salida acoplaremos al sistema de ecuaciones representadas por (4.108) que nos dar´a las velocidades {U, V }del problema tridimensional. Debido a estructura funcional hiperb´olica del sistema de ecuaciones diferenciales a resolver, puede proponerse como soluci´on una funci´on compleja definida a partir de las velocidad horizontales, donde la velocidad Ues la parte real y la velocidad Ves la parte imaginaria D≡U+iV . (4.107) Agrupando la parte inhomog´enea en la funci´on G–dependiente de {Az,p a,ζ,ρ}–e independiente de {U, V },quesonlasinc´ognitas que nos ocupan ahora, queda la ecuaci´on que sigue: ∂2D ∂z2−p2 1D=G(4.108) donde, G=1 ρoAv∂pa ∂x +i∂pa ∂y + +g Av∂ζ ∂x +i∂ζ ∂y+g ρoAv0 z∂ρ ∂x +i∂ρ ∂ydη (4.109) p2 1=f Av ii≡√−1 (4.110) Las condiciones de contorno son las siguientes: ρAv ∂D ∂z =Tspara z =ζ(4.111) ρAv ∂D ∂z =Tbpara z =−H(4.112) 134 CAP´ ITULO 4 donde TsyTbson representaciones complejas de los vectores −→ Ts,−→ Tbcuyas componentes son las tensiones tangenciales que sufre el fluido en la superficie y en el fondo de la columna de l´ıquido respectivamente: −→ Ts≡Ts xz Ts yz−→ Tb≡Tb xz Tb yz(4.113) ydemaneraan´aloga se definen las representaciones complejas: Ts≡Ts xz +iTs yz Tb≡Tb xz +iTb yz (4.114) donde TsyTbya fueron parametrizadas en el apartado 3.3 de este cap´ıtulo: las tensiones en la superficie del oc´eano {Ts xz,T yz}, veremos que son debidas a la acci´on del viento; y las tensiones en el fondo del mar Tb xz,Tb yzson modeladas en funci´on del coeficiente de fricci´on l´ıquido–s´olido y del flujo de masa {Mx,M y}que circula por encima en la columna de l´ıquido. Integrando (4.108), teniendo en cuenta (4.109),(4.111) y (4.112), y la metodolog´ıa de la funci´on de Green, se llega a: D=Ts ρop1Az cosh [p1(H+z)] sinh p1H+ −cosh p1z sinh p1HTb ρop1Az +2B1(x, y, −H)coshp1H+ +B1(x, y, z)ep1z+B2(x, y, z)e−p1z+G1(4.115) donde: B1(x, y, z)= 1 2p1z 0G2(x, y, η)e−p1ηdη (4.116) B2(x, y, z)=−1 2p1z 0G2(x, y, η)ep1ηdη (4.117) G1=1 ρofi∂pa ∂x +i∂pa ∂y +g fi∂ζ ∂x +i∂ζ ∂y(4.118) G2=g ρoAvz 0∂ρ ∂x +i∂ρ ∂ydη (4.119) La Circulaci´on en el Mar B´altico y su Modelizaci´on 135 M´as adelante retomaremos esta soluci´on para terminar de plantear el esquema de corrientes estacionarias tridimensionales que usaremos en esta tesis doctoral. Esto ser´aen el apartado 5 del cap´ıtulo 4, donde tras adaptar las soluciones de Jankowski [80] a nuestra malla de c´alculo, se plantear´a un algoritmo de modulaci´on para extender esta informaci´on (relativa a los 8 vientos principales con Ua= 10m s–1 ) a un escenario meteorol´ogico algo m´as general. 4.ADAPTACI´ ON DE LOS MODELOS PREVIOS A NUESTRA MALLA DE C´ ALCULO Se ha hecho un trabajo arduo en la adaptaci´on de los resultados de estos modelos hidrodin´amicos previos a nuestra malla de c´alculo de resoluci´on horizontal 20 km×20 km. A modo resumido estos han sido los c´alculos y los procesos de adaptaci´onquesehan realizado: (a) Se han practicado interpolaciones lineales eligiendo muestras representativas de los campos de velocidades horizontales tanto al modelo 3D como al 2D. (b) Se ha impuesto la ecuaci´on de continuidad para el c´alculo de la corrientes verticales, con algoritmos que empiezan por los nudos de la capa inferior y suben hasta llegar a su nudo hom´ologo en la capa superficial. Una vez all´ısea˜nade el flujo correspondiente horizontal residual que, respetando la conservaci´on de la materia, m´as se aproxima a los valores de los modelos previos que sirven de fuente. (c) Todos estos c´alculos se han realizado mediante algoritmos MATLAB. Debido a la dificultad de la batimetr´ıa tridimensional, y a la menor resoluci´on de nuestro modelo en relaci´on a los modelos hidrodin´amicos de base, esta sistem´atica de c´alculo programado ha tenido que ser completada en su etapa final por peque˜nos ajustes manuales en puntos de especial dificultad, especialmente del contorno, hasta encajar perfectamente la ecuaci´on de continuidad en el sistema con el m´aximo respeto a las publicaciones hidrodin´amicas fuente de este trabajo. (d) Al hacer las interpolaciones sobre las cartas de corrientes originales aparecen en algunos puntos de la frontera velocidades con componente normal. Estas componentes han sido autom´atica anuladas en los c´alculos posteriores. 136 CAP´ ITULO 4 5.GENERALIDADES DE LOS DOS MODELOS DE CIRCULACI´ ON ADAPTADOS A ESTA TESIS En esta tesis doctoral vamos a trabajar con la hidrodin´amica de dos modelos tridimensionales, a los que llamaremos: modelo HD-3D; modelo HD-C3D. 5.1. Modelo HD–3D Este modelo es la adaptaci´on total a nuestra malla de c´alculo del modelo del SMHI [43] descrito anteriormente. Se trata de un modelo de circulaci´on media, basado en los eventos meteorol´ogicos m´as probable. La labor de adaptaci´on ha sido larga por toda la metodolog´ıa de interpolaci´on y determinaci´on de corrientes verticales manteniendo la inviolabilidad de la ecuaci´on de continuidad. A partir de ahora nos referiremos a ´el como: Modelo HD–3D. A continuaci´on se muestran las cartas de corrientes para las seis capas en que se ha dividido el sistema y que han sido fruto de este proceso de adaptaci´on. Ser´an usadas en nuestro modelo de transporte. El predominio de las corrientes Ekman en las capas m´as superficiales est´an acompa˜nadas por zonas que abandonan esta tendencia dominante para presentar desviaciones e incluso remolinos locales en determinadas regiones. Algunas zonas que se salen de la norma antes descrita son: •El noroeste de la isla de Bornholm, el norte de Polonia, puntos de la zona cercana a la isla de Gotland y la entrada al Golfo de Finlandia. •Tambi´en, en el norte y sur del Mar de Bothnia: circulaci´on anticicl´onica en el sur y cicl´onica en el norte. Y en el norte y sur de la Bah´ıa de Bothnia: circulaci´on anticicl´onica en el sur y cicl´onica en el norte. •Aestohayquea˜nadir, l´ogicamente, el alineamiento costero de las corrientes. Se puede decir que las corrientes geostr´oficas inducidas por los vientos dominantes del oeste y suroeste hacen que en la circulaci´on media del sistema est´em´as dominada por el sentido cicl´onico que por el anticicl´onico. En las capas no superficiales la circulaci´on est´a dominada por la topograf´ıa y la fuerza de Coriolis. Cabe destacar la corriente de retorno hacia el noroeste en las costas de Polonia; y hacia el oeste en Golfo de Finlandia. Son importantes las canalizaciones de corrientes en capas profundas que se producen dos zonas concretas; una en el B´altico Sur y otra en el Mar de Bothnia, y que explican f´acilmente analizando el mapa batim´etrico de profundidades. La informaci´on de estas La Circulaci´on en el Mar B´altico y su Modelizaci´on 137 corrientes medias anuales en las seis capas del B´altico se ha usado como input en diversas publicaciones del doctorando as´ı como en varias partes de esta tesis doctoral. Por otro lado, las corrientes medias tanto estacionales como anuales que aparecen en el modelo del SMHI [43] y que tienen como input los vientos del a˜no meteorol´ogico que se ha expuesto aqu´ı, ser´an muy ´utiles para simular el transporte en condiciones meteorol´ogicas gen´ericas. 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 1400 km km SWEDEN FINLAND BALTIC REPUBLICS 20 cm s−1 Figura 4.14: Corrientes medias anuales en las profundidades 0–10 m obtenidas por interpolaci´on de las corrientes del SMHI [43]. 144 CAP´ ITULO 4 adaptado a nuestra malla de c´alculo de 6 capas, present´andose ahora las corrientes asociadasalascapas1(0–5m)y4(20–40m). Estassolucionesser´an parte de la hidrodin´amica para una importante simulaci´on del cap´ıtulo 7, centrada en la dispersi´on de una mancha radiactiva en el Golfo de Finlandia, varias semanas despu´es del accidente de Chernobyl. Hemos apoyado la hidrodin´amica de esta tesis doctoral fundamentalmente en dos modelos previos (HD – 3D y HD – C3D), que son respectivamente el de corrientes medias anuales [43] y el de ocho vientos [80], y queremos a continuaci´on hacer una peque˜na evaluaci´on de la coherencia mutua de sus soluciones. Dada la diferente naturaleza de los resultados num´ericos disponibles, el ejercicio de comparaci´on de soluciones de ambos modelos ser´am´as cualitativo que cuantitativo: (a) Fundamentos. El lector puede consultar los innumerables puntos de coherencia intr´ınsecos de estos modelos que ya se detallaron en las secciones introductorias (2.1 – 2.10 y 3.1 – 3.6) de este cap´ıtulo. Como consecuencia de todo ello, los resultados ser´an compatibles, como se ver´a a continuaci´on. (b) Resultados. Pueden compararse las figuras 4.21 y 4.22 asociadas al modelo de 8 vientos, con las figuras 4.14 y 4.16, respectivamente. Las primeras corresponden a las soluciones del SW tales que Ua=10ms –1, y las segundas, con las corrientes anuales asociadas al a˜no meteorol´ogico t´ıpico cuya variada gama de vientos, con dominio del SW, ya se ha discuti´o en el apartado 2.10 de este cap´ıtulo. Son razonablemente compatibles los perfiles hidrodin´amicos de ambos estudios: i. Las velocidades en superficie est´an dominadas en ambos casos por las corrientes de Ekman, del orden un dos por ciento de la velocidad del viento en zonas interiores superficiales, con una adaptaci´on adecuada a las condiciones de frontera. ii. La profundidad de inversi´on de las velocidades (posici´on en la columna de l´ıquido donde se anula la velocidad horizontal, por debajo de la cual se inician velocidades opuestas a las de superficie), se sit´ua t´ıpicamente para ambos modelos (3D y C3D) en el entorno de los 20 metros: zinv(3D) ∼zinv(C3D) ∼−20 m .(4.124) La Circulaci´on en el Mar B´altico y su Modelizaci´on 145 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 1400 km km SWEDEN FINLAND BALTIC REPUBLICS 50 cm s−1 Figura 4.21: Corrientes inducidas por vientos del SW con Ua=10ms –1, adaptados de la salida del modelo del IOPAN [80] a la capa 1 (0 – 5 m) de nuestra malla de c´alculo. 146 CAP´ ITULO 4 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 1400 km km SWEDEN FINLAND BALTIC REPUBLICS 20 cm s−1 Figura 4.22: Corrientes inducidas por vientos del SW con Ua=10ms –1, adaptados de la salida del modelo del IOPAN [80] a la capa 4 (20 – 40 m) de nuestra malla de c´alculo. Recordemos, que en las soluciones de U0(x, y, z), V0(x, y, z)yD0(x, y, z) provenientes del modelo del IOPAN [80], no se tienen en cuenta las corrientes por densidad, aproxima- La Circulaci´on en el Mar B´altico y su Modelizaci´on 147 ci´on que vamos a procurar matizar y delimitar a continuaci´on. No son por lo general importantes en este sistema las corrientes asociadas a los gradientes medios horizontales de densidad. Puede comprobarse que si hacemos una estimaci´on en la f´ormula de G2introduciendo el gradiente horizontal medio de salinidad en el Mar B´altico, que es del orden de |∇S|∼ 10 PSU 1000 km =⇒|∇ρ|∼1kgm –3 100 km (4.125) obtenemos velocidades peque˜nas (∼1cms –1) tanto comparadas con las velocidades t´ıpicas asociadas al viento (∼10 cm s–1) como con las geostr´oficas, que en las zonas de mayor desarrollo compiten (∼10 cm s–1)conestas´ultimas. Existen excepciones [81], como la zona suroeste del mar B´altico, donde los gradientes horizontales de salinidad pueden ser importantes en las corrientes de fondo, factor que s´ıestenidoencuentaenelmodelo de Simons [170] y Funkquist [43], en que se basan las corrientes medias anuales de la otra aproximaci´on hidrodin´amica de esta tesis doctoral aplicada al mar B´altico. En otro trabajo de Jankowski [81] – a diferencia del trabajo de los 8 vientos [80] que trabaja con densidad constante – se intenta ver la influencia de las corrientes de densidad en relaci´on con las corrientes Ekman para diferentes velocidades de viento. Algunas conclusiones generales de este trabajo [81], como primera aproximaci´on para las aplicaciones de esta tesis doctoral, son las siguientes: (a) Las corrientes Ekman son dominantes en las capas superficiales del mar B´altico. (b) Respecto a las capas profundas (donde en presencia de tensiones insuficientes para la resuspensi´on de materia en suspensi´on, transcurren los procesos de acumulaci´on hacia los sedimentos), el mar B´altico profundo, puede ser dividido en dos zonas: i. Zona profunda A: gobernada por corrientes Ekman. B´altico Estricto norte, Golfo de Finlandia y Mar de Botnia. ii. Zona profunda B: gobernada por corrientes por salinidad. B´altico Estricto centro y sur, y Bah´ıa de Botnia. El acercamiento de las soluciones de ambos modelos (comparaci´on de figuras 4.14 y 4.16 con 4.21 y 4.22) es a´un mayor en las partes del mar B´altico con zonas profundas tipo 148 CAP´ ITULO 4 A(B´altico Estricto norte, Golfo de Finlandia, Mar de Botnia), donde la influencia de las corrientes por densidad es menor. Se puede decir que el modelo HD–C3D ser´a una buena alternativa tanto para el modelado del transporte hidrodin´amico en las capas superficiales del sistema, como en la zona profunda A. En la zona profunda B puede trabajarse con el modelo HD–3D de corrientes medias anuales del SMHI que tiene en cuenta el transporte por gradientes de densidad. A continuaci´on plantearemos una modulaci´on del modelo para vientos con velocidades distintas a Ua=10ms –1 en las ocho direcciones principales de la rosa de los vientos. 5.2.1. Algoritmo de modulaci´on seg´un la velocidad del viento Vamos a descomponer la expresi´on de D(4.121) en diferentes t´erminos. Los dos primeros son los correspondientes a la capa Ekman superficial y de fondo respectivamente: DE,s ≡Ts ρop1Az cosh [p1(H+z)] sinh p1H≡Ts ρop1Az fs(z) (4.126) DE,b ≡− Tb ρop1Az cosh p1z sinh p1H≡− Tb ρop1Az fb(z) (4.127) donde se han definido los t´erminos fs(z)yfb(z) que dependen de la profundidad. Integrando en la columna de l´ıquido se obtienen los valores medios siguientes: fs(z)= 1 p1H=fb(z).(4.128) Si definimos la longitud dE≡1/|p1|se comprueba que dEes la zona de influencia del viento en las capas de Ekman, pudiendo ser usada para caracterizar sus espesores de influencia. Tras hacer unas estimaciones, mediante el modelo de Avpropuesto en (4.82), introduciendo diferentes valores t´ıpicos de profundidad y velocidad de viento, se obtiene que en la mayor´ıa de los casos, el valor de dEest´a en rango 10 −20 metros. Los t´erminos restantes, que dependen de la presi´on, est´an asociados respectivamente alapresi´on atmosf´ericayalafunci´on profundidad, siendo notadas las velocidades con sub´ındices alusivos al origen atmosf´erico y geostr´ofico: La Circulaci´on en el Mar B´altico y su Modelizaci´on 149 Datm =1 ρofi∂pa ∂x +i∂pa ∂y (4.129) Dgeo =g fi∂ζ ∂x +i∂ζ ∂y.(4.130) Para el algoritmo de modulaci´on nos interesa hacer un breve an´alisis de la dependencia de algunas velocidades parciales estacionarias con la velocidad del viento Ua. Si trabajamos con el valor medio de CD Ts∝U2 a(4.131) Seg´un la teor´ıa de Fenselbaum, expresada en modelos de la literatura [155, 80], y seg´un las condiciones, Azpuede ser modelada proporcionalmente tanto a Ua,comoasu cuadrado U2 a Az∝Ua,obien ,A z∝U2 a(4.132) en funci´on de ciertas condiciones de la profundidad y velocidad de viento, estando ambas aproximaciones ligadas respectivamente a zonas de poca o mucha profundidad. En esta discusi´on trabajaremos con la expresi´on Az∝Uα a,1≤α≤2 (4.133) Por definici´on se cumple que p1∝A−1/2 z∝U−α/2 a(4.134) Si tomamos valor medio en DE,s y evaluamos su dependencia respecto del viento llegamos a DE,s =Ts ρop1Az 1 p1H∝Ts p2 1Az∝Ts∝U2 a(4.135) Evaluando la velocidad de Ekman en la superficie (z=0) 150 CAP´ ITULO 4 DE,s(0) ∝TsA−α/2 z∝U2−α/2 a(4.136) La relaci´on entre la velocidad del viento y el gradiente de presi´on atmosf´erica puede ser modelado [45] as´ı Ua=CrBT∇npa(4.137) donde Cres un factor de reducci´on que tiene en cuenta el debilitamiento del viento por fricci´on sobre la superficie del mar y BTun coeficiente de estratificaci´on para la interfase mar–atm´osfera. Unos valores de trabajo, seg´un [45], son Cr=0.7yBT=4.7,donde el valor del gradiente vendr´a dado en las unidades Pa/km . Por tanto ∇xpa∝U1 a,∇ypa∝U1 a,(4.138) lo que nos lleva a suponer que Dgeo ∝∇ npa∝U1 a.(4.139) Hagamos ahora la discusi´on sobre DE,b yDgeo .De manera an´aloga a la estimaci´on de DE,s se llega a DE,b ∝Tb∝−→ M2(4.140) La velocidad geostr´ofica no es f´acilmente modelable en funci´on de la velocidad del viento. Peri´a˜nez [147] estim´o una relaci´on l´ımite ideal bajo condiciones de poca influencia de la tensi´on de fondo –situaci´on t´ıpica en zonas relativamente profundas– que es la que sigue Dgeo ∝U2 a(4.141) siendo por tanto razonable trabajar con la relaci´on Dgeo ∝Uβ a,1≤β≤2.(4.142) La Circulaci´on en el Mar B´altico y su Modelizaci´on 151 Teniendo en cuenta toda esta discusi´on vamos a trabajar con una relaci´on global para la velocidad total que ser´a D∝Uγ a,1≤γ≤2.(4.143) As´ı, una primera aproximaci´on para un algoritmo global de modulaci´on ser´ıa el siguiente D(x, y, z)=D0(x, y, z)Ua 10γ ,1≤γ≤2.(4.144) Para la simulaciones, elegiremos para este par´ametro, un valor de trabajo intermedio, que por defecto ser´aγ=3/2. Respecto a la modulaci´on de la tensi´on de fondo Tb, llamaremos Tb 0alatensi´on producida por cada uno de los 8 vientos del modelo de referencia [80]. Si suponemos que dada una columna de l´ıquido situada en el punto P(x, y) se produce una proporcionalidad entre Ts(P,UA)yTb(P,UA) Tb(P,UA)∝Ts(P,UA)∝U2 a(4.145) se deduce que Tb(x, y)=Tb 0(x, y, )Ua 102 (4.146) si bien dejaremos una propuesta m´as abierta Tb(x, y)=Tb 0(x, y, )Ua 10λT ,1≤λT≤2.(4.147) aunque trabajemos con λT= 2 como valor por defecto. Este algoritmo de modulaci´on es ampliable a otro m´as detallado que consiste en establecer estas relaciones en funci´on de la profundidad: ser´ıa un algoritmo de modulaci´on diferencial que se completar´ıa finalmente con el c´alculo de las corrientes verticales a partir de la ecuaci´on de continuidad empezando por la capa inferior y siguiendo hasta acabar en capa superficial. 152 CAP´ ITULO 4 L´ogicamente, todos los par´ametros que aparecen en el modelo, podr´an hacerse variar para ver la sensibilidad del modelo a los mismos, y buscar tambi´en el mejor de los acuerdos entre los resultados del modelo y la experiencia. 5.2.2. Correcciones por el efecto atenuante de la picnoclina La energ´ıa que absorbe el oc´eano por su contacto con la atm´osfera encuentra cierta barrera en la picnoclina permanente de cara a transmitirse [127, 86] a las capas m´as profundas. Esta es una limitaci´on del modelo del IOPAN [80], y por tanto del HD– C3D, para la soluci´on que ofrece en las capas profundas bajo la picnoclina permanente del mar B´altico. Este efecto atenuante, llevado al extremo, invalidar´ıa la soluci´on en la parte m´as profunda de estas columnas de l´ıquido. Podemos plantearnos dos correcciones aesteproblema. Unaser´ıa trabajar con un coeficiente de atenuaci´on a las corrientes de Ekman bent´onicas (de fondo). Hay autores que exponen que estas ´ultimas pueden despreciarse [147] si la picnoclina permanente es suficientemente marcada, por lo que una segunda opci´on, para estas capas profundas bajo picnoclina muy desarrollada, podr´ıa ser, no trabajar con las corrientes Ekman de fondo, sino con la soluci´on de circulaci´on media se propone el SMHI [43] para capas profundas, que tiene en cuenta el factor de la estratificaci´on vertical. La obtenci´on de las corrientes verticales se har´ıa al final, en base a la ecuaci´on de continuidad. La limitaci´on del modelo del IOPAN [80] sobre las aguas que est´an bajo una picnoclina permanente afecta tambi´en al valor que este modelo ofrece sobre las tensiones de fondo bajo las columnas de l´ıquido afectadas. Esto afectar´ıa te´oricamente a los valores de la tensi´on de fondo para nuestro modelo de transporte de materia en suspensi´on del cap´ıtulo 8, que se nutren de la modulaci´on del modelo del IOPAN para obtener Tb,pero no tendr´a consecuencias pr´acticas, ya que los fen´omenos altamente dependientes de Tb, como son los episodios de resuspensi´on, son de car´acter despreciable [86] bajo la capa picnoclina. Dicho de otro modo, s´ı trabajaremos con el valor que el modelo ofrecer´apara Tben la base de capa Ekman bent´onica asociada a columnas de l´ıquido que no atraviesen una capa picnoclina permanente, aquellas en las que precisamente s´ısedanlosfen´omenos de resuspensi´on, importantes para din´amica de la materia en suspensi´on as´ı como para el transporte de radion´uclidos no conservativos. Cap´ıtulo 5 La Difusi´on en el Mar B´altico y su Modelizaci´on 1.LA ENERG´ IA DE LAS CORRIENTES EN EL MAR B´ ALTICO 1.1. Introducci´on Ya se vio en el cap´ıtulo 3 que las diferentes formas de modelar un fluido oce´anico dependen doblemente de la escala del problema que estemos estudiando y de la resoluci´on espacial y temporal de nuestro modelo num´erico. Remolinos de diferentes escalas se superponen en el espacio y el tiempo contribuyendo a la dispersi´on de las sustancias que forman parte del fluido. El filtrado matem´atico ayuda a separar y estudiar los distintos tipos de oscilaciones que entran en juego en el problema de la difusi´on, por lo aqu´ı seguiremos un esquema que empezar´aconladifusi´on a menor escala para ir posteriormente estudiando las escalas mayores. Nuestro modelo de dispersi´on trabajar´a con una circulaci´on media mensual o anual, cuya difusi´on dominante es la de mesoescala. A´un as´ı, interesa ver tambi´en aqu´ı otras fuentes de difusi´on de escala menor, como la diaria, aunque sus coeficientes de intercambio est´en impl´ıcitamente incluidos en la difusi´on de escala mayor, ya que as´ı entenderemos mejor el problema que nos ocupa y puede sernos de utilidad en futuros modelos de una resoluci´on mayor. Para poder desarrollar esta parte del trabajo, contribuir´a positivamente el an´alisis de los espectros de corrientes y la distribuci´on de energ´ıa en funci´on de los diferentes fen´omenos f´ısicos que entran en juego en las distintas escalas espaciales y temporales. 153 256 CAP´ ITULO 7 150 m2s-1 ≤Kx≤1000 m2s-1 .(7.19) El modelo lagrangiano usa un n´umero de part´ıculas trazadoras proporcional a la salinidad, dada en unidades PSU. La aproximaci´on 1D parece razonable para este sistema, relativamente bien canalizado, debido a las diferentes intervenciones (canalizaciones, cortas y dragados) a que ha sido sometido en el ´ultimo siglo. La discretizaci´on del modelo se hace con Δx=1km,with Δt=1h.Comoprimera aproximaci´on, la circulaci´on residual ha sido estimada asumiendo un flujo medio en cada secci´on igual al aporte de agua dulce en el contorno con la presa. Como condiciones inicales, un n´umero de part´ıculas N(x) = 100 S×Δxfue asignado en cada celda, siendo Sla salinidad. Como condici´on de contorno en la desembocadura del estuario, se ha incluido un compartimento que reproduce el mar abierto mediante una condici´on de mezcla r´apida y homog´enea. Tras imponer las condiciones descritas, y tras un periodo transitorio de unos pocos d´ıas, se alcanza el estado estacionario. Un acuerdo razonable entre el modelo y las observaciones se alcanza cuando se introduce la siguiente parametrizaci´on del coeficiente efectivo de difusi´on horizontal: Kx(x)=K01−γx lE2 (7.20) donde K0es el correspondiente valor del coeficiente cuando x=0,γes un par´ametro libre que ser´a calibrado, y lEes la longitud del estuario. La siguiente figura (7.19) muestra los resultados del estado estacionario de la distribuci´on de salinidad (en unidades PSU) obtenida para un aporte de agua dulce Q=25m 3s–1 (t´ıpico durante peridos secos), y se compara con valores experimentales [159]. An´alogamente, la otra figura (7.20) muestra el correspondiente resultado para un aporte dulce igual a Q=60m 3s–1,t´ıpico durante los periodos lluviosos. Estas curvas de salinidad son una primera aproximaci´on a la problem´atica de la intrusi´on salina en el estuario, conocido como el problema del ”tap´on de sal” por los habitantes de los n´ucleos Modelo Lagrangiano de Transporte Conservativo 257 de poblaci´on lim´ıtrofes con el estuario, conocedores de unas condiciones medioambientales tan pr´oximas a las marinas como prolongado sea el avance del ”tap´on” r´ıo arriba. 0 20 40 60 80 100 0 5 10 15 20 25 30 35 km Salinity(PSU) Model Experiment Figura 7.19: Comparaci´on entre los valores de la salinidad experimentales (promediados sobre ciclos semidiurnos) y los calculados con el modelo, como funci´on de la distancia a la desembocadura. El aporte dulce fue igual a Q=25m 3s–1 (t´ıpico durante periodos secos), K0= 300 m2s–1,γ=3/4ylE= 110 km. Hay puntos muy sensibles (ver figuras 7.19 y 7.20) a estas variaciones de la curva, como el punto x= 50 km , donde la salinidad en ´epoca seca es unas cinco veces superior a la de los periodos lluviosos. El acuerdo razonable entre las predicciones del modelo y las medidas le da apoyo a esta versi´on preliminar de nuestro modelo 1D de dispersi´on lagrangiana para el estuario del r´ıo Guadalquivir. Futuras mejoras ser´an incluidas en este modelo, como las oscilaciones de marea, coeficientes de difusi´on dependientes del tiempo y una m´as detallada descripci´on de la hidrodin´amica del estuario. 258 CAP´ ITULO 7 0 20 40 60 80 100 0 5 10 15 20 25 30 35 km Salinity(PSU) Model Experiment Figura 7.20: Iguales condiciones que en la figura anterior, para un aporte dulce igual a Q=60m 3s-1 (t´ıpico durante periodos lluviosos), K0= 600 m2s-1,γ=1/2ylE= 110 km . Un conocimiento m´as detallado del estuario del r´ıo Guadalquivir viene siendo demandado en los ´ultimos a˜nos por razones socioecon´omicas y ecol´ogicas: labores de dragado peri´odicas para mantenimiento de la navegabilidad, acuicultura, agricultura intensiva de arroz, turismo, y el impacto ecol´ogico en el Parque Nacional de Do˜nana, entre otros. Por este motivo, un modelo de transporte m´as avanzado y ramificado ser´a necesario en el futuro. Cap´ıtulo 8 Materia en Suspensi´on y su Modelizaci´on 1.ASPECTOS GENERALES DE LA MATERIA EN SUSPENSI´ ON EN EL OC´ EANO En el desarrollo de esta tesis doctoral estamos interesados en el estudio y modelizaci´on del transporte de part´ıculas en general y de radion´uclidos en particular cuando se encuentran inmersos en el oc´eano. Nuestro inter´es se centra en radion´uclidos que se pueden encontrar tanto en fase disuelta como en fase s´olida, bien agregados a la materia en suspensi´on o formando parte del sedimento marino. Los radion´uclidos que se encuentran s´oloenlafasedisueltase denominan conservativos, mientras que los que se presentan ambas fases son conocidos como radion´uclidos no conservativos. Tanto la fase disuelta como la fase s´olida suspendida participan de todos los procesos hidrodin´amicos que hemos expuesto hasta ahora. Los radion´uclidos disueltos forman parte del movimiento del fluido como si de mol´eculas de agua se tratara, mientras que los radion´uclidos en la fase s´olida requieren un tratamiento aparte ya que previamente hay que estudiar la din´amica de las part´ıculas s´olidas a las que van adheridos. Como escenario de partida, habr´a que introducir el modelo hidrodin´amico previo -ya que las part´ıculas que estamos estudiando son pasivas y no intervienen en la din´amica del fluido que las soporta-, la informaci´on meteorol´ogica necesaria, los datos sobre las fuentes de materia en suspensi´on y la din´amica espec´ıfica de la fase s´olida en el seno del fluido. 259 260 CAP´ ITULO 8 A partir de ah´ı, levantaremos un sencillo y suficiente por ahora Modelo de la Materia en Suspensi´on –o Modelo de Part´ıculas Suspendidas–, que nos dar´a la concentraci´on de part´ıculas s´olidas a lo largo del tiempo en el sistema, a las que ir´a agregada la parte no disuelta de los radion´uclidos que estamos estudiando. En el Modelo de la Materia en Suspensi´on se habr´an de tener en cuenta dos procesos fundamentales de las part´ıculas s´olidas en el oc´eano: la deposici´on y la resuspensi´on. Posteriormente habr´a que tener en cuenta la posibilidad de que radion´uclidos de la fase disuelta pasen a la fase s´olida y viceversa. La primera aproximaci´on que podemos hacer a este problema es suponer que se encuentran en equilibrio, condici´onquesealcanza al cabo de un tiempo dependiente de cada radion´uclido, pero que en nuestro caso ser´a corto comparado con las escalas temporales que vamos a manejar para nuestro estudio del sistema. El intervalo transitorio cin´etico o de no equilibrio, ser´atenidoencuentaen futuras ampliaciones del modelo cuando se haga un escalado m´as fino del problema. Por convenio, se llama materia en suspensi´on al material retenido en un filtro para part´ıculasconuntama˜no de poro entre 0.4μmy0.5μm . El material que no es retenido en este sistema y por tanto con un di´ametro φ<0.4μm es considerado como materia disuelta. Se considera materia suspendida fina [3, 29] a aquella cuyo di´ametro φ<62.5μm. Laspart´ıculas con di´ametro superior se depositan mucho m´as r´apidamente incorpor´andose a sedimentos de una textura mucho m´as gruesa que los que van a protagonizar nuestro estudio. La raz´on fundamental de esto ´ultimo es que las part´ıculas finas son mucho m´as eficaces que las gruesas a la hora de transportar radion´uclidos, por la sencilla explicaci´on geom´etrica de que una masa mde part´ıculas finas presenta una superficie de adherencia potencial para los radion´uclidos mucho mayor que la misma cantidad mde part´ıculas gruesas. Basta considerar el siguiente sencillo experimento mental de romper una part´ıculagruesahomog´ enea en diferentes finas: sup´ongase por sencillez que es c´ubica, de lado ay que se le practican Ncortes paralelos a cada uno de los lados del cubo, apareciendo N3cubitos de lado a/N. Si llamamos SAySBa las superficies de ambos sistemas de igual masa mllegamos a: SA=6a2SB=N3×6(a/N)2=6a2N(8.1) La idea es que la capacidad de absorber radion´uclidos por una part´ıcula aumenta Materia en Suspensi´on y su Modelizaci´on 261 cuando decrece [4, 5, 29] su di´ametro caracter´ıstico φ,raz´on por la que nos interesamos ´unicamente en la part´ıculas suspendidas finas. La superficie SyelvolumenVde una part´ıcula s´olida est´an ligados a trav´es de una constante kydeldi´ametro φ, V=kφS (8.2) aunque, seg´un las investigaciones de Abril [7, 8, 9] y otros autores, la relaci´on puede ser m´as compleja, ya que la superficie es rugosa, con poros y bordes libres, que pueden aumentar la superficie efectiva de adsorci´on. Lamasadelapart´ıcula est´a ligada con la densidad media ρ, m=ρV (8.3) luego si calculamos la relaci´on entre la superficie S´util por unidad de masa para absorber radion´uclidos, y la denominamos σm,llegamos a σm≡S m=1 kρφ (8.4) donde se nota su dependencia inversamente proporcional con el di´ametro de las part´ıculas s´olidas φ. En el caso de part´ıculas esf´ericas, k=6,y un rango frecuente [204] para ρen algunos de nuestros problemas es ρ∼2−3gcm –3 . Aunque la materia en suspensi´on est´a formada por part´ıculas finas (φ<62.5μm) y gruesas, dada la importancia de las primeras en los problemas que vamos a abordar, se sobrentender´a a partir de ahora y mientras no se diga lo contrario que la materia en suspensi´on son las part´ıculas finas suspendidas. Se presentar´a a continuaci´on la tabla 8.1, donde se muestran el nombre y los tama˜nos de grano del amplio abanico de s´olidos que pueden formar parte de un sedimento. Para esta tabla, y por su utilidad en futuros an´alisis, se va introducir la escala logar´ıtmica de Wentworth [107], que toma como referencia el grano grueso de arena de un mil´ımetro de di´ametro, y que coloca a los granos m´as peque˜nos en la parte positiva de la escala como premio a su car´acter din´amico: 262 CAP´ ITULO 8 Φ=log 2 φ0 φ⇐⇒ φ=φ02−Φ,φ 0=1mm (8.5) Unida a esta escala usaremos la norma DIN 4022 para la nomenclatura y segmento de los diferentes tipos de s´olidos: arcilla, limo, arena, gravas, y piedras. Se dibuja una l´ınea de separaci´on para distinguir las part´ıculas finas (φ<62.5μm : limo y arcilla) de las dem´as. TABLA 8.1: Clasificaci´on de s´olidos en suspensi´on en funci´on de su di´ametro Tipo de s´olido φΦ Piedras [64 ,∞]mm [−∞,−6] Gravas [2 ,64]mm [−6,−1] Arena [1/16 ,2] mm = [0.0625 ,2]mm [−1,4] Limo [1/512 ,1/16] mm ∼ =[2 ,62.5]μm[4 ,9] Arcilla [0 ,1/512]mm ∼ =[0 ,2]μm[9 ,∞] 2.DIN´ AMICA GENERAL DE PART´ ICULAS SUSPENDIDAS EN EL OC´ EANO Ya se ha explicado que las part´ıculas suspendidas participan pasivamente del movimiento del fluido cuyo modelo y soluci´on conocemos de los bloques previos. Lo que ahora toca es poner sobre la mesa las condiciones bajo las cuales las part´ıculas rompen la disciplina del fluido y se mueven respecto a ´el apareciendo un problema s´olido–fluido que es representado principalmente por los siguientes procesos: (a) Deposici´on hacia el fondo marino. (b) Resuspensi´on desde el fondo marino. Una expresi´on muy usada en la literatura [75, 117, 119, 208, 150] para el flujo por deposici´on es la siguiente: DEP = wsms1−Tb Tcd (8.6) donde wses la velocidad de deposici´on de las part´ıculas suspendidas, mses la funci´on concentraci´on de part´ıculas suspendidas (materia en suspensi´on), Tbes la tensi´on en el Materia en Suspensi´on y su Modelizaci´on 263 fondo marino, y Tcd la tensi´on cr´ıtica para la deposici´on de las part´ıculas suspendidas. La tensi´on Tcd es un valor por encima del cual el proceso de deposici´on –o sedimentaci´on– no tiene lugar. La velocidad de deposici´on ser´a calculada a partir de la ley de Stokes (Skreen, Irlanda, 1819 – Cambridge, R.U., 1903) para el movimiento de ca´ıda –o sedimentaci´on– de un s´olido esf´erico en el seno de un fluido ws=ρ−ρw ρw gφ2 18ν(8.7) donde {ρ, φ}son respectivamente la densidad y el di´ametro de la part´ıcula suspendida, as´ıcomo{ρw,ν}son la densidad y la viscosidad cinem´atica del agua. A efectos de c´alculos tomaremos la viscosidad cinem´atica a temperatura ambiente [21] como ν≃1.10−6m2s–1 . Y una f´ormula estimativa para wsen la regi´on que nos interesa es: ws(m s–1)∼φ2(μm2)×10−6(8.8) que evaluada para tres di´ametros concretos en el intervalo de inter´es nos da: ws(1μm) ∼1×10−6ms –1 ,w s(2μm) ∼4×10−6ms –1 ,w s(20μm) ∼4×10−4ms –1 .(8.9) Una expresi´on muy usada en la literatura [117, 118, 150, 75, 129] para el flujo por resuspensi´on –o erosi´on–eslasiguiente: ER = EfTb Tce −1(8.10) donde Ees la constante de erosi´on sobre el sedimento, fes la fracci´on de part´ıculas finas en la superficie del sedimento, Tbes la tensi´on en el fondo marino, y Tce es la tensi´on cr´ıtica para la erosi´on del sedimento. La tensi´on Tce es un valor por debajo del cual el proceso de resuspensi´on -o erosi´onno tiene lugar. En la literatura [143] pueden encontrarse valores para esta tensi´on en el siguiente intervalo Tce ∈[0.1,1.5] N m–2 .(8.11) 264 CAP´ ITULO 8 Para sedimentos no cohesivos [143] se trabaja admitiendo iguales las dos tensiones cr´ıticas, bajo la hip´otesis: Tcd =Tce (8.12) Y la constante de erosi´on Ealcanza valores [152] como el que sigue: E=0.04 g m2s–1 .(8.13) En el art´ıculo publicado por Kuhrts [102] sobre una parte del B´altico se trabaja con los siguientes valores: ws=4×10−4cm s–1 ,T cd =Tce =0.4Nm –2 .(8.14) En esta parte, nuestro modelo va a trabajar con las siguientes aproximaciones: (a) Las fuentes de materia en suspensi´on se han separado en org´anica e inorg´anica. (b) Para la primera se us´o una velocidad de sedimentaci´on ws=4×10−4cm s–1 y unas tensiones Tcd =Tce =0.4Nm –2 , tomadas de un estudio del transporte de materia en suspensi´on [102] en el sur del Mar B´altico donde domina la producci´on primaria, de lo que, usando la f´ormula de Stokes, se deducen densidades para los s´olidos org´anicos ligeramente (del orden de un dos por ciento) superiores a la del agua. (c) La materia inorg´anica se ha separado en tres tama˜nos de grano siguiendo un esquema similar al publicado en [143], con di´ametros {φ1=2μm,φ 2=8μm,φ 3=32μm}y tensiones cr´ıticas {0.3,0.6,1.2}Nm –2 respectivamente. 3.C´ ALCULO DE LA TENSI ´ ON DE FONDO −→ TbAPARTIRDELMODELO HIDRODIN´ AMICO 3.1 Modelo conceptual de −→ Tb Nuestro modelo hidrodin´amico, basado a su vez en otros modelos previos de los institutos IOPAN y SMHI respectivamente [80, 43], plantea una parametrizaci´on [101, 155, 168] de las tensiones Tb x,Tb ysobre el fondo marino an´aloga a la que se hace para las tensiones en la superficie: el vector de tensiones debido al fondo oce´anico −→ Tblleva la direcci´on del vector transporte de masa −→ Mysum´odulo Tbes proporcional al cuadrado del Materia en Suspensi´on y su Modelizaci´on 265 vector −→ My a la inversa del cuadrado de la profundidad total H+ζ.Aqu´ı se muestra su dependencia respecto a −→ M,sum´odulo y sus componentes, la profundidad y el coeficiente de fricci´on r: −→ Tbub≡−→ M −→ M,T b∝M21 (H+ζ)2⇒Tb=rM21 (H+ζ)2 ⇒Tb x=rM (H+ζ)2Mx,T b y=rM (H+ζ)2My, (8.15) donde ubes el unitario en la direcci´on del vector −→ Mcuyo m´odulo se puede expresar en funci´on de U(x, y), V (x, y). Recordemos que −→ Mes el vector transporte de masa −→ Mdefinido como la cantidad de masa por unidad de tiempo que atraviesa un rect´angulo imaginario de base la unidad y de altura H+ζcolocada perpendicularmente a al vector −→ M=[MxMy]: Mx=ζ −HρUdz My=ζ −HρV dz (8.16) siendo {U, V }las componentes horizontales del campo de velocidades del modelo del IOPAN [80]. En los c´alculos se tom´o [155, 168] un valor para el coeficiente de fricci´on r=2.5×10−3(Sistema Internacional). Es de vital importancia recordar que del modelo del IOPAN [80] se obtienen unas gr´aficas (cap´ıtulo 4) sobre las velocidades horizontales promediadas a lo largo de la columna de l´ıquido U0(x, y), V 0(x, y)obtenidas para las 8 direcciones principales y una velocidad del viento tal que Ua=10ms –1 . Lo importante es que estas velocidades U(x, y), V (x, y)est´an ´ıntimamente ligadas, como ahora veremos, a las componentes del vector −→ M, funciones principales de las que depende nuestro objetivo ahora, Tb: Mx=ζ −HρUdz =ζH,ρ∼ =cte =ρHU, (8.17) My=ζ −HρV dz ζH,ρ∼ =cte =ρHV, (8.18) y sustituyendo en Tb: Tb(x, y)=rρ2U2+V2(x, y)r∼3.10−3(S.I.),ρ∼103kg m–3.(8.19) 272 CAP´ ITULO 8 p0,1≡p(x0,x 1)=x1 x0 pW(x)dx =e−x2 0−e−x2 1,con p(0,∞)=1.(8.26) Generemos sucesos de acuerdo con esa distribuci´on. Si tenemos un total de Nsucesos y llamamos N0,1al n´umero de ellos en el intervalo [x0,x 1], se cumplir´aque: N0,1=Ne−x2 0−e−x2 1(8.27) Si asignamos el primer suceso a x0= 0 y el segundo a x1: N0,1=1=N1−e−x2 1=⇒x1=ln N N−1(8.28) Si asignamos el tercero a x2y generalizamos: x2=ln N N−2,x n=ln N N−n,n=0, ..., N −1.(8.29) Para obtener xnseg´un la distribuci´on de probabilidad pWbastar´a generar n´umeros aleatorios enteros nentre 0 y N−1. Si usamos la funci´on de RANWde MATLAB que los genera en el intervalo [0,1] y redondeamos al entero m´as cercano, tendremos nyxn n=(N−1) RANW:xn=ln N N−(N−1) RANW .(8.30) Para pasar al continuo, hacemos tender el n´umero de sucesos Na infinito, con lo que N−1≈Nylaexpresi´on se simplifica: x=ln 1 1−RANW ,(8.31) donde RANWson n´umeros aleatorios entre 0 y 1 que se generan en un proceso de Monte Carlo. Si Uaes el valor medio de la distribuci´on de Weibull de las velocidades, y llamamos Uaa las velocidades que queremos generar, podremos dimensionar el esquema anterior x x=Ua Uax=√π 2(8.32) por lo que el algoritmo de generaci´on de velocidades Uaser´a Materia en Suspensi´on y su Modelizaci´on 273 Ua=2Ua √πln 1 1−RANW ,(8.33) donde Uavan a ser las t´ıpicas velocidades estacionales medias {8,6,5,7}ms –1 . N´otese que estamos impl´ıcitamente asumiendo iguales tiempos de permanencia para los diferentes vientos de la distribuci´on. Una forma m´as directa de llegar al resultado anterior es trabajar con una funci´on generatriz G(x): la funci´on probabilidad entre 0 y un punto gen´erico xresultante de integrar la funci´on densidad de probabilidad, G(x)=1−e−x2, cuyo recorrido va desde 0 hasta 1, como el generador de n´umeros aleatorios RANW. Si generamos n´umeros aleatorios entre 0 y 1 vamos obteniendo rectas horizontales cuyos cortes con G(x) originan los puntos cuyas abscisas nos van a dar las velocidades de vientos deseadas. As´ı funciona nuestro algoritmo estoc´astico para generaci´on de velocidades de viento weibullianas: G(x)=1−e−x2=RAN W=⇒x=ln 1 1−RANW ,(8.34) (b) Generaci´on de ua: Rosas de vientos. Una vez generado Uapara un determinado d´ıa del a˜no, pueden darse dos casos: i. Ua<6ms –1 . En este caso generaremos uaaleatoriamente seg´un lo visto en el cap´ıtulo 7. Dividiremos el segmento 0–1 en 8 trozos iguales, uno para cada tipo de viento y generaremos n´umeros aleatorios entre 0 y 1. ii. Ua≥6ms –1 . En esta situaci´on utilizaremos las tablas 7.2 y 7.3 sobre probabilidades direccionales cuya representaciones geom´etricas constituyen las rosas de los vientos. Dividiremos el segmento 0–1 en 8 trozos proporcionales a la probabilidad de cada viento y finalmente, generaremos n´umeros aleatorios entre 0 y 1 para obtener el unitario ua. (c) Generaci´on de tiempos de permanencia tw: Espectro de energ´ıa. De los espectros de energ´ıa del cap´ıtulo 5, dedicado a la difusi´on y la energ´ıa en el Mar B´altico, se desprende que los picos energ´eticos asociados a los vientos se encuentran entre los periodos T=3d´ıas y T=4d´ıas, con ca´ıda hacia ambos lados del eje horizontal. Por ello se propone a continuaci´on una funci´on simple de 274 CAP´ ITULO 8 probabilidades cuyo perfil sea an´alogo al perfil de energ´ıas. Se trata de una funci´on discreta de probabilidades con simetr´ıa en el intervalo T∈[1,6] d´ıas, concentrada en los valores enteros del intervalo, y nula en el resto de la recta real positiva pw(Ti)= 1 12Ti,T i=1,2,3; pw(Ti)= 1 12 (7 −Ti),T i=4,5,6. pw(Ti)=0 ,en el resto de casos . (8.35) Dividiremos el segmento 0–1 en 6 trozos proporcionales a la probabilidad de cada periodo y finalmente, generaremos n´umeros aleatorios entre 0 y 1 para obtener tw. 4.MODELO DE TRANSPORTE DE LAMATERIAENSUSPENSI ´ ON 4.1 Fuentes de materia en suspensi´on en el Mar B´altico Ya que no pod´ıa ser de otro modo, las fuentes que vamos a considerar en el modelo de materia en suspensi´on, son tambi´en fuentes principales [29] en la inmensa mayor´ıa de los ecosistemas del planeta: (a) Los r´ıos. (b) La producci´on primaria (materia org´anica). Como primera aproximaci´on y omitiendo muchos detalles sobre determinados puntos del sistema e instantes concretos a lo largo del a˜no, se puede adelantar que ambas contribuciones van a estar normalmente en el mismo orden de magnitud, y que la aportaci´on org´anica global ser´a ligeramente mayor. 4.1.1 Materia en suspensi´on por r´ıos Seg´un estudios publicados sobre esta materia [68] el aporte de los r´ıos al Mar B´altico en forma de materia en suspensi´on es del orden de 1010 kg anuales . El volumen anual de agua suministrado [209, 34, 35, 125] por los r´ıos al sistema es de unos 15000 m3s–1 ,que en otras unidades ser´ıan unos 500 km3a–1 ,loquesuponeun 2.5 % del volumen del sistema. El aporte al Mar B´altico proviene del drenaje en tierra llevado a cabo en todos los pa´ıses de la cuenca b´altica y suma un ´area unos 1.7×106km2, unas cuatro veces Materia en Suspensi´on y su Modelizaci´on 275 la superficie del Mar B´altico, m´as de tres veces la superficie de un pa´ıs como Espa˜na o Suecia. Como curiosidad [29] el volumen de agua recogido por todos r´ıos b´alticos es similar al de uno s´olo de los grandes r´ıos del planeta, como el Mississippi. En la figura [68] se aprecia que son nueve los pa´ıses costeros con el B´altico: Suecia, Finlandia, Rusia (territorios de Leningrado y Kaliningrado), Estonia, Letonia, Lituania, Polonia, Alemania y Dinamarca. Figura 8.5: Mapa ilustrativo de las cuencas hidrogr´aficas del Mar B´altico. 276 CAP´ ITULO 8 Para completar los de la cuenca b´altica hay que a˜nadir otros cuatro pa´ıses, alguno de cuyos r´ıos termina vertiendo en el sistema que nos ocupa: Rep´ublica Checa, Eslovaquia, Ucrania y Bielorrusia. Ucrania y Bielorrusia tambi´en tienen tierras que dan a la cuenca del Mar Negro, entre las que se incluyen las afectadas por el desastre de Chernobyl, que drenaron primero hacia el r´ıo Pripyat, afluente del r´ıo Dnieper, perteneciente a la cuenca del Mar Negro. Conectando el flujo de agua Φriv wy el de masa suspendida Φriv mpodemos hacer una estimaci´on sobre la concentraci´on de la materia suministrada por los r´ıos: Φriv w∼1.5×104m3s–1 ,(8.36) Φriv m=ρriv mΦriv w∼1010 kg a–1 =⇒ρriv m∼20 g m–3 .(8.37) Otro hecho destacable es que existen un grupo reducido de unos cinco r´ıos catalogables como medianos-grandes y un elevado n´umero (unos sesenta) r´ıos medianos-peque˜nos que riegan de una manera m´as o menos uniforme [125] las costas b´alticas. Dentro del grupo de los r´ıos mayores los que m´as van a afectar a nuestro estudio son: Neva, V´ıstula, ´ Oder, Neman, Kemijoki y Angermanalven. Para recoger esto en el modelo, trabajaremos teniendo en cuenta la aportaci´on y flujos de los r´ıos mayores como fuentes puntuales, y trataremos el resto de las aportaciones como un continuo costero ponderando las descargas de los r´ıos peque˜nos en cada zona. A partir de ahora llamaremos r´ıos grandes a los que modelemos como fuente puntual y peque˜nos a los que lo hagamos con fuente continua. Esta superposici´on de una distribuci´on discreta de fuentes m´as una distribuci´on continua para modelar los r´ıos peque˜nos es razonable dado el n´umero de r´ıos peque˜nos, la resoluci´on espacial de nuestro modelo y las escalas temporales con las que vamos a trabajar. Tendremos unos 60 r´ıos peque˜nos y un contorno costero de unos 3000 km, lo que supone un r´ıo peque˜no cada dos o tres nudos de nuestro mallado, que tiene 20 km de resoluci´on espacial horizontal. Para acercarnos un poco m´asalarealidaddeestosr´ıos, podemos hacer una comparaci´on entre ellos y los caudales medios anuales de algunos r´ıos espa˜noles, anotando que Materia en Suspensi´on y su Modelizaci´on 277 lugar ocupar´ıan estos ´ultimos en el ranking de los 65 r´ıos b´alticos: el r´ıo Tajo ocupar´ıa el lugar 9o, el Guadalquivir estar´ıa en el puesto 25o, mientras que J´ucar y el Segura estar´ıan en los puestos 47oy58 orespectivamente. En la tabla siguiente aparecen los cinco r´ıos que m´as influencia van a tener en nuestro estudio. TABLA 8.2: Caudal de los grandes r´ıos de la cuenca b´altica R´ıo Neva V´ıstula Neman ´ Oder Kemijoki ˚ Angerman¨alven Aporte medio ( m3s–1 )2460 1065 632 573 562 489 Si distribuimos los r´ıos por zonas compartimentales veremos una mayor aportaci´on relativa mayor en las zonas del norte (Mar de Botnia, Bah´ıa de Botnia) y Golfo de Finlandia respecto al B´altico Estricto (grueso del B´altico Central y Sur), con una ponderaci´on por metro de l´ınea de costa que puede ser entre dos y cinco veces superior, como puede comprobarse en la siguiente tabla. TABLA 8.3: Aporte de agua dulce a los compartimentos b´alticos Compartimento Aporte fluvial ( m3s–1 ) Bah´ıa de Botnia 3104 Mar de Botnia 2860 Golfo de Finlandia 3556 Golfo de Riga 1020 B´altico Estricto 3610 Estrechos de Dinamarca, Kattegat 1159 TOTAL 15310 Si ponderamos la aportaci´on de los r´ıos con los vol´umenes de cada zona la desigualdad anterior es mayor, debido a que la zona B´altico Estricto es la m´as profunda y voluminosa, como puede comprobarse a continuaci´on. Para tener presente la variaci´on mensual del caudal de los r´ıos, veamos el siguiente el perfil (figura 8.6), que es una curva media calculada a partir de datos experimentales [209, 125]. 278 CAP´ ITULO 8 TABLA 8.4: Aporte de agua dulce en relaci´on al volumen de los compartimentos b´alticos Compartimento Incremento de volumen anual por r´ıos (%) Mar de Botnia y Bah´ıa de Botnia 2.8 Golfo de Finlandia 10.7 Golfo de Riga 6.5 B´altico Estricto 0.9 TOTAL 2.5 Aunque se representa un perfil medio (figura 8.6), existen peque˜nas diferencias nortesur, debido a una mayor y m´as tard´ıa subida del caudal de los r´ıos del norte, cuyas aguas primaverales provienen en gran proporci´on de los deshielos de estas zonas fr´ıas. 0 2 4 6 8 10 12 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 month of the year Monthly inflow / Annual inflow Monthly mean river inflow Annual mean inflow Figura 8.6: Caudal medio mensual de los r´ıos b´alticos en unidades relativas de caudal medio anual. Materia en Suspensi´on y su Modelizaci´on 279 4.1.2 Fuentes de producci´on primaria La informaci´on sobre las fuentes de producci´on primaria en el Mar B´altico se basan fundamentalmente en dos trabajos previos: (a) Un estudio sobre materia de suspensi´on en el Mar B´altico de Hakanson [58] donde se analizan y modelan con bastante detalle peque˜nos sistemas (2 −20 km2) repartidos por diferentes latitudes de este mar. En este trabajo se establece una estrecha correlaci´on matem´atica entre la tasa de producci´on primaria y la densidad de clorofila en el B´altico. (b) Un modelo sobre la din´amica anual de la clorofila en el Mar B´altico [37] de donde salen unos mapas con la densidad de clorofila mes a mes a lo largo del a˜no. 0 2 4 6 8 10 12 −10 −5 0 5 10 15 20 month of the year Monthly mean temperature ( o C ) Baltic Proper south Baltic Proper north Gulf of Finland Bothnian Sea Bothnian Bay Figura 8.7: Curvas de temperatura media mensual [108] en diferentes zonas del Mar B´altico. Aunque la relaci´on es m´as complicada, se ha simplificando sustituyendo variables por valores medios que aparecen en el propio estudio ya citado [37], llegando a la siguiente ecuaci´on para la tasa producci´on org´anica ,Φorg,en funci´on de la temperatura y la densidad de clorofila ρchl: 280 CAP´ ITULO 8 Φorg =cTaρb chl (8.38) donde: Φorg =Φ org(g m–2 mes–1),ρ chl =ρchl(mg m–3),(8.39) a=1 ,b=0.93 ,c=0.35 ,(8.40) T=1 10 +θ×max (0,sign θ) siendo θ=θ(oC) .(8.41) La temperatura θ(oC) es la temperatura media mensual de la superficie del agua. Como indicadores de estas temperaturas se pueden usar los anuarios de temperaturas de los institutos Suecia y Finlandia o el trabajo recopilatorio de Landsberg [108]. En la figura 8.7 se muestran unas gr´aficas sobre las temperaturas medias mensuales en una estaci´on meteorol´ogica representativa de cada zona del B´altico: B´altico Estricto Sur, B´altico Estricto Norte, Golfo de Finlandia, Mar de Botnia, Bah´ıa de Botnia. Respecto a las densidades de clorofila (mg m–3) a lo largo del a˜no, se van a utilizar unas cartas mensuales existentes en la literatura [37] tomando valores medios en cada compartimento del B´altico y teniendo en cuenta los m´aximos locales en las desembocaduras de los r´ıos m´as importantes. 4.2. Ajuste de la tensi´on cr´ıtica de resuspensi´on al tipo de fondo marino Para acercarnos a esta cuesti´on es necesario describir aunque sea de una manera somera las caracter´ısticas principales del fondo marino b´altico. La clasificaci´on de la superficie sedimentaria m´as usada en este sistema es la llamada aproximaci´on gen´etica [78] basada en las caracter´ısticas ac´usticas del sedimento. Los tipos b´asicos de suelo marino seg´un esta clasificaci´on son los siguientes: (a) Suelo duro (Hard bottom): formado por piedras y gravas gruesas. (b) Suelo arenoso (Sand bottom): formado por arena y gravas finas. (c) Suelo suave (Soft bottom): formado por limos y arcillas. Cuando se hacen las mediciones sobre una zona concreta del fondo del mar se encuentra en muchos casos que la se˜nal contiene informaci´on combinada de dos tipos distintos Materia en Suspensi´on y su Modelizaci´on 281 de suelos. As´ı, podemos hablar de suelos binarios. Los suelos binarios pueden presentarse compensados (igual proporci´on entre sus partes) o descompensados (uno de los dos suelos es mayoritario). De ah´ı, que el mapa detallado del fondo marino b´altico [115, 78] contenga doce variedades de suelos, que se muestran en la siguiente tabla: TABLA 8.5: Tipos b´asicos de fondos marinos b´alticos y sus abreviaturas Tabla de suelos SOFT SAND HARD SOFT SF SF-SN , SF-sn , SN-sf SF-H , SF-h , sf-H SAND -SN SN-H , SN-h , sn-H HARD - - H Abreviaturas: SOFT = SF , soft = sf ; SAND = SN , sand = sn ; HARD = H , hard = h. La letra min´uscula indica que se encuentra en una proporci´on menor. A partir de esta bibliograf´ıa se ha elaborado, tras un largo trabajo de adaptaci´on manual y digital a nuestro mallado, un mapa con los cinco tipos de suelos m´as abundantes, recogido todo ello en la figura 8.8. Para llegar a este mapa, han sido despreciados los suelos muy minoritarios, a la vez que se han agrupado entre s´ı(porejemplo,SN-h−→ SN) algunas clases de suelos por su muy similar comportamiento respecto a la din´amica de la materia en suspensi´on, que es lo que aqu´ı nos convoca. Los cinco tipos de suelos que van a aparecer en nuestro modelo son: SOFT SOFT-HARD SAND HARD-soft HARD . El origen estos suelos binarios con una estructura tan irregular a peque˜na escala, donde conviven texturas tan distintas, obedece al pasado glacial relativamente reciente del Mar B´altico, del que ya se habl´oenelcap´ıtulo 3. La ´ultima glaciaci´on afect´o severamente al B´altico Norte, y el sistema no termin´o de formarse como tal hasta hace unos siete mil a˜nos (figura 3.3 del cap´ıtulo 3). Antes de adquirir su actual condici´on de mar salobre (7000 B.P.) pas´o incluso por una fase de lago (Ancylus Lake; [41]) durante unos dos mil a˜nos, hasta que se abrieron los estrechos de Dinamarca, gracias a la recuperaci´on de la forma de los suelos B´altico adentro. 288 CAP´ ITULO 8 3 32.5 3 2.5 32.5 2.5 3 4 4 3 3 2.5 2.5 3 455 6 10 5 5 4 5 4 4 4 5 6 6 6 5 3 32.5 3 2.5 32.5 2.5 3 4 4 3 3 2.5 2.5 3455 6 10 5 5 4 5 4 4 4 5 6 6 6 5 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 km km 5 66 6 5 4 4 4 5 4 5 510 6 5 5 4 3 2.5 2.5 3 3 44 3 2.5 2.5 3 2.5 3 2.5 3 3 0 2 4 6 8 10 12 Figura 8.10: Salida del modelo de transporte de materia en suspensi´on para la funci´on ρS(r,t) expresada en g m–3 y evaluada para el instante del ciclo con fecha 1 de febrero. El n´umero de part´ıculas trazadoras inicialmente colocadas en el sistemas ha sido N= 2096780 .Con muy buena aproximaci´on, los resultados de la funci´on de concentraci´on se repiten con periodicidad anual a partir del segundo ciclo anual de simulaci´on: Materia en Suspensi´on y su Modelizaci´on 289 ρS(r,t)≈ρS(r, t +T),t≥1a˜no ,T=1a˜no .(8.54) Se muestran a continuaci´on el resultado de las simulaciones para la evoluci´on de la materia en suspensi´on. 3 3.5 5 3 3.5 3 6 4.5 4.5 3.5 3 3 3 3.5 3.5 3 3.5 4 3.5 4 4 4 4 3 3.5 5 3 3.5 3 6 4.5 4.5 3.5 3 3 3 3.5 3.5 3 3.5 4 3.5 4 4 4 4 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 km km 4 4 4 4 3.5 4 3.5 33.5 3.5 3 3 3 3.5 4.5 4.5 6 3 3.5 3 53.5 3 0 2 4 6 8 10 12 Figura 8.11: Salida del modelo de transporte de materia en suspensi´on para la funci´on ρS(r,t) expresada en g m–3 y evaluada para el instante del ciclo con fecha 1 de agosto. Versiones anteriores de este modelo de materia en suspensi´on elaboradas por el autor de esta memoria se pueden encuentran en la bibliograf´ıa [191, 192]. 290 CAP´ ITULO 8 Los resultados del modelo son muy consistentes con valores de materia en suspensi´on obtenidos en el B´altico, tanto por medidas directas [58], como indirectamente a trav´es de relaciones emp´ıricas con la profundidad Secchi medida en el sistema [1, 162, 58, 106]. Podemos hacer uso de una relaci´on emp´ırica [58] entre la profundidad Secchi SEC (m) y la densidad de materia en suspensi´on SPM (g m–3) en superficie: SPM (g m–3)=10 A+B ×SEC (m) ,A=17 20 ,B=−1 20 ,(8.55) cuya confianza estad´ıstica es del 90 % . Se muestra a continuaci´on la tabla 8.7, que se ha calculado a partir de los datos de SEC (m) en diferentes compartimentos del Mar B´altico donde aparecen el valor medio y el rango de variaci´on de la densidad de materia en suspensi´on . TABLA 8.7: Profundidad Secchi y densidades medias de materia en suspensi´on en el Mar B´altico por compartimentos Zona B´altico Estricto Golfo de Finlandia Mar de Botnia Bah´ıa de Botnia SEC (m) 7.5 ±2.5 2.0 ±1.0 4.0 ±2.0 3.0 ±1.5 SPM (g m–3)3.0 ±1.5 5.5 ±1.5 4.5 ±1.5 5.0 ±1.5 Hay que mostrar tambi´en que en los datos experimentales se observa un claro gradiente negativo desde la costa al interior, o dicho de una manera m´as din´amica, desde las zonas de erosi´on a las zonas de sedimentaci´on, gobernadas ambas por tensiones de fondo cualitativamente distintas en relaci´on a las tensiones cr´ıticas. TABLA 8.8: Gradientes costa–interior e incrementos invierno–verano t´ıpicos de la densidad de materia en suspensi´on en el Mar B´altico por compartimentos Zona B´altico Estricto Golfo de Finlandia Mar de Botnia Bah´ıa de Botnia ΔxSPM (g m–3)[0.5 , 1.5] [1.0 , 3.0] [1.0 , 2.0] [1.0 , 2.5] ΔtSPM (g m–3)[–0.5 , 0.0] [–1.5 , –0.5] [–1.0 , –0.5] [–1.0 , –0.5] A partir de las publicaciones mencionadas [1, 162, 58, 106] se ha elaborado la tabla 8.8 donde se han estimado los rangos de los gradientes de SPM costa-interior, ΔxSPM (g m–3). Materia en Suspensi´on y su Modelizaci´on 291 Tambi´en se muestran valores t´ıpicos de la variaci´on invierno-verano, ΔtSPM (g m–3), que suele ser generalmente negativa por la mayor energ´ıa oce´anica durante el invierno aunque hay factores atenuantes como la producci´on primaria, mayor en las ´epocas suaves. Otra propiedad relacionada con la tabla anterior es que el gradiente espacial es menor en verano y m´as acusado en invierno. La zona del B´altico Estricto (compartimento m´as al sur) presenta un comportamiento distinto al resto del sistema. Por un lado, tiene menor aporte fluvial (favorecedor del gradiente espacial de SPM) y por otro, debido a sus condiciones t´ermicas m´as suaves, una producci´on primaria que se convierte en la fuente dominante de la zona, generando la materia de una manera m´as homog´enea que el aporte fluvial. Los dos m´aximos relativos del sur corresponden a la desembocadura del Oder y del V´ıstula, respectivamente de oeste a este. El m´aximo absoluto, en la parte m´as oriental del Golfo de Finlandia, es debido a las grandes descargas del r´ıo Neva, el mayor con diferencia de los r´ıos b´alticos. Se han mostrado los resultados de dos simulaciones del modelo de transporte de materia en suspensi´on en el Mar B´altico; la primera –figura 8.10– es una salida del modelo para invierno (1 de febrero), mientras que la segunda –figura 8.11– es otra salida, correspondiente a verano (1 de agosto). Por lo general, los valores de la funci´on densidad y los gradientes presentan un acuerdo aceptable con los valores experimentales mostrados en las tablas 8.7 y 8.8, m´axime teniendo en cuenta la naturaleza de baja resoluci´on de nuestro modelo num´erico: (a) Se alcanza un buen acuerdo en las densidades materia en suspensi´on si comparamos la figura 8.10 con los rangos de la tabla 8.7 acompa˜nada de la variaci´on espacial expresada en la tabla 8.8. (b) La figura 8.11 es una simulaci´on que representa valores alcanzados durante el verano y reproduce cualitativamente la variaci´on estacional de la materia en suspensi´on contemplada en la literatura y resumida en la tabla 8.8. (c) La escasa variaci´on estacional de materia en la zona del B´altico Estricto as´ıcomosu distribuci´on pr´acticamente homog´enea en verano se˜nalada en la literatura tambi´en se reproduce con aceptable correcci´on en el modelo. 292 CAP´ ITULO 8 Encuantoalaresoluci´on vertical de funci´on ρS(z),tambi´en se ha comprobado que hay un buen acuerdo, tomando para comparar la capa pr´oxima al sedimento de erosi´on y transporte, donde se aprecia una densidad 1–2 g m–3 por encima de la densidad en la capa superficial, de acuerdo con valores t´ıpicos de la literatura [58]. Tambi´en hay un buen acuerdo del modelo con los valores de ρSen la capa de fondos de acumulaci´on, que est´an [59] t´ıpicamente en el rango 0.5−1.5gm –3 . Veamos ahora la validaci´on del modelo con medidas de materia en suspensi´on (tabla 8.9) acumulada en el sedimento en diferentes puntos del sistema. A tenor de las caracter´ısticas de nuestro modelo num´erico de baja resoluci´on, las medidas que se van a usar para la validaci´on han sido seleccionadas descartando puntos experimentales [120] fuertemente influenciados por cuestiones locales como cercan´ıa a desembocaduras de r´ıos o pertenencia a peque˜nas fosas, factores ambos que elevan la tasa de sedimentaci´on por encima de los valores promedios zonales representados por algunas estaciones de medida encontradas en la literatura cient´ıfica y mostradas a continuaci´on. TABLA 8.9: Tasa de sedimentaci´on experimental [120, 202] en diferentes puntos del Mar B´altico elegidos para la validaci´on del modelo. Compartimento Posici´on Tasa (g m–2 a–1) Bah´ıa de Botnia 64.70oN , 22.06oE320 Bah´ıa de Botnia 65.00oN , 22.25oE220 Bah´ıa de Botnia 64.31oN , 22.36oE220 Mar de Botnia 62.60oN , 20.26oE150 Mar de Botnia 61.08oN , 20.60oE270 Golfo de Finlandia 59.58oN , 23.63oE400 Golfo de Finlandia 60.13oN , 27.46oE440 B´altico Estricto 56.94oN , 19.81oE160 B´altico Estricto 56.95oN , 20.26oE230 B´altico Estricto 57.31oN , 20.08oE140 El error experimental de los valores que aparecen en la tabla 8.9 es muy variado, si bien el error medio relativo es del orden de un treinta por ciento. Las tasas de sedimenta- Materia en Suspensi´on y su Modelizaci´on 293 ci´on han sido obtenidas mediante el modelo CRS (Constant Rate of Supply) de dataci´on por el is´otopo 210Pb. 200 200 200 100 100 200 200 100 100 200 400 200 300 300 100 100 100 200 200 200 100 100 200 200 100 100 200 400 200300 300 100 100 100 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 km km 100 100 100 300 300 200 400 200 100 100 200 200 100 100 200 200 200 50 100 150 200 250 300 350 400 450 458 / 400 196 / 220 235 / 150 180 / 230 192 / 150 419 / 440 188 / 270 203 / 220 284 / 320 191 / 140 Figura 8.12: Salida del modelo de transporte de materia en suspensi´on de la tasa de sedimentaci´on en el estado estacionario RS(g m–2 a–1).En el numerador se presenta el valor del modelo y en el denominador el de la medidas experimentales, Rmod S/Rexp S. En la figura 8.12 se aprecia que el acuerdo entre el modelo y las medidas es muy aceptable, pues se han descartado para validar puntos de fuerte influencia por r´ıos (ver figura 8.5) donde la tasa de sedimentaci´on puede llegar a ser superior a la tasa normal 294 CAP´ ITULO 8 en un factor 3 o superior, superando en algunos casos el nivel RS= 1000 g m–2 a–1 .Las posiciones de los puntos validadores del modelo son m´as interiores de lo que parece al ver la figura en la primera impresi´on, ya que es la punta rigurosa de la flecha la que los define; de hecho, todos los puntos est´an en la zona de acumulaci´on. 0 100 200 300 400 500 600 0 100 200 300 400 500 600 Experimental sedimentation rate (g m2 yr−1) Modelled sedimentation rate (g m2 yr−1) Figura 8.13: Gr´afica comparativa de la tasa de sedimentaci´on RS(g m–2 a–1)delmodelocon la experimental, donde se representan los diez puntos de validaci´on de la figura anterior 8.12. Estas diferencias entre los valores normales y los extremos pueden verse amplificadas por el efecto local “focusing” de peque˜na escala –no tenido en cuenta en este modelo de baja resoluci´on– que afecta a depresiones en forma de “U” o en forma de “V” y que seg´un Hakanson [58] puede aportar amplificaciones de RSen el rango 1–3. En la figura 8.13 se representan los 10 puntos de validaci´on de la figura anterior, 8.12, junto a la recta de acuerdo total entre modelo y experiencia, apreci´andose aceptables discrepancias modelo–experiencia en los dos sentidos. Se muestra a continuaci´on un resultado de la tasa de sedimentaci´on no estacionaria, soluci´on del problema que es interesante para identificar las zonas donde la materia en Materia en Suspensi´on y su Modelizaci´on 295 suspensi´on se almacena provisionalmente. Se trata de la salida del modelo para la tasa de sedimentaci´on RS(g m–2 a–1) no estacionaria, calculada entre el d´ıa 1 de agosto del ciclo y el 1 de febrero del ciclo siguiente. 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 km km 0 00 0 00 0 0 0 0 −400 −300 −200 −100 0 100 200 300 400 Figura 8.14: Salida del modelo de transporte de materia en suspensi´on de la tasa de sedimentaci´on RS(g m–2 a–1) no estacionaria, calculada entre el 1 de agosto del ciclo y el 1 de febrero del ciclo siguiente. Se tienen zonas tanto de signo positivo como negativo, pudiendo deducirse, por comparaci´on con la gr´afica anterior, el almacenamiento provisional en algunas zonas de mate- 296 CAP´ ITULO 8 ria durante primavera y verano que posteriormente, con la llegada del oto˜no y el invierno, es resuspendido y almacenado finalmente las zonas m´as profundas del sistema. En la figura 8.14 se representa la l´ınea de nivel 0, criterio elegido para dividir el fondo marino en dos partes con comportamientos din´amicos distintos durante el ciclo anual de la materia en suspensi´on en el mar B´altico. Los resultados que ofrecen las simulaciones de nuestro modelo son bastante congruentes con los datos experimentales, a pesar de la complejidad del sistema con el que estamos trabajando. Otros trabajos relacionados con la modelizaci´on de la din´amica de la materia en suspensi´on han sido publicados anteriormente [102, 59, 214] pero siempre estudiando una parte del Mar B´altico y no el sistema como un todo. Un objetivo cumplido de este trabajo es haber conseguido, con todas las limitaciones expuestas aqu´ı, un primer modelo, l´ogicamente de baja resoluci´on, de la evoluci´on de la materia en suspensi´on para todo el Mar B´altico, como sistema completo. Futuras ampliaciones del modelo basadas en aproximaciones m´as realistas y de mayor resoluci´on como las ya comentadas en esta tesis, apoyadas incluso en un sub–modelo fino de las corrientes en el fondo marino [72], har´an posible un acercamiento mayor del modelo a la realidad compleja del sistema. Cap´ıtulo 9 Modelo del Transporte No Conservativo de Radion´uclidos 1. INTRODUCCI´ ON En esta parte de la tesis se va a presentar un modelo de transporte de radion´uclidos no conservativos del Mar B´altico que ser´a validado con datos provenientes del accidente Chernobyl. Este estudio se apoya en anteriores partes de esta tesis como son, por un lado, el modelo de transporte conservativo que tiene en cuenta la evoluci´on de todo aquello que sigue fielmente los movimientos del fluido oce´anico, incluida la fase disuelta del radion´uclido, y por otro, el modelo de transporte materia en suspensi´on a la que se encuentra ligada la fase suspendida del mismo. Los radion´uclidos pueden cambiar de fase, pasando de la fase disuelta a una s´olida (adsorci´on), o al contrario (desorci´on), mecanismos de los que a continuaci´on daremos m´as detalles. 2. FASES: DISUELTA Y S´ OLIDA EN SUSPENSI´ ON O SEDIMENTO Los mecanismos de adsorci´on y desorci´on se desarrollan en unos procesos cin´eticos [157, 8, 9] con unos tiempos de relajaci´on que conducen al equilibrio, y que en el is´otopo que nos ocupa, 137Cs, presentan una fase de transici´on r´apida (minutos), una intermedia (horas) y una lenta (d´ıas). Si trabajamos con un paso de tiempo suficientemente grande, podremos apoyarnos en la hip´otesis de equilibrio entre la fase disuelta de 137Cs y las fases s´olidas (suspendida y en sedimento) a trav´es del coeficiente de distribuci´on kd: 297 304 CAP´ ITULO 9 De este modo, las probabilidades restantes de transici´on se calculan: p12 =max[0,−ε1,2],p 13 =max[0,−ε1,3].(9.35) (d) Desintegraci´on radiactiva [143]: pi4=1−e−λΔt,λ=ln 2 T1/2 ,i=1,2,3.(9.36) Esta forma matem´atica de expresar la probabilidad de transici´on es extensible [141] de una manera an´aloga al esquema cin´etico de no equilibrio, donde el papel jugado por la constante de desintegraci´on λlo jugar´ıan los coeficientes de adsorci´on (k1)y de desorci´on (k2). (e) Otras probabilidades simples: p4i=0,i=1,2,3; pii =1− 4  j=1 (1 −δij)pij ,i=1,2,3,4.(9.37) 4. MODELO DE TRANSPORTE,RESULTADOSYVALIDACI ´ ON El transporte advectivo de este modelo se rige por el esquema aplicado en el cap´ıtulo 7 – dedicado al transporte conservativo – tomando el modelo HD–3D de circulaci´on anual detallado en el cap´ıtulo 4. Debido al uso de un esquema de corrientes advectivas anuales, trabajaremos con la tabla de coeficientes de difusi´on (tabla 9.1.) elaborada para estos efectos en el cap´ıtulo 7, donde se muestran rangos de referencia para la difusi´on horizontal anual Kanu h(m2s–1)y para la vertical Kv(cm2s–1) en la superficie del mar, variables seg´un las estaciones el a˜no y de implementaci´on inmediata en las simulaciones de este cap´ıtulo. TABLA 9.1: Rangos de referencia para los coeficientes de difusi´on horizontal y vertical Tabla de Difusi´on Invierno Primavera Verano Oto˜no Kanu h(m 2s–1)400–800 200–400 150–300 300–600 Kv(cm 2s–1)15–30 5–10 1–3 10–20 Los coeficientes de difusi´on que se han usado para las capas no superficiales se han calculado a partir de los valores superficiales de la tabla de arriba teniendo en cuenta los Modelo del Transporte No Conservativo de Radion´uclidos 305 perfiles verticales de Kh(m2s–1)yKv(cm2s–1) que aparecen en las tablas del cap´ıtulo 5, que es la parte de esta tesis doctoral dedicada por entero a la difusi´on. En relaci´on a las datos de actividad de 137Cs en el sistema, tenemos que en Agosto de 1987 se realiz´o una expedici´on cient´ıfica llamada SWRV ATTAIR para analizar la evoluci´on de la actividad por comparaci´on con las medidas de 1986, gran parte de las cuales provienen del crucero cient´ıfico RV GAUSS. Los datos de este ´ultimo han sido muy ´utiles para elaborar las condiciones iniciales de la actividad en octubre de 1986, mientras que los datos de la SWRV ATTAIR junto con otros forman un conjunto [53, 54, 131, 66] de datos experimentales cuyas caracter´ısticas fueron expuestas en el cap´ıtulo 7, conjuntamente a su adaptaci´on a nuestra malla de c´alculo cuyos resultados han sido validados por publicaciones internacionales relevantes. Respecto a la actividad previa en agua y en sedimentos, el dato importante es que la actividad de 137Cs anterior a Chernobyl es [76, 66] aproximadamente uniforme en el B´altico y pr´acticamente independiente de la profundidad. El valor de fondo de referencia se sit´ua en ∼10–20 Bq m–3 . Es destacable que las medidas de 137Cs realizadas en oto˜no de 1986 en puntos por debajo de la picnoclina permamente indican unos valores del mismo orden que los niveles de fondo antes citados, lo que indica un efecto barrera a considerar, a la vez que la relevancia en estas columnas de l´ıquido del transporte de 137Cs a trav´es de la materia en suspensi´on como principal mecanismo de dispersi´on. El gradiente de actividad disuelta es grande entre los puntos a ambos lados de la picnoclina y la mayor´ıa de los radion´uclidos atraviesan la barrera ligados a la materia en suspensi´on; una vez en las capas m´as profundas una parte de ellos va pasando a la disoluci´on hasta que se alcanza el equilibrio, mientras que el resto queda ligado a la materia suspendida que puede acabar definitivamente en el sedimento. Se ha analizado abundante informaci´on de diferentes cores de sedimentos [76, 66] en el B´altico y la actividad integrada –anterior al accidente de Chernobyl– en el espesor de mezcla es muy congruente con este valor de fondo de ∼10–20 Bq m–3 . Otro valor de referencia [65, 66, 67] para la actividad de 137Cs contenida en un core de 10 cm y anterior al accidente de Chernobyl es ∼500–2000 Bq m–2 . Una ratio de referencia para la sedimentaci´on (medida con 210Pb) en el Mar B´altico [202] es ∼0.5–1 cm a–1 .Todaesta informaci´on se ha utilizado para calcular la condiciones iniciales de actividad de 137Cs en los sedimentos. 306 CAP´ ITULO 9 La asignaci´on de part´ıculas representativas del proceso de Monte Carlo se ha hecho asignando inicialmente 1 part´ıcula por cada Bq m–3 contenido cada unidad de volumen de dimensiones hx×hy×5m , distribuidas uniformemente por las celdas del modelo. El n´umero necesario de part´ıculas trazadoras de 137Cs inicialmente en disoluci´on ha sido Nsol(0) = 1786746 . An´alogamente y tal como se vio en cap´ıtulo 8, se han distribuido part´ıculas trazadoras de la materia en suspensi´on tanto en agua como en sedimento, asignando posteriormente las part´ıculas trazadoras no conservativas de 137Cs a partir de la hip´otesis inicial de equilibrio dada por kd. Se ha usado la metodolog´ıa num´erica lagrangiana expuesta en el cap´ıtulo 6, tomando para las simulaciones los siguientes pasos de tiempo para la circulaci´on y la difusi´on respectivamente: Δtcir =6h ,Δtdif =1h.(9.38) Los intercambios de radion´uclidos no conservativos se han computado bajo la hip´otesis de equilibrio con un paso de tiempo ΔtNC =2d´ıas . El resto de par´ametros f´ısicos y num´ericos usados son los que se exponen en la simulaci´on tridimensional sobre el Mar B´altico efectuada en el apartado 3 del cap´ıtulo 7, y en los detalles del modelo de transporte de materia en suspensi´on del cap´ıtulo 8. El tiempo de CPU necesario –en el ordenador personal especificado en el cap´ıtulo 6– para la simulaci´on de 300 d´ıas de transporte ha sido del orden de 6 horas. El resultado del modelo, simulando desde el 15 de Octubre de 1986 hasta el 15 de Agosto de 1987, tras computar el transporte de las part´ıculasytodalainformaci´on detallada aqu´ı, se muestra en la figura 9.1. Con todas estas consideraciones se ha hecho la simulaci´on y calculado la figura 9.1, que ha de ser comparada, a efectos de validaci´on, con la figura experimental que le sigue inmediatamente (figura 9.2). Modelo del Transporte No Conservativo de Radion´uclidos 307 50 50 100 200 200 200 200 300 300 400 200 100 50 50 100 200 200 200 200 300 300 400 200 100 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 km km 100 200 400 300 300 200 200 200 200 100 50 50 0 100 200 300 400 500 600 700 800 900 1000 Figura 9.1: Actividad (Bq m–3)de137Cs en la superficie del Mar B´altico resultado de la simulaci´on del modelo, para el 15 de Agosto de 1987. Esta gr´afica ha de utilizarse para validar el modelo con la siguiente gr´afica experimental (figura 9.2). 308 CAP´ ITULO 9 100 100 200 300 300 200 200 300 300 200 100 100 100 200 300 300 200 200 300 300 200 100 km km SWEDEN FINLAND BALTIC REPUBLICS 100 200 300 300 200 200 300 300 200 100 100 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 0 100 200 300 400 500 600 700 800 900 1000 Figura 9.2: Actividad (Bq m–3)de137Cs en la superficie del Mar B´altico durante Agosto de 1987, interpolada a partir de medidas experimentales contenidas en la literatura [53, 54, 131], ´util para la validaci´on del modelo la superficie del sistema. Modelo del Transporte No Conservativo de Radion´uclidos 309 El acuerdo entre el modelo y los valores experimentales es muy aceptable. Vuelven a aparecer discrepancias en la zona sureste, como ha ocurrido aqu´ı en la validaci´on en superficie, motivadas por las limitaciones de nuestro modelo de corrientes medias anuales. A continuaci´on se expondr´aaqu´ı otro proceso de validaci´on mediante la tasa de actividad de 137Cs, ΦA(Bq m–2 a–1),sedimentada en ocho puntos del B´altico correspondientes a zonas de acumulaci´on, seleccionados adecuadamente a la naturaleza de baja resoluci´on de nuestro modelo. Se ha seguido una metodolog´ıa similar a la del cap´ıtulo 8 a la hora de seleccionar los puntos de simulaci´on, descartando [120, 58] localizaciones sensibles a efectos locales como la cercan´ıa a desembocaduras de r´ıo o la pertenencia a fondos marinos con ”focusing costero” [59] relativamente frecuentes en muchos suelos mixtos del Mar B´altico. Estos datos experimentales est´an contenidos en la literatura [120, 66, 67] y van a ser posteriormente comparados con el resultado de nuestro modelo cuyo intervalo de simulaci´on va desde el a˜no 1986 al 1987. La mayor´ıa de los cores elegidos para la validaci´on coinciden con los del modelo de transporte de materia en suspensi´on del cap´ıtulo anterior y muestran la radiactividad acumulada durante una d´ecada (1986–1996) a partir de la que se ha calculado la parte correspondiente al a˜no 1986–1987 teniendo en cuenta los siguientes factores: (a) Estimaci´on de la distribuci´on de la radiactividad a lo largo de los a˜nos de la d´ecada en cuesti´on utilizando como variable de ponderaci´on los niveles t´ıpicos de 137Cs (Bq m–3) disuelto en la capa de fondo en contacto con el sedimento, informaci´on disponible en la literatura [69]. (b) Decaimiento radiactivo del is´otopo 137Cs en el intervalo considerado. En lo que sigue (tabla 9.2), se presenta la informaci´on relacionada con datos experimentales [120, 66, 67] de actividad de 137Cs sedimentada, ΦA(Bq m–2 a–1), en diferentes puntos del Mar B´altico fruto de la selecci´on descrita anteriormente. El error experimental de los valores que aparecen en la tabla 9.2. oscila t´ıpicamente entre el 10 y el 20 por ciento. 310 CAP´ ITULO 9 TABLA 9.2: Tasa de actividad de 137Cs sedimentada, ΦA(Bq m–2 a–1), durante el intervalo de simulaci´on, en diferentes puntos del Mar B´altico [120, 66, 67] elegidos para la validaci´on del modelo. Compartimento Posici´on ΦA(Bq m–2 a–1) Bah´ıa de Botnia 64.70oN , 22.06oE1180 Bah´ıa de Botnia 65.00oN , 22.25oE745 Bah´ıa de Botnia 65.13oN , 23.56oE520 Mar de Botnia 62.60oN , 20.26oE1120 Mar de Botnia 61.08oN , 20.60oE1420 Golfo de Finlandia 59.58oN , 23.63oE1125 Golfo de Finlandia 60.24oN , 27.72oE1255 B´altico Estricto 54.70oN , 19.35oE130 B´altico Estricto 57.31oN , 20.08oE52 B´altico Estricto 58.65oN , 18.50oE280 Cabe destacar los valores tan bajos alcanzados en los sedimentos de la zona B´altico Estricto, en alg´un caso por debajo de 1 orden de magnitud respecto a otros puntos de la tabla. Son tres los factores que ayudan a crear esta diferencia: (a) A excepci´on de la zona norte que limita con Mar de Botnia, los valores de 137Cs en agua son inferiores o muy inferiores al resto. (b) El efecto barrera de la picnoclina, que sobre todo en zonas muy profundas, hace que el aporte de 137Cs a las capas inferiores en un periodo de un a˜no sea del mismo orden que la actividad anterior al accidente de Chernobyl. (c) La salinidad media de este compartimento es superior al resto por su influencia atl´antica, lo que sube su coeficiente kdy aumenta la capacidad de almacenar radion´uclidos en disoluci´on a costa de disminuir la proporci´on suspendida y su futura sedimentaci´on. En general, el acuerdo entre el modelo no conservativo y las medidas experimentales es bastante aceptable a tenor de la complejidad del problema que estamos estudiando. Modelo del Transporte No Conservativo de Radion´uclidos 311 1000 500 1500 2000 100 100 100 500 1000 1500 1000 500 1000 500 1000 500 1500 2000 100 100 100 500 1000 1500 1000 500 1000 500 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 km km 500 1000 500 1000 1500 1000 500 100 100 100 2000 1500 500 1000 500 1000 1500 2000 2500 708 / 745 1112 / 1180 1482 / 1255 1274 / 1420 276 / 280 77 / 52 69 / 130 1416 / 1125 732 / 520 943 /1120 Figura 9.3: Salida del modelo de transporte de radion´uclidos. La gr´afica representa la actividad de 137Cs sedimentada, Φmod A(Bq m–2 a–1) en el intervalo de simulaci´on. En los or´ıgenes de las flechas el numerador representa el valor del modelo y el denominador el de la medidas experimentales, Φmod A/Φexp A.El extremo de la flecha est´ase˜nalando la posici´on del punto elegido para la validaci´on, simbolizado por una estrella roja de ocho puntas. 312 CAP´ ITULO 9 0 500 1000 1500 2000 0 200 400 600 800 1000 1200 1400 1600 1800 2000 Actividad sedimentada (Bq m−2 a−1) experimental Actividad sedimentada (Bq m−2 a−1) modelada Figura 9.4: Comparaci´on entre la actividad experimental sedimentada de 137Cs, medida en diferentes estaciones del B´altico y la salida del modelo. Se dibuja en azul la bisectriz del cuadrante materializando el acuerdo total entre el modelo y la experiencia. A pesar de este buen acuerdo general, la mayor discrepancia se produce en el punto m´as al sur, 69/130 , debido fundamentalmente, seg´un vimos en el cap´ıtulo 7, a una desviaci´on fruto de un flujo de corrientes advectivas del norte al sur, no contemplado en nuestro modelo de corrientes medias anuales, que pudo aumentar la actividad en esta zona coloc´andola muy por encima del modelo. En la zona sur, en general, la actividad sedimentada es menor por las tres razones expuestas anteriormente, siendo por tanto m´as sensible num´ericamente a las fluctuaciones de nuestro m´etodo de Monte Carlo, por lo que las l´ıneas de nivel de la figura 9.3 tienden a quebrarse m´as en esta parte del sistema. Hay comentar el resultado del modelo en el Golfo de Finlandia, zona tambi´en influenciada por una picnoclina permamente (en la cota aproximada de los cincuenta metros), Modelo del Transporte No Conservativo de Radion´uclidos 313 aunque no tan fuerte como la de B´altico Estricto debido sobre todo a su menor profundidad con espesores en la capa inferior de unos 20–40 m tales que la densidad de actividad varios meses despu´es del accidente de Chernobyl fue notable aunque inferior a la de otras zonas del golfo no sometidas a esta influencia. Recordemos que la estructura de la actividad en superficie en oto˜no de 1986 era una mancha central en el golfo. El efecto del modelo de transporte en esta parte del sistema sobre esta funci´on de actividad, con forma aproximada de campana, es una distribuci´on de actividad en el sedimento cuya caracter´ıstica m´as destacable, si miramos con atenci´on, es una especie de anillo de m´ınimos relativos que nos recuerda un poco la forma de las funciones potenciales con forma de sombrero mexicano. De los dos puntos de validaci´on del Golfo de Finlandia, el punto m´as al este, experiment´om´as sedimentaci´on (1255 Bq m–2 a–1) que el del oeste (1125 Bq m–2 a–1) a pesar de tener una actividad inicial en superficie del orden de la mitad y una similar tasa de sedimentaci´on de materia en suspensi´on Rs(g m–2 a–1). Nuestra conclusi´on es que la picnoclina es la responsable de esta diferencia cualitativa y nuestro modelo, a pesar de su naturaleza de baja resoluci´on, consigue reproducir este efecto razonablemente bien, resolvi´endolo y estableciendo diferencias contrastadas dentro de un mismo compartimento. Por tanto, subrayemos finalmente la capacidad que tiene nuestro trabajo, dentro de las limitaciones expuestas, de aportar esta herramienta relativamente poco presente en la literatura cient´ıfica como es la modelizaci´on del transporte de radion´uclidos bajo condiciones dominantes de picnoclina permamente. El siguiente paso del modelo en el futuro ser´ıa abordar simulaciones para periodos de varios a˜nos, del orden de una d´ecada e incluso m´as all´aeneltiempo. Estoser´ıa factible a corto plazo mejorando la eficiencia de los algoritmos para no tener problemas de memoria en las simulaciones, teniendo definitivamente que abrir el sistema al Oc´eano Atl´antico para periodos de varios a˜nos en que su influencia empieza a no ser despreciable, sobre todo en la zona sur, habida cuenta de los tiempos de residencia del Mar B´altico que para una part´ıcula media van de dos a tres d´ecadas. Ser´ıa interesante llegar a una soluci´on estacionaria de la distribuci´on de 137Cs bajo la hip´otesis de mar cerrado, e incluso, compararla con la misma para el caso de sistema abierto. Otra opci´on de trabajo futuro consiste en la posibilidad, si lo permite la red de estaciones de monitorizaci´on que coordina la HELCOM (Helsiki Commision; Baltic Marine Environment Protection Commission), de sacar una rama de nuestro modelo que permita hacer predicciones a˜no a a˜no como 320 CAP´ ITULO 10 mente [102, 58, 214] pero siempre estudiando una parte del Mar B´altico y no el sistema como un todo. Un objetivo cumplido de este trabajo es haber conseguido, con todas las limitaciones expuestas aqu´ı, un primer modelo, l´ogicamente de baja resoluci´on, de la evoluci´on de la materia en suspensi´on para todo el Mar B´altico, como sistema completo. 11. El modelo reproduce bien las densidades de materia en suspensi´on en la superficie del mar B´altico, habiendo un buen acuerdo con los valores medios experimentales, con los gradientes t´ıpicos costa–interior y con las variaciones a lo largo de las diferentes estaciones del a˜no. 12. La validaci´on en los sedimentos, teniendo en cuenta que se trata de un modelo de baja resoluci´on y que siempre se han elegido puntos no sometidos a fluctuaciones locales, es muy aceptable, si tenemos en cuenta la dificultad del problema que estamos estudiando. 13.Enelcap´ıtulo 9, el acuerdo entre el modelo no conservativo y las medidas experimentales es bastante aceptable. El mayor desacuerdo se produce en el punto de validaci´on m´as al sur, debido fundamentalmente, seg´un vimos en el cap´ıtulo 7, a una desviaci´on fruto de un flujo de corrientes advectivas del norte al sur, no contemplado en nuestro modelo de corrientes medias anuales, que pudo aumentar la actividad en esta zona coloc´andola muy por encima del modelo. 14. Hay que comentar el resultado del modelo en el Golfo de Finlandia, zona tambi´en influenciada por una picnoclina permamente, aunque no tan fuerte como la del B´altico Estricto. Recordemos que la estructura de la actividad en superficie en oto˜no de 1986 era una mancha central en el golfo. El efecto del modelo de transporte en esta parte del sistema sobre esta funci´on de actividad, con forma aproximada de campana, es una distribuci´on de actividad en el sedimento cuya caracter´ıstica m´as destacable, si miramos con atenci´on, es una especie de anillo de m´ınimos relativos que nos recuerda un poco la forma de las funciones potenciales con forma de sombrero mexicano. De los dos puntos de validaci´on del Golfo de Finlandia, el punto m´as al este, experiment´om´as sedimentaci´on de radiactividad que el del oeste a pesar de tener una actividad inicial en superficie del orden de la mitad y una similar tasa de sedimentaci´on de materia en suspensi´on Rs(g m–2 a–1). Nuestra conclusi´on es que la picnoclina es la responsable de esta diferencia cualitativa y nuestro modelo, a pesar de su naturaleza de baja resoluci´on, consigue reproducir este Sumario, Conclusiones y Trabajo Futuro 321 efecto razonablemente bien, resolvi´endolo y estableciendo diferencias contrastadas dentro de un mismo compartimento. 15. Por tanto, subrayemos finalmente la capacidad que tiene nuestro trabajo, dentro de las limitaciones expuestas, para aportar esta herramienta relativamente poco presente en la literatura cient´ıfica como es la modelizaci´on del transporte de radion´uclidos bajo condiciones dominantes de picnoclina permamente, adem´as de las ventajas expuestas en los puntos anteriores. 10.2. TRABAJO FUTURO EN RELACI´ ONAESTATESISDOCTORAL La tesis doctoral cuya memoria se presenta aqu´ı es una especie de ´arbol con una serie de ramas por varias de las cuales podr´ıa avanzarse tanto como para desarrollar, con empe˜no, paciencia y curiosidad suficiente, otra tesis doctoral. Nos encontramos por tanto, como no pod´ıa ser de otro modo, con una estructura casi fractal, como ocurre tambi´en en una escala superior con la estructura reticular de la ciencia como superestructura que abarca todo el conocimiento humano. A continuaci´on, se ir´an apuntando las l´ıneas de trabajo futuro conforme han ido surgiendo en el desarrollo de los cap´ıtulos que vertebran esta obra. 1.Comol´ınea general de trabajo futuro, es interesante la idea de centrar todos los esfuerzos y las mejoras de la metodolog´ıa de esta tesis doctoral en una parte del Mar B´altico, el GolfodeFinlandia,porunaseriederazonesquetambi´en se exponen a continuaci´on. Aqu´ısere´unen una buena parte de los trabajos futuros relacionados con ampliaciones y refinamientos de las metodolog´ıas contenidas en los cap´ıtulos 4–7: (a) El Golfo de Finlandia es una parte del sistema que representa como ninguna otra, los procesos f´ısicos y medioambientales del Mar B´altico: alta concentraci´on de materia en suspensi´on, picnoclina permanente y mar generalmente helado en invierno, entre otros. (b)SielMarB´altico es uno de los mares m´as contaminados del mundo por su naturaleza estrechamente semicerrada, el Golfo de Finlandia es su compartimento m´as representativo por los aportes del mayor r´ıo, con diferencia, de la cuenca b´altica: el r´ıo Neva. 322 CAP´ ITULO 10 (c) El retorno de la inversi´on en trabajo nuevo de campo proveniente de la mejora de la resoluci´on espacial y temporal del modelo se producir´ıa en un tiempo razonable ypermitir´ıa pasar del actual modelo de baja resoluci´on a uno de una generaci´on posterior, con mayor resoluci´on, similar al HIROMB del SMHI. (d) Se podr´ıa dise˜nar en un tiempo aceptable un modelo hidrodin´amico propio de corrientes inducidas por vientos que nos dar´ıa independencia respecto al SMHI. Una opci´on ser´ıa un modelo POM (Princeton Ocean Model) como el modelo interesante que Zhurbas [213] ha desarrollado para el sur del B´altico, con circulaci´on inducida por los ocho vientos principales. (e) Riqueza de aplicaciones del m´etododeMonteCarloenlazonadelestuariodelr´ıo Neva. (f) El Golfo de Finlandia es una zona abierta al B´altico pero con un contorno no dif´ıcil de modelar num´ericamente. La dispersi´on de la pluma contaminante del r´ıo Neva se podr´ıa tratar, como primera aproximaci´on, mediante un sistema cerrado. (g) El Golfo de Finlandia es el eje marino y comercial m´as importante del Mar B´altico, rodeado de ciudades notables (Estocolmo, San Petersburgo, Tallin y Helsinki), con intenso tr´afico mar´ıtimo y problemas derivados de los vertidos petr´oleo, susceptibles de ser modelados. (h) Importantes centros de investigaci´on de la regi´on –SMHI, Finnish Institute of Marine Research, Uppsala University, entre otros– podr´ıan prestar y solicitar colaboraci´on yapoyo. 2. Un conocimiento m´as profundo sobre el estuario del R´ıo Guadalquivir es demandado para aplicaciones tanto socioecon´omicas como medioambientales: labores peri´odicas de dragado para mantener su navegabilidad, agricultura intensiva, turismo, y el impacto ecol´ogico sobre el Parque Nacional de Do˜nana, entre otros. Para este complejo sistema, una ampliaci´on futura de esta primera versi´on del modelo de transporte ser´a necesaria. 3. Son prometedoras las posibilidades de una l´ınea de trabajo futura en ETSI de la Universidad de Sevilla –situada en el punto kilom´etrico 94 del estuario– con el dise˜no de una l´ınea de proyectos fin de carrera que en el futuro pueden ir evolucionando hacia un proyecto de investigaci´on mayor que todos creemos necesario. Ya hay en marcha un proyecto fin de carrera dirigido por m´ıenestal´ınea y la direcci´on de la escuela apoya la idea de trabajar Sumario, Conclusiones y Trabajo Futuro 323 en este terreno, sobre todo pensando en las futuras promociones de alumnos de ingenier´ıa civil que ya se est´an formando en el centro, que empezar´an a trabajar en sus proyectos fin de carrera a partir del a˜no que viene y que en el futuro podr´an aportar y recibir mucho en relaci´on a este campo. Adem´as, existe un acuerdo con la profesora coordinadora de las asignaturas de hidr´aulica de esta titulaci´on para trabajar en una l´ınea com´un. 4. Otras colaboraciones que podr´ıan ser interesantes en el futuro en relaci´on con la investigaci´on del estuario del r´ıo Guadalquivir son las relacionadas con los grupos andaluces implicados en el informe sobre el estuario [159]. 5. Futuras ampliaciones del modelo de materia en suspensi´on del B´altico, basadas en aproximaciones m´as realistas y de mayor resoluci´on como las ya comentadas en esta tesis, apoyadas incluso en un sub–modelo fino de las corrientes en el fondo marino (Hille et al., 2006), har´an posible un acercamiento mayor del modelo a la realidad compleja del sistema. 6. Otra versi´on futura del modelo de materia en suspensi´on bajo diferencias finitas ser´ıa someterlo a advecci´on media anual (cap´ıtulo 7) corregida a su vez por el m´etodo de limitaci´on de flujo MSOU 3D expuesto en el cap´ıtulo 6. 7.Adem´as, ser´atambi´en interesante en el futuro poder hacer ejercicios de comparaci´on de metodolog´ıas. Primero, se podr´ıa hacer la comparaci´on entre las dos versiones de diferencias finitas; la primera corrige los coeficientes de difusi´on y la segunda limita los flujos advectivos (MSOU 3D). Y segundo, se podr´ıan comparar los resultados de las metodolog´ıas de diferencias finitas con los del modelo lagrangiano de Monte Carlo. 8. El siguiente paso en el futuro para el modelo de transporte de radion´uclidos no conservativos ser´ıa abordar simulaciones para periodos de varios a˜nos, del orden de una d´ecada e incluso ir m´as all´aeneltiempo.Estoser´ıa factible a corto plazo mejorando la eficiencia de los algoritmos para no tener problemas de memoria en las simulaciones, teniendo definitivamente que abrir el sistema al Oc´eano Atl´antico ya que para periodos de varios a˜nos su influencia empieza a no ser despreciable, sobre todo en la zona sur, habida cuenta de los tiempos medios de residencia del Mar B´altico que para una part´ıcula media oscilan entrelasdosylastresd´ecadas. 9.Ser´ıa interesante llegar a una soluci´on estacionaria de la distribuci´on de 137Cs bajo la 324 CAP´ ITULO 10 hip´otesis de mar cerrado, e incluso, compararla con la misma para el caso de sistema abierto provenientes de los valores experimentales de la literatura cient´ıfica. 10. Otra opci´on de trabajo futuro consiste en la posibilidad, si lo permite la red de estaciones de monitorizaci´on que coordina la HELCOM (Helsiki Commision; Baltic Marine Environment Protection Commission), de sacar una rama de nuestro modelo que permita hacer predicciones a˜no a a˜no como una herramienta m´as dentro de los amplios servicios de seguridad nuclear del norte de Europa. 11. En la ciencia de la modelizaci´on medioambiental y dentro de la geof´ısica aplicada se trabaja fundamentalmente en dos problemas, el directo y el inverso. El primero utiliza como base los modelos m´as fundamentales de la geof´ısica (atmosf´ericos y oceanogr´aficos) que son adaptados a su malla de c´alculo concreta y convertidos en el pilar sobre el que se sostiene un determinado modelo de transporte de materia en suspensi´on, de radion´uclidos o de cualquier materia pasiva natural o contaminante presente en el planeta. El problema inverso consiste en usar tanto informaci´on experimental como de modelos de transporte de una determinada sustancia –aqu´ı destacan notablemente los radion´uclidos– para avanzar cient´ıficamente en los modelos geof´ısicos m´as fundamentales de la meteorolog´ıa y la oceanograf´ıa. 12. El problema directo es que el que fundamentalmente ha protagonizado esta tesis doctoral, enmarcada claramente en el campo de la f´ısica aplicada, tan amplia en planteamiento y desarrollo, que no ha dejado mucho hueco para el an´alisis inverso, pero s´ı que pone los cimientos s´olidos para futuros estudios de esta naturaleza. 13. Una vez optimizado este problema, habr´ıa que, parte a parte, avanzar en el problema inverso para refinar el conocimiento de la f´ısica del Mar B´altico mediante el uso de la informaci´on radioecol´ogica disponible y de las herramientas matem´aticas desarrolladas en esta tesis doctoral. De este modo, al iniciarse una l´ınea de trabajo en el problema inverso, se cerrar´ıa de alg´un modo el c´ırculo entre estos dos itinerarios cient´ıficos, lo que nos lleva alareuni´on inicial que se mantuvo con los ocean´ografos expertos del SMHI en N¨orkoping (Suecia) donde se se˜nalaba al problema inverso como l´ınea de inter´es particular para ellos, m´as si cabe que el problema directo, con todas las ventajas y aplicaciones futuras que esta herramienta pudiera tener en eventuales episodios con material radiactivo presente en otros sistemas oce´anicos. El autor de esta memoria era entonces un joven investigador Sumario, Conclusiones y Trabajo Futuro 325 inexperto caminando pr´acticamente a solas en este terreno, y entendi´oqueeram´as natural que el problema inverso fuera una segunda fase de trabajo, es decir, el camino de vuelta; por eso decidi´o emprender el camino de ida trabajando, desde la lejan´ıa del sur de Espa˜na, en el problema directo del transporte de la contaminaci´on nuclear en el Mar B´altico. 14. Por todo lo anterior, creo que ha llegado el momento de definir claramente una l´ınea de trabajo futuro en el problema inverso para que all´ı retorne una parte del conocimiento aprendido y generado con la metodolog´ıa directa. Llegado a este punto, probablemente lo m´as productivo sea trabajar paralelamente con ambas metodolog´ıas para aprovechar realimentaciones. 15. Estas conclusiones vienen de nuevo a se˜nalarnos la conveniencia pr´acticamente ineludible de fortalecer o crear las colaboraciones y alianzas necesarias con las instituciones y departamentos dedicados total o parcialmente al estudio del Mar B´altico u otros sistemas naturales con problemas afines pendientes de ser abordados, cuya resoluci´on beneficie a todas las partes implicadas. Bibliograf´ıa [1] Aarup T., 2002. Transparency of North Sea and the Baltic Sea –a Secchi depth data mining study–. Oceanologia 44(3), 2002, pp.323–337. [2] Abascal A.J., Castanedo S., Medina R., Losada I.J., and Alvarez–Fanjul E.; 2009. Application of HF radar currents to oill spill modelling. Marine Pollution Bulletin, vol. 8, pp. 238–248. 2009. [3] Abril J.M., Hernandez, Garcia–Leon M., 1991. A Mathematical Approach for Modelling Radionuclide Dispersion Along the Marine Environment. Journal of Environmental Radioactivity, Vol. 13(1), pp. 39–54, 1991. [4] Abril J.M. and M. Garcia–Leon, 1992. A Marine Dispersion Model for Rasdionuclides and its Calibration from Non–radiological Information. Journal of Environmental Radioactivity, Vol. 16(2), pp. 127–46, 1992. [5] Abril J.M. and M. Garcia–Leon, 1993. A 2D 4–Phases marine dispersion model for non–conservative radionuclides. Part 1: Conceptual and computational model. Journal of Environmental Radioactivity, Vol. 20(2), pp. 71–88, 1993. [6] Abril J.M. and M. Garcia–Leon, 1993. A 2D 4–Phases marine dispersion model for non–conservative radionuclides. Part 2: Applications. Journal of Environmental Radioactivity, Vol. 20(2), pp. 89–115, 1993. [7] Abril, J. M. and Fraga, E.; 1996. Some physical and chemical features of the variability of kd distribution coefficients for radionuclides. Journal of Environmental Radioactivity 30, 253–270, 1996. 326 Bibliograf´ıa 327 [8] Abril, J. M., 1998. Basic microscopic theory of the distribution, transfer and uptake kinetics of dissolved radionuclides by suspended particulate matter. Part I: theory development. Journal of Environmental Radioactivity 41, 307–324, 1998. [9] Abril, J. M., 1998. Basic microscopic theory of the distribution, transfer and uptake kinetics of dissolved radionuclides by suspended particulate matter. Part II: applications. Journal of Environmental Radioactivity 41, 325–342, 1998. [10] Abril J.M., Abdel–Aal M.M., 2000. Marine Radioactivity Studies in the Suez Canal, Part I: Hydrodynamics and Transit Times. Estuarine, Coastal and Shelf Science, 50, pp. 489–502, 2000. [11] Abril J.M., Abdel–Aal M.M., Al–Gamal S.A., Abdel–Hay F.A., Zahar H.M.; 2000. Marine Radioactivity Studies in the Suez Canal, Part II: Field Experiments and a Modelling Study of Dispersion. Estuarine, Coastal and Shelf Science, 50, pp. 503–514, 2000. [12] Albergel A., Martin D., Strauss B. and Gros J.M.; 1988. The chernobyl accident: Modelling of dispersion over Europe of the radioactive plume and comparison with air activity measurements. Atmospheric Environment Volume 22, Issue 11, 2431–2444, 1988. [13] Baretta–Bekker J.G., Duursma E.K., Kuipers B.R.; 1998. Encyclopedia of Marine Sciences. Springer–Verlag Berlin Heidelberg, 1998. [14] Barrero, A., P´erez–Saborid M., 2005. Fundamentos y Aplicaciones de la Mec´anica de Fluidos. McGraw–Hill. [15] Bergstr¨om H., 1999 Offshore wind measurements for wind energy aplpications”. Proceedings of the European Wind Energy Applications Conference. Nice, France, pp. 1212–1215. [16] Bobertz B. and Harff J., 2004. Sediment facies and hydrodynamica setting: a study of the Baltic Sea. Ocean Dynamics, 54, 39–48. [17] Bryan, K., 1969. A numerical method for the study of ocean circulation. J. Comput. Phys., 4, pp. 400–414. [18] Chechko V.A., 2002. Analysis of space and time variations in the suspended matter distribution in Kaliningrad Bay of the Baltic Sea. Water Resources 29, 4, pp. 388–395, 2002. 328 ANEXO [19] Chen C.S., Reid R.O., Nowlin W.D., 1996. Near–inertial oscillations over the Texas Louisiana shelf. Journal of Geophysical Research–Oceans 101 (C2), 3509–3524, 1996. [20] Clark, M.J., Webb, G.A.M., 1980. A model to assess exposure from releases of radioactivity into the seas of northern Europe. In: Impacts of radionuclide releases into the marine environment. Proc. IAEA Symp., Vienna, 6–10 October 1980, pp. 629–48. [21] Colling A., Brown E., Park D., Phillips J., Rothery D., Wright J., 2004. Ocean Circulation. The Open University, Butterworth–Heinemann, 2004. [22] Csanady G.T., 1975. Hydrodymamics of large lakes. Annual Review of Fluid Mechanics, Vol. 72, No. 16, pp. 357–386. [23] Csanady G.T., 1982. Circulation in the Coastal Ocean. D. Reidel, Dordrecht, Holland. [24] Davis R.E., De Szoeke R., and Niiler P., 1981. Variability in the upper ocean during MILE. Part II: Modeling the mixed layer response. Deep–Sea Research 28A (12): 1453–1475. [25] D’Asaro E.A., 1989. The decay of wind–forced mixed layer inertial oscillations due to the βeffect. J. Geophys. Res., 94 (C2), 2045–2056. [26] Defant A.I., 1961. Physical Oceanography, vol. 1, Pergamon Press, Oxford, 729 pp. [27] Douglas R. Hazell a, Matthew H. England, 2003. Prediction of the fate of radioactive material in the South Pacific Ocean using a global high resolution ocean model. Journal of Environmental Radioactivity 65, 329–355, 2003. [28] Eilola K., Stigebrandt A., 1998. Spreading of juvenile freshwater in the Baltic Proper. J. Geophys. Res., 103 (C12), 27795–27807. [29] Eisma D., 1993. Suspended Matter in the Aquatic Environment. Springer–Verlag Berlin Heidelberg, 1993. [30] Ekman V.W., 1905. On the influence of the Earth’s rotation on ocean currents. Arkiv for Matematik, Astronomi, och Fysik, 2 (11). [31] Elliot A.J., Clarke S., 1998. Shallow water tides in the Firth of Forth. Hydrographic Journal, 87, 19–24. Bibliograf´ıa 329 [32] Emeis K., Christiansen C., Edelvang K., et. al., 2002. Material transport from the near shore to the basinal environmental in the southern Baltic Sea. II: Synthesis of data on origin and properties of material. Journal of Marine Systems 35, 2002, pp.151–168. [33] Evans, S., 1985. A box model for calculation of collective dose commitment from radioactive waterborne releases to the Baltic Sea. Journal of Environmental Radioactivity 2, 41–57. [34] Falkenmark, M. ; Mikulski, Z.; 1975. The Baltic Sea: A semi–enclosed sea, as seen by the hydrologist. Nordic Hydrology, 6, 115–136, 1975. [35] Falkenmark, M. ; Mikulski, Z; 1988. The Baltic Sea as an example: hydrological conditions as a determinant of the ecosystem in a semi–enclosed. Nature and Resources, 24, pp.14–25, 1988. [36] Feistel R., Naush G., Wasmund N., 2008. State and evolution of the Baltic Sea, 1952–2005. John Wiley and Sons, Hobolken, NJ, 704 pp. [37] Fennel W., Neumann T., 2006. Introduccion to the Modelling of Marine Ecosystems. Elsevier Oceanography Series, 72; 2006. [38] Ferrari G.M., Bo F.G. and Babin M., 2003. Geo–Chemical and optical charaterizations of suspended matter in European coastal waters. Estuarine, coastal and shelf science 57, pp.17–24, 2003. [39] Forsythe G.E. , Malcom M.A., Moler C.N., 1977. Computer Methods for Mathematical Computations, Prentice Hall, USA. [40] Franke, E., 1981. A contribution to the investigations of the current conditions in the surface layer in the area of Darss sill. Proceedings of the 12th Baltic Oceanography Conference, Leningrad. [41] Freden C., 1967. A historical review of the Ancylus Lake and the Svea River. Geol. F¨oren. Stockholm F¨orh., 89, pp. 239–267. [42] Frischmuth K., Jankowski A., 1996. Estimation of period of free oscillations (seiches) in the Kirrbucht. Oceanologia, No. 38 (4), 503–526. 336 ANEXO [110] Lehmann, A., 1995. A three–dimensional baroclinic eddy–resolving model of the Baltic Sea. Tellus 47A, 1013–1031. [111] Lehmann, A., Krauss, W., Hinrichsen, H.–H., 2002. Effects of remote and local atmospheric forcing on circulation and upwelling in the Baltic Sea. Tellus A, 54, 299–316. [112] Leonard, B. P., 1988. Simple High–Accuracy Resolution Program for Convective Modelling of Discontinuities. International Journal for Numerical Methods in Fluids, vol. 8, 1291–1318, 1988. [113] Leonard, B. P., 1991. The ULTIMATE conservative difference scheme applied to unsteady one–dimensional advection. Computer Methods in Applied Mechanics and Engineering, 88, 17–74, 1991. [114] Lepp¨aranta, M., Zhang, Z.H., 1992. A viscous–plastic ice dynamics test model for the Baltic Sea. Internal report, Vol. 3, 14 pp., Finnish Institute of Marine Research. [115] Lepp¨aranta, M., Myberg, K., 2009. Physical Oceanography of the Baltic Sea. Springer–Praxis, 2009. [116] Li˜n´an A., Clavin P., 1987. Premixed flames with non–branching chain reactions. Combustion and Flame, 70, 137–159, 1987. [117] Liu J.T., Chao S., Hsu R.T., 2002. Numerical modeling study of sediment dispersal by a river plume. Continental Shelf Research, Vol.22, Issues 11–13, pp.1745–1773, 2002. [118] Liu W.C., Hsu M.H., Kuo A.Y., 2002. Modelling of hydrodynamics and cohesive sediment transport in Tanshui River estuarine system, Taiwan. Marine Pollution Bulletin, 44(10), pp. 1076–88, 2002. [119] Lumborg U., Windelin A., 2003. Hydrography and cohesive sediment modelling: application to the Romo Dyb tidal area. Journal of Marine Systems, Volume 38, N.3, pp. 287–303, 2003. [120] Mattila J., Kankaanpaa H., Ilus E., 2006. Estimation of recent sediment accumulation rates in the Baltic Sea using artificial radionuclides Cs–137 and Pu–239,240 as time markers. Boreal Environmental Research 11 (2), 95–107, 2006. Bibliograf´ıa 337 [121] McNally G.J., Patzert W.C., Kirwan A.D., and Vastano A.C., 1983. The near– surface circulation of the North Pacific using satellite tracked drifting buoys. Journal of Geophysical Research 88 (C12): 7507–7518. [122] Meier, H. E. M., 2001. On the parameterization of mixing in three–dimensional Baltic Sea models. Journal of Geophysical Research, 106(C12), 30997–31016. [123] Meier, H. E. M., 2002. Regional ocean climate simulations with a 3D ice–ocean model for the Baltic Sea. Part 1: Model experiments and results for temperature and salinity. Climate Dynamics, 19, 237–253. [124] Meier, H. E. M., 2002b. Regional ocean climate simulations with a 3D ice–ocean model for the Baltic Sea. Part 2: Results for sea ice. Climate Dynamics, 19, 255–266. [125] Mikulski Z., 1970. Inflow of river water to the Baltic Sea in the period 1951–1960. Nordic Hydrology 4, pp.216–227, 1970. [126] National Academy of Sciences, 1972. Numerical models of ocean circulation. Proc. Symp. Durham, N.H., 1972. [127] Nerheim S., 2004. Shear–generating motions at various length scales and frequencies in the Baltic Sea – an attempt to narrow down the problem of horizontal dispersion. Oceanologia, 46 (4), 2004. pp. 477–503. [128] Nakano M., Povinec P.P., 2003. Oceanic general circulation model for the assessment of the distribution of 137Cs in the world ocean. Deep Sea Research, Part II, 50, pp. 2803–2816. [129] Nicholson J., O’Connor B.A., 1986. Cohesive Sediment Transport Model. Journal of Hydraulic Engineering, Vol. 112, 7, pp.621–640, 1986. [130] Nielsen S.P., 1995. A box model for North–East Atlantic coastal waters compared with radioactive tracers. Journal of Marine Systems, 6, 545–560, 1995. [131] Nies, H., 1989. The distribution of Chernobyl fallout over the Baltic Sea and its change during 1987 and 1988. In: Three years observations of the levels of some radionuclides in the Baltic Sea after the Chernobyl accident. Baltic Sea environment proceedings, vol. 31, Helsinki, Finland, pp. 31–51. 338 ANEXO [132] Oh I.S., Zhurbas V., Park W., 2000. Estimating horizontal diffusivity in the East Sea (Sea of Japan) and the northwest Pacific from satellite–tracked drifter data. Journal Of Geophysical Research, vol. 105, No. C3, pp. 6483–6492, 2000. [133] Okubo A., 1971 Oceanic diffusion diagrams. Deep–Sea Res., 18, 789–802. [134] Omstedt, A., 1990. Modelling the Baltic Sea as thirteen sub–basins with vertical resolution. Tellus 42A, 286–301. [135] Omstedt,A., Nyberg, L., 1991. Sea level variations during ice–covered periods in the Baltic Sea. Geophysica, 27(1–2), 41–61. [136] Omstedt, A.; Nyberg L., 1996. Response of Baltic Sea ice to seasonal, interannual forcing and climate change. Tellus 48A, 644–662. [137] Omstedt, A., Axell, L.B., 1998. Modeling the seasonal, interannual, and long–term variations of salinity andtemperature in the Baltic Proper. Tellus 50A, 637–652. [138] Omstedt A., Elken J., Lehmann A., Piechura J., 2004. Knowledge of the Baltic Sea physics gained during the BALTEX and related programmes. Progress in Oceanography, 63 (1–2), 1–28, 2004. [139] Osinski R., Rak D., Walczowski W., Piechura J., 2010. Baroclinic Rossby radius of deformation in the southern Baltic Sea. Oceanologia, 52 (3), 2010. pp. 417–429, 2010. [140] Peria˜nez P., Abril J.M., Garcia–Leon M.; 1994. A Modeling Study of Ra–226 Dispersion in an Estuarine System in South–West Spain. Journal of Environmental Radioactivity, vol. 24, no. 2, pp. 159–179; 1994. [141] Peria˜nez, R. and Elliott, A.J., 2002. A particle tracking method for simulating the dispersion of non–conservative radionuclides in coastal waters. Journal of Environmental Radioactivity 58, 13–33. [142] Peria˜nez R., 2004. A particle–tracking model for simulating pollutant dispersion in the Strait of Gibraltar. Marine Pollution Bulletin 49, 613–623. [143] Peria˜nez R., 2005. Modelling the Dispersion of Radionuclides in the Marine Environment. An Introduction. Springer–Verlag Berlin Heidelberg, 2005. Bibliograf´ıa 339 [144] Peria˜nez R., 2005b. GISPART: a numerical model to simulate the dispersion of contaminants in the Strait of Gibraltar. Environmental Modelling and Software 20, 797–802. [145] Peria˜nez R., 2005c. An operative lagrangian model for simulating radioactivity dispersion in the Strait of Gibraltar. Journal of Environmental Radioactivity 84, 95–101. [146] Peri´a˜nez R. and Pascual–Granged A., 2008. Modelling surface radioactive, chemical and oil spills in the Strait of Gibraltar. Computers and Geosciences 34, 163–180. [147] Peri´a˜nez R., 2010. Fundamentos de Oceanograf´ıa Din´amica. Secretariado de Publicaciones, Universidad de Sevilla, 2010. [148] Persson, C., Rodhe, H. De Geer, L. E.; 1987. The Chernobyl accident. A meteorological analysis of how radionuclides reached and were deposited in Sweden. Ambio 16, 20–31, 1987. [149] Prandle D., 1984. Phil. Trans. R. Soc. Lond. A. 310, 407–36. [150] PrandleD.,HargreavesJ.C.,McManusJ.P.,CampbellA.R.,DuweK.,Lane A., Mahnke P., Shimwell S., Wolf J., 2000. Tide, wave and suspended sediment modelling on an open coast, Holderness. Coastal Engineering, Vol. 41, N.1, pp. 237–267, 2000. [151] Preller R.H., Cheng A., 1999. Modeling the Transport of Radioactive Contaminants in the Arctic. Marine Pollution Bulletin, Vol. 38, No. 2, pp. 71–91, 1999. [152] Puls W. and Sundermann J., 1990. Simulation of suspended sediment dispersion in the North Sea. In: Residual currents and long term transport. pp.356–372. [153] Radziejewska T., Chabior M., 2004. Climatic and hydrological controls over the zoobenthos in a southern Baltic coastal lagoon. Hydrobiologia, 514, pp.171–181, 2004. [154] Ralph E.A., Niiler P.P., 2000. Wind–driven currents in the tropical Pacific. Journal of Physical Oceanography 29 (9): 2121–2129. [155] Ramming H.G., Kowalik Z., 1980. Numerical modelling of marine hydrodynamics, Elsevier Oceanography Series, 26, 1980, 369 pp. 340 ANEXO [156] Reynolds O., 1883. An experimental investigation of the circumstances which determine whether the motion of water will be direct of sinuous, and the law of resistance in parallel channels. Philosophical Transactions, Royal Society London 174: 935. [157] Ribbe, J., M¨uller–Navarra, S.H., Nies, H., 1991. A one–dimensional dispersion model for radionuclides in the marine environment appliedto the Chernobyl fallout over the Northern Baltic Sea. Journal of Environmental Radioactivity 14, 55–72. [158] Richardson E.G., 1961. 2nd edition Dynamics of Real Fluids. London: Edwards Arnolds. [159] Ruiz–Segura J., Losada–Rodr´ıguez M.A., Polo–G´omez M.J., Bramato S., D´ıez–Minguito, M., 2010. Propuesta metodol´ogica para diagnosticar y pronosticar las consecuencias de las actuaciones humanas en el Estuario del Guadalquivir. 2010. [160] Sanchez–Cabeza J.A., Ortega M., Merino J., Masque P.; 2002. Long–term box modelling of 137Cs in the Mediterranean Sea. Journal of Marine Systems 33–34, 457– 472, 2002. [161] Sankaranarayanan S., Shankan N.J., Cheong H.F., 1998. Three–Dimensional Finite Difference Model for Transport of Conservative Pollutants. Ocean Engineering, 25, 6, 425–442, 1998. [162] Savchuk O. and Rodriguez Medina M., 2003. Secchi depth and nutrient concentrations in the Baltic Sea. Technical Report, Dept. of Systems Ecology, Stockholm University. [163] Schrum, C., Hubner, U., Jacob, D. and Podzun, R., 2003. A coupled atmosphere– ice–ocean model for the North Sea and the Baltic Sea. Climate Dynamics, 21, 131–151. [164] Scott E.M., Gurbuttb P., Harms I., Heling R., Nielsen S.P., Osvat I., Preller R., Sazykina T., Wada A., Sjoeblom K.L; 1997. Benchmarking of numerical models describing the dispersion of radionuclides in the Arctic Seas. The Science of the Total Environment, 202, 123–134, 1997. [165] Skamarock W.C., 2006. Positive–Definite and Monotonic Limiters for Unrestricted– Time–Step Transport Schemes. Monthly Weather Review, 134, 2241–2250, 2006. Bibliograf´ıa 341 [166] Shearman R.K., 2005. Observations of near–inertial current variability on the New England shelf. Journal of Geophysical Researc–Oceans 110 (C2), Art. No. C02012, 2005. [167] Sheu W.H., Lee S–M, Wang M.T., 1997. Multi–dimensional monotone flux discretization scheme for convection dominated flows. International Journal for Numerical Methods for Heat and Fluid Flow, 7, 8, pp. 814–842, 1997. [168] Simons, T.J., 1973. Development of three–dimensional numerical model of the Great Lakes. Inland Water Directoriate Scientific Series No.12, Canada Centre for Inland Waters, Burlington, Ontario, 1973, 26pp. [169] Simons, T.J., 1974. Verification of numerical models of Lake Ontario. Part 1: Circulation in spring and early summer. J. Phys. Oceanogr., 4, 507–523. [170] Simons, T.J., 1978. Wind–driven circulation in sotheast Baltic. Tellus, 30, 272–283. [171] SMHI (Swedish Meteorological and Hydrological Institute), 2010. http://www.smhi.se/en/services/professional–services/data–and–statistics/wind–statistics– 1.7664. Wind Statistics Service, SMHI. [172] Smolarkiewicz P.K., 1984. A Fully Multidimensional Positive Definite Advection Transport Algorithm with small Implicit Diffusion. Journal of Computational Physics, 54, 325–362, 1984. [173] Soomere T., Keevallik S., 2001. Anisotropy Of Moderate And Strong Winds In The Baltic Proper. Proc. Estonian Acad. Sci. Eng., 7, 1, 35–49; 2001. [174] Stanev E.V., Buesseler K.O., Staneva J.V., Livingston H.D., 1999. Acomparison of modelled and measured Chernobyl 90Sr distributions in the Black Sea. Journal of Environmental Radioactivity, 43(2), pp. 187–203. [175] Stewart R.H., 2006. Introduction To Physical Oceanography. Texas A &M University. [176] Stigebrandt, A., 1983. A model for the exchange of water and salt between the Baltic Skagerrak. Journal of Physical Oceanography 13, 411–427. [177] Stigebrandt A., Lass H.U., Liljebladh B., Alenius P., Piechura J., Hietala R., Beszczynska A., 2002. DIAMIX – An experimental study of diapycnal deepwater mixing in the virtually tideless Baltic Sea. Boreal Environment Research, 7 (4), 363–369, 2002. 342 ANEXO [178] Stokes G.G., 1847. On the theory of oscillatory waves. Transactions of the Cambridge Philosophical Society, vol. 8, part 4, 441–455, 1847. [179] Svansson A., 1972. Canal models of sea level and salinity variations in the Baltic and adjacents waters, Series Hydrography, Rap. 26, Fishery Board of Sweden, 72 pp. [180] Svansson, A., 1972b. The water exchange of the Baltic. Ambio Special Report 1, 15–19. [181] Svansson, A., 1980. Exchange of water and salt in the Baltic and adjacent seas. Oceanologica Acta 3, 431–440. [182] Taylor, G.I., 1921 Diffusion by continuous movements. Proceedings of the London Mathematical Society (2), 20, 196. [183] Toscano–Jimenez M., R. Garcia–Tenorio, 2001. Proceedings of I International Congress On Marine Science And Technology, Oceans III Millennium. Pontevedra, 2001. [184] Toscano–Jimenez M., Garcia–Tenorio R., 2002. Proceedings of 5th International Topical Meeting on Industrial Radiation and Radioisotope Measurement Applications, Bologna, Italy, 2002. [185] Toscano–Jimenez M., Garc´ıa–Tenorio R., 2002b. Proceedings of Computational and Mathematical Methods on Science and Engineering (CMMSE–2002 Conference), Alicante, Spain, 2002. Vol. II, p. 342. 2002. [186] Toscano–Jimenez M., Garcia–Tenorio R., 2002c. Proceedings of OCEANS 2002 MTS/IEEE, Biloxi (Mississippi, USA), pp. 2472–2476. 2002. [187] Toscano–Jimenez M. and R. Garc´ıa–Tenorio, 2003. Modelling The Dispersion Of 137Cs In Marine Ecosystems With Monte Carlo Methods. Nuclear Instruments and Methods B, vol.213, pp. 789–793. [188] Toscano–Jimenez M., Garcia–Tenorio R. and Abril J.M., 2003b. Proceedings of the First International Meeting on Applied Physics, Badajoz, Spain, 2003. [189] Toscano–Jimenez M. and Garc´ıa–Tenorio R., 2004. A three–dimensional model for the dispersion of radioactive substances in marine ecosystems. Application to the Baltic Sea after the Chernobyl disaster. Ocean Engineering 31, pp. 999–1018. 2004 Bibliograf´ıa 343 [190] Toscano–Jimenez M., Abril J.M. and Garcia–Tenorio R., 2004b. Proceedings of OCEANS 2004 MTS/IEEE Conference. Kobe (Japan), pp. 2362–2365, 2004. [191] Toscano–Jimenez M., Abril J.M. and Garcia–Tenorio R., 2005. Proceedings of OCEANS 2005, MTS/IEEE Conference, Brest (France), pp. 475–480, 2005. [192] Toscano–Jimenez M., Abril J.M., 2007 Proceedings of OCEANS 2007 MTS/IEEE Conference, Aberdeen (Scotland, U.K.), pp. 521–526, 2007. [193] Toscano–Jimenez M., Abril J.M., 2010. Modelling the transport of particles released from nuclear accidents. Validation in the Baltic Sea after Chernobyl disaster. Abstract Book, International Conference on Engineering Failures Analysis, University of Cambridge, pp. 28–29, 2010. [194] Toscano–Jimenez M., Abril J.M., 2011. A Lagrangian Transport Model Applied to two Different Brackish Systems: the Baltic Sea and the Guadalquivir River. Proceedings of the OCEANS 2011 IEEE Conference, Santander (Spain), pp.1–8, 2011. [195] Tsanis I.K., Wu J., 2000. Application and verification of a three–dimensional hydrodynamic model to Hamilton Harbour, Canada. Global Nest: the International Journal, Vol. 2, No. 1, pp 77–89. [196] Udin I., Omsted A., 1976. Sea Ice–75, Dynamical Report. Winter Navigation Research. 16(8), 1–64. [197] Van Meurs P., 1998. Interactions between near–inertial mixed layer currents and the mesoscale: The importance of spatial variability in the vorticity field. Journal of Physical Oceanography 28 (7): 1363–1388. [198] Valli A., Lepp¨aranta M., 1975. Calculations of ice drift in the Bothnian Bay and the Quark. Winter Navigation Research. 13, 1–14. [199] Vested H.J., Baretta J.W., Ekebjaerg L.C., Labrosse A., 1996. Coupling of hydrodynamical transport and ecological models for 2D horizontal flow. Journal of Marine Systems, 8, 255–267, 1996. [200] Voipio, A., 1981. TheBalticSea.Elsevier, Amsterdam, pp. 162–167. 344 ANEXO [201] Weller R.A., Plueddmann A.J., 1996. Observations of the vertical structure of the oceanic boundary layer. Journal of Geophysical Research 101 (C4): 8, 789–806. [202] Weiss, D., 1989. The distribution of radionuclides in bottom sediments of the Open Baltic Sea and Greifswald Bodden. In: Three years observations of levels of some radionuclides in the Baltic Sea after the Chernobyl accident. Baltic Sea environment proceedings, vol. 31, Helsinki, Finland, pp. 94–122. [203] Witting R., 1911. Tidvatten i Osterj¨on och Finska viken. Fennia, 29(2), 84 pp. [204] Wright J., Colling A., Park D., Brown E., Phillips J., Rothery D., 2005. Waves, Tides and Shallow–Water Processes. The Open University, Butterworth–Heinemann, 2005. [205] Wu J., 1975. Wind–induced drift currents. Journal of Fluid Mechanics, Vol. 68, No. 1, pp.49–70. [206] Wu J., 1993. An Integrated Hydrodynamic and Pollutant Transport Model for the Nearshore Areas of the Great Lakes and Their Tributaries. Open Access Dissertations and Theses. pp.1–336, Paper 3117, McMaster University, Ontario, Canada. [207] Wu J., and Tsanis I.K., 1995. VHI3D: A vertical/horizontal integrating wind–induced circulation model. Advances in Water Research, Vol. 18, No. 2, pp.77–87. [208] Wu Y., Falconer R.A. and Uncles R.J., 1999. Modelling of Water Flows and Cohesive Sediment Fluxes in the Humber Estuary. Marine Pollution Bulletin, Vol. 37, Issues 3–7, pp. 182–189, 1999. [209] Wulff F., Rahm L., Larsson P.; 2001. A Systems Analysis of the Baltic Sea. Springer– Verlag Berlin Heidelberg, 2001. [210] Zhang, Z., Lepp¨aranta, M., 1995. Modelling the influence of ice on sea level variations in the Baltic Sea. Geophysica, 31(2), 31–45. [211] Zhurbas V.M., Laanemets J., Kuzmina N.P., Muraviev S.S., Elken J., 2008. Direct Estimates of the Lateral Eddy Diffusivity in the Gulf of Finland of the Baltic Sea (Based on the Results of Numerical Experiments with an Eddy Resolving Model). Oceanology, 2008, Vol. 48, No. 2, pp. 175–181, 2008. Bibliograf´ıa 345 [212] Zhurbas V., Laanemets J., Vahtera E., 2008 b. Modeling of the mesoscale structure of coupled upwelling/downwelling events and the related input of nutrients to the upper mixed layer in the Gulf of Finland, Baltic Sea. Journal of Geophysical Research, vol. 113, C05004, 2008. [213] Zhurbas, V. M., Elken, J., V¨ali, G., Kuzmina, N.P., Paka, V.T.; 2010. Pathways of suspended particles transport in the bottom layer of the southern Baltic Sea depending on the wind forcing (Numerical Simulation). Oceanology, Volume 50, Issue 6, pp.841–854, 2010. [214] Ziervogel K., and Bohling B., 2003. Sedimentological parameters and erosion behaviour of submarine coastal sediments in the south–western Baltic Sea. Geo–Marine Letters, Vol. 23, No. 1, pp. 43–52, 2003. [215] Zubov V.A., Rozanov E.V., Schlesinger M.E., 1999. Hybrid Scheme for Three– Dimensional Advective Transport. Monthly Weather Review, 127, 1335–1346, 1999.