Un modelo 1D NPZ de acoplamiento entre la hidrodinámica y los flujos biogeoquímicos en estrechos bicapa. Aplicación a la dinámica mareal en el Estrecho de Gibraltar

XXIII Congreso de Ecuaciones Diferenciales y Aplicaciones XIII Congreso de Matem´ atica Aplicada Castell´on, 9-13 septiembre 2013 (pp. 1–8) Un modelo 1D NPZ de acoplamiento entre la hidrodin´amica y los flujos biogeoqu´ımicos en estrechos bicapa. Aplicaci´on a la din´amica mareal en el Estrecho de Gibraltar J. Mac´ ıas1, M.J. Castro1, E. Fern´ andez-Nieto2 1Universidad de M´alaga, Departamento de An´alisis Matem´atico, Facultad de Ciencias, Campus de Teatinos, s/n, 29080-M´alaga. E-mails: [email protected], [email protected]. 2Universidad de Sevilla Departamento de Matem´atica Aplicada I, ETS Arquitectura, Avda. Reina Mercedes, 2, 41012-Sevilla. E-mail: [email protected]. Palabras clave: Modelo NPZ, shallow water bicapa, acoplamiento, flujos biogeoqu´ımicos, Estrecho de Gibraltar Resumen Los modelos NPZ (siglas para Nutrientes-Fitoplancton-Zooplancton en ingl´es) son com´unmente utilizados en estudios de biolog´ıa marina. Este tipo de modelos utiliza un conjunto de ecuaciones diferenciales muy sencillo para definir la din´amica del plancton marino. Las variables de estado de las cuales se modeliza su evoluci´on son los nutrientes, el fitoplancton y el zooplankton, esto se hace en t´erminos de su contenido de nitr´ogeno, ya que este compuesto el que normalmente limita la producci´on primaria en el oc´eano. En este trabajo se implementa el acoplado de un modelo NPZ para el modelado de los flujos biogeoqu´ımicos con un modelo de aguas poco profundas bicapa para la hidrodin´amica. El objetivo es aplicar este modelo a la simulaci´on de flujos biogeoqu´ımicos en el Estrecho de Gibraltar con imposici´on de la din´amica mareal. 1. Introducci´on El Estrecho de Gibraltar es un ´area de nuestros mares y oc´eanos fascinante desde muchos puntos de vista. Su din´amica es muy complicada y energ´etica, y los procesos biol´olgicos asociados de un enorme inter´es socio-econ´omico y ambiental. En la actualidad no existe ning´un modelo f´ısico-biol´ogico completamente acoplado para el Estrecho de Gibraltar. Un primer intento se recoge en [5] donde en primer lugar se integra el modelo hidrodin´amico y tras la simulaci´on se extraen los valores de las variables necesarias para la integraci´on del modelo biol´ogico. En este trabajo proponemos un modelo unidimensional 1 J. Mac´ıas, M.J. Castro, E. Fern´andez-Nieto que acopla la hidrodin´amica y los procesos biol´ogicos en estrechos de secci´on arbitraria y con un flujo bicapa, se pretende as´ı representar la din´amica que se observa en el Estrecho de Gibraltar, que es esencialmente bicapa. La estrategia tanto num´erica como computacional difiere completamente del trabajo anteriormente mencionado. En [6] aplicamos este tipo de estrategias para modelos acoplados con una componente biol´ogica de mucha mayor complejidad que la que aqu´ı se presenta. 2. La componente hidrodin´amica La parte hidrodin´amica del modelo acoplado consiste en un modelo de aguas poco profundas unidimensional de dos capas inmiscibles con superficie libre en un canal de secci´on arbitraria. La deducci´on de este modelo puede encontrarse en [1]. En nuestro caso la primera capa representa la capa de agua atl´antica que penetra en superficie en el Mar de Albor´an, en la cual tienen lugar los procesos biol´ogicos, mientras que la segunda capa est´a formada por el agua mediterr´anea m´as salada y densa que fluye en profundidad en su camino a trav´es del Estrecho hacia en Atl´antico Norte. El sistema de ecuaciones diferenciales que resulta es el siguiente:                                                          ∂A1 ∂t +∂Q1 ∂x = 0, ∂Q1 ∂t +∂ ∂x Q2 1 A1 +g 2σ1 A2 1=−g σ1 A1 ∂A2 ∂x +g 21 σ1x A2 1−gdb dx σb σ1 A1+ g σ1 (I3,1+I3,2)A1, ∂A2 ∂t +∂Q2 ∂x = 0, ∂Q2 ∂t +∂ ∂x Q2 2 A2 +g 2σ2 A2 2=−rg σ1 A2 ∂A1 ∂x +g 21 σ2x A2 2−gdb dx σb σ2 A2+ gr σ1 (I3,1+I3,2) + (1 −r) σ3 I3,2A2, (1) donde el ´ındice 1 denota la capa superficial y el ´ındice 2 la capa m´as profunda. La coordenada xrepresenta la posici´on a lo largo de un eje tranversal al canal, tes el tiempo, gla aceleraci´on de la gravedad, b(x) es la funci´on que define la batimetr´ıa del fondo y σ(x, z) la funci´on anchura que depende de la coordenada vertical. Cada una de las capas se supone que tiene una densidad constante, ρi,i= 1,2 (ρ1< ρ2), siendo r=ρ1 ρ2 la raz´on entre estas densidades. Las variables Ai(x, t), i= 1,2 representan las secciones h´umedas para cada capa y se relacionan con el espesor de la capa, hi(x, t), i= 1,2 mediante las expresiones: A1(x, t) = Zb+h2+h1 b(x)+h2(x,t) σ(x, z)dz, A2(x, t) = Zb(x)+h2(x,t) b(x) σ(x, z)dz. (2) 2 Un modelo 1D de acoplamiento hidrodin´amica-biolog´ıa La inc´ognitas Qi(x, t), i= 1,2 representan el flujo de masa en cada secci´on de coordenada xen el tiempo ty se relacionan con la velocidad media de cada capa mediante la expresi´on: Qi(x, t) = ui(x, t)Ai(x, t), i = 1,2. Los t´erminos I3,1eI3,2que aparecen en las ecuaciones de conservaci´on de momento vienen dados por: I3,1(x, t) = Zb(x)+h1(x,t)+h2(x,t) b(x)+h2(x,t) ∂σ ∂x (x, z)dz, I3,2(x, t) = Zb(x)+h2(x,t) b(x) ∂σ ∂x (x, z)dz. (3) Finalmente, σ1,σ2,σ3yσbest´an definidas por σ1(x, t) = σ(x, b(x) + h2(x, t) + h1(x, t)), σ3(x, t) = σ(x, b(x) + h2(x, t)), σ2(x, t) = σ1σ3 σ1(1 −r) + rσ3 , σb(x) = σ(x, b(x)). El sistema (1) se puede escribir como: ∂wH ∂t +A(σ, wH)∂wH ∂x =S(x, σ, wH, b) (4) donde wH=A1(x, t)Q1(x, t)A2(x, t)Q2(x, t)T, σ=σ1(x, t)σ2(x, t)σ3(x, t)σb(x)T, A(σ, wH) =        0100 g¯ h1−u2 12u1g¯ h10 0001 rgˆ h20g¯ h2−u2 22u2        (5) con ¯ hi=Ai σi , i = 1,2ˆ h2=A2 σ1 (6) y S(x, σ, wH, b) =            0 gA1 σ1Zb+h1+h2 b ∂σ ∂x dz −gA1 σb σ1 db dx 0 −gσb σ2 A2 db dx +gA2r σ1Zb+h1+h2 b ∂σ ∂x dz +1−r σ3Zb+h2 b ∂σ ∂x             .(7) Teniendo en cuenta las definiciones de A1yA2y usando la regla de Leibnitz se tiene que: ∂ ∂x(A1+A2) = Zb+h1+h2 b ∂σ ∂x dz +db dx +∂h1 ∂x +∂h2 ∂x σ1−db dxσb(8) 3 J. Mac´ıas, M.J. Castro, E. Fern´andez-Nieto ∂A2 ∂x =Zb+h2 b ∂σ ∂x dz +db dx +∂h2 ∂x σ3−db dxσb.(9) Ahora, usando (8) y (9), el t´ermino fuente S(x, σ, w, b) se puede reescribir formalmente como: S(x, σ, wH, b) =           0 g¯ h1 ∂ ∂x (A1+A2)−gA1 ∂ ∂x (b+h1+h2) 0 g¯ h2 ∂A2 ∂x +rgˆ h2 ∂A1 ∂x −gA2 ∂ ∂x (b+h2+rh1)           .(10) Los autovalores de A(σ, w) se clasifican en externos (2) e internos (otros 2). Los autovalores externos, λ± ext, est´an relacionados con la velocidad de propagaci´on de las perturbaciones barotr´opicas y los autovalores internos, λ± int, con la velocidad de propagaci´on de la perturbaciones barocl´ınicas. Desafortunadamente no existen expresiones anal´ıticas para estos autovalores, sin embargo, cuando r∼ =1, se tienen las siguientes aproximaciones de primer orden (v´ease [7]): λ± ext =u1¯ h1+u2¯ h2 ¯ h1+¯ h2 ±qg(¯ h1+¯ h2) (11) λ± int =u1¯ h2+u2¯ h1 ¯ h1+¯ h2 ±sg0¯ h1¯ h2 ¯ h1+¯ h21−(u1−u2)2 g0(¯ h1+¯ h2),(12) con g0=g(1 −rˆ h2 ¯ h2). Obs´ervese como, de acuerdo con la expresi´on (12), el sistema (1) es hiperb´olico si F2 I=(u1−u2)2 g0(¯ h1+¯ h2)≤1, siendo F2 Iconocido como n´umero de estabilidad de Froude. En los casos que abordamos en las aplicaciones reales, donde r∼ =1, se tienen las siguientes desigualdades: λ− ext < λ− int < λ+ int < λ+ ext. Dependiendo del valor de los autovalores, el flujo se dice que es subcr´ıtico si los autovalores internos tienen signos distintos, cr´ıtico si uno de estos autovalores es cero y, en otro caso, el flujo se dice supercr´ıtico. De la expresi´on de la matriz A(σ, w) se deduce que las secciones cr´ıtica, subcr´ıtica y supercr´ıtica se caracterizan por ser secciones donde G2= 1, G2>1, G2<1, respectivamente, donde G2=F2 1+F2 2−g0 gF2 1F2 2,(13) y F2 1=u2 1 g0¯ h1 , F2 2=u2 2 g0¯ h2 .(14) En estas expresiones GyFi,i= 1,2 son las definiciones adecuadas para este caso del n´umero de Froude compuesto y el n´umero de Froude interno, respectivamente. 4 Un modelo 1D de acoplamiento hidrodin´amica-biolog´ıa 3. La componente biol´ogica: un modelo NPZ La componente biol´ogica del modelo acoplado es un modelo sencillo NPZ que, b´asicamente, es el mismo que el descrito en [5], pero formulado en variables conservadas. Las tres variables de estado del modelos, N,P,Z, se modelizan usualmente en funci´on de su contenido de nitr´ogeno, al ser el nitr´ogeno el componente que normalmente limita la producci´on primaria. Los cambios temporales en la concentraci´on de los contituyentes vienen dadas por las ecuaciones: ∂(A1N) ∂t +∂(u1A1N) ∂x =−f(I)g(N)A1P+ (1 −β)h(P)A1Z+γ j(Z)A1Z. ∂(A1P) ∂t +∂(u1A1P) ∂x =f(I)g(N)A1P−h(P)A1Z−i(P)A1P, ∂(A1Z) ∂t +∂(u1A1Z) ∂x =β h(P)A1Z−j(Z)A1Z, (15) El primer paso para formular un modelo NPZ consiste en determinar las expresiones para las funciones de transferencia que aparecen en las ecuaciones anteriores. La elecci´on de la forma funcional para estos t´erminos es cr´ıtica para la din´amica del modelo y puede tambi´en restringuir la parametrizaci´on del modelo debido a los coeficientes necesarios para su escritura. En este trabajo consideramos las mismas parametrizaciones para las funciones de transferencia que en [5]. Sin embargo, hemos implementado una versi´on del modelo extremadamente vers´atil que permite a cualquier potencial usuario del modelo, que lo utilice a trav´es de la Plataforma HySEA, escribir y usar cualquier otra expresi´on funcional para las ecuaciones del modelo NPZ particular que desee considerar de forma sencilla. Tambi´en es posible realizar la implementaci´on de otros tipos de modelos biol´ogicos m´as complicados, como el que se presenta en [6]. Las cinco funciones de transferencia que aparecen en un modelo NPZ han sido parametrizadas, en nuestro caso, como se indica a continuaci´on: f(I) = Jes la respuesta del fitoplancton a la luz, siendo una funci´on de la irradiaci´on solar, h1,Py otros par´ametros como latitud, nubosidad, ... g(N) = Qp=µp N k1+N,es el consumo de nutrientes por el fitoplancton donde µp es la tasa m´axima de crecimiento para el fitoplancton, y k1es el valor medio para la saturaci´on del consumo de nitrato. h(P) = Gp=µgP kg+P,parametriza la ingesti´on de fitoplancton por el zooplancton, donde µges la tasa m´axima de consumo de fitoplancton por el zooplancton, y kgel valor medio de saturaci´on. i(P) = mp=µ1P kg+P,es la p´erdida de fitoplancton, donde µ1es la tasa de mortalidad natural m´axima del fitoplancton. 5 J. Mac´ıas, M.J. Castro, E. Fern´andez-Nieto j(Z) es la p´erdida de zooplankton debido a muerte, excreci´on y predaci´on: mz=µzZ kz+Z, donde kzes la mitad del valor de saturaci´on para la mortalidad del zooplancton y µzla tasa m´axima de mortalidad del zooplancton. Finalmente, βes la asimilaci´on de zooplancton y γes la fracci´on de p´erdida de zooplancton en relaci´on a los nutrientes. 3.1. Producci´on debida a los procesos de mezcla Adem´as de la advecci´on, los procesos de mezcla entre la capa superficial de agua atl´antica y la capa profunda mediterr´anea son fundamentales para la compresi´on de los patrones biogeoqu´ımicos en este ´area marina (Mac´ıas et al., 2007). Con el objetivo de evitar el uso de parametrizaciones complejas para los procesos de mezcla interfacial, que adem´as dependen fuertemente en la elecci´on de diversos par´ametros a los que muestran gran sensibilidad, resultando dif´ıciles de ajustar, nosotros hemos considerado el n´umero de estabilidad de Froude, F2 I, para determinar cuando los procesos de mezcla tienen lugar y cuando no, a la vez que para modular su intensidad. Cuando F2 I>1 esto est´a relacionado con la aparici´on de las conocidas como inestabilidades de Kelvin-Helmholtz y, en tal caso, la mezcla juega un papel importante. As´ı pues, podremos parametrizar la mezcla mediante un coeficiente de mezcla a trav´es de la picnoclina que sea constante o bien modularlo en funci´on de este n´umero de estabilidad de Froude. Teniendo todo esto en cuenta, la expresi´on final para el modelo para los flujos biogeoqu´ımico considerado es la siguiente:                    ∂ A1N ∂t +∂(u1A1N) ∂x =kmix (Ndeep −N)−JQpA1P+ (1 −β)GpA1Z+γ mzA1Z, ∂ A1P ∂t +∂(u1A1P) ∂x =−kmix P+JQpA1P−GpA1Z−mpA1P, ∂ A1Z ∂t +∂(u1A1Z) ∂x =−kmix Z+β GpA1Z−mzA1Z, (16) donde Ndeep es la concentraci´on de nitrato en la capa inferior y kmix es el coeficiente o funci´on de mezcla a trav´es de la picnoclina. En nuestro caso, kmix es una funci´on que depende de F2 IyA1. Por ´ultimo el sistema (16) se puede escribir como: ∂wB ∂t +∂FB ∂x (wB, u1) = SB(wB),(17) donde wB=  A1N A1P A1Z , FB(wB, u1) =   u1A1N u1A1P u1A1Z  6 Un modelo 1D de acoplamiento hidrodin´amica-biolog´ıa y SB(wB) =     kmix (Ndeep −N)−JQpA1P+ (1 −β)GpA1Z+γ mzA1Z, −kmix P+JQpA1P−GpA1Z−mpA1P, −kmix Z+β GpA1Z−mzA1Z.     4. Esquema num´erico Las soluciones de (1) y (17) pueden desarrollar discontinuidades, incluso partiendo de datos iniciales suaves. En este trabajo, proponemos un m´etodo de vol´umenes finitos de segundo orden con limitador de flujo, basado en una combinaci´on no lineal adecuada del esquema IFCP (de “Intermediate Field Capturing Parabola”) introducido en [2] y el esquema Lax-Wendroff, con un tratamiento especial para las ecuaciones de transporte del sistema (17). El esquema que resulta se reduce al IFCP para la componente hidrodin´amica y a un esquema tipo upwind con un tratamiento descentrado upwind para las ecuaciones de transporte en las regiones de soluci´on no suave y al m´etodo de Lax-Wendroff en las zonas regulares, consigui´endose orden dos en las regiones suaves. Para ver una descripci´on completa del esquema num´erico v´ease [4]. 5. Tests num´ericos Hemos realizado multitud de test num´ericos para distintas geometr´ıas del canal, imponiendo una diversidad de condiciones de contorno y, en particular, imponiendo el efecto de la luz solar constante o dependiente del tiempo. Como ejemplo m´as complejo hemos considerado la geometr´ıa del Estrecho de Gibraltar (extremadamente abrupta y complicada), sobre la cual hemos hecho una imposici´on realista del efecto de las mareas, al igual que en [1], teniendo en cuenta sus cuatro componentes principales. Hemos realizado varias simulaciones de 10 d´ıas en diferentes ´epocas del a˜no y con una parametrizaci´on adecuada del t´ermino de irradiaci´on solar. Figura 1: Instant´aneas del modelo NPZ en el Estrecho de Gibraltar. 7 J. Mac´ıas, M.J. Castro, E. Fern´andez-Nieto Figura 2: Instant´aneas del modelo NPZ en el Estrecho de Gibraltar. Agradecimientos Este trabajo ha sido parcialmente financiado por el Gobierno de Espa˜na y los fondos FEDER a trav´es del proyecto de investigaci´on DAIFLUID (referencia MTM2012-38383- C02-02) y por la Junta de Andaluc´ıa a trav´es del proyecto de excelencia TESELA (referencia P11- RNM7069). Referencias [1] M.J. Castro, J.A. Garc´ıa-Rodr´ıguez, J.M. Gonz´alez-Vida, J. Mac´ıas y C. Par´es. Numerical simulation of two-layer shallow-water flows through channels with irregular geometry J. Comp. Physics., Vol. 195, (2004), pp. 202-235. [2] E.D. Fern´andez-Nieto, M.J. Castro y C. Par´es. On an intermediate field capturing Riemann solver based on a parabolic viscosity matrix for the two-layer shallow water system J. Sci. Comput., 48(1), (2011), pp. 117-140. [3] P.J. Franks. NPZ models of plankton dynamics: Their construction, coupling to physics, and application. J. Oceanogr., Vol. 58, (2002), pp. 379-387. [4] J. Mac´ıas, M.J. Castro y E.D. Fern´andez-Nieto. A 1D coupled hydrodynamical biological NPZ model for channels with arbitrary sections. En preparaci´on, 2013. [5] D. Mac´ıas, A.P. Mart´ın, J. Garc´ıa-Lafuente, C:M: Garc´ıa, A. Yool, M. Bruno, A. V´azquez-Escobar, A. Izquierdo, D.V. Sein y F. Echevarr´ıa. Analysis of mixing and biogeochemical effects induced by tides on the Atlantic-Mediterranean floew in the Strait of Gibraltar throuhg a physical-biological coupled model. Prog. Oceanogr., Vol. 74, (2007), pp. 252-272. [6] E. Ram´ırez-Romero , M. Vichi, M.J. Castro, J. Mac´ıas, D. Mac´ıas, C.M. Garc´ıa y M. Bruno. Modeling biogeochemical seasonal cycle in the Strait of Gibraltar. Enviado (2013). [7] J.B. Schijf y J.C. Schonfeld. Theoretical considerations on the motion of salt and fresh water. En Proc. of the Minn. Int. Hydraulics Conv., 321-333. Joint meeting IAHR and Hyd. Div. ASCE., Sept. 1953. 8