Approximation of the transport of pollutants with reaction terms in shallow waters using finite elements

<!-- metadata commented in wiki content
==Aproximación del transporte de contaminantes con términos de reacción en aguas someras mediante elementos finitos==

'''Angel Villota<math>^{a,b}</math>, Ramon Codina<math>^{a}</math>'''
<math>^{a}</math>Universidad Politécnica de Cataluña, Jordi Girona 1-3, Edificio C1, 08034 Barcelona, Spain 

<math>^{b}</math>Escuela Politécnica Nacional, Ladrón de Guevara E 11-253, Quito, Ecuador
-->
==Resumen==

En este artículo presentamos la aproximación del modelo acoplado de las ecuaciones del movimiento de un fluido en aguas poco profundas con la ecuación convección-difusión-reacción (CDR) del transporte de contaminantes. Dicha aproximación se realiza mediante elementos finitos de alto orden y usando métodos variacionales estabilizados de subescalas. El sistema acoplado de ecuaciones, previamente discretizado en el tiempo y linealizado, lo escribimos como una ecuación vectorial transitoria de CDR. Los métodos estabilizados de elementos finitos utilizados son los conocidos métodos de subescalas ASGS y OSS, los mismos que nos permiten usar igual interpolación para todas las incógnitas, así como tratar con flujos de convección y reacción dominantes. Consideramos la posibilidad de no linealidad tanto en el término convectivo como en el de reacción. No consideraremos el posible desarrollo de choques en la solución. Con el fin de examinar la precisión y robustez de los métodos ASGS y OSS, presentamos cuatro casos de prueba: convergencia en malla, transporte de un contaminante en una cavidad cuadrada, transporte de un contaminante en el golfo de Roses y en la desembocadura del río Guadalquivir, y el modelo depredador-presa, que puede escribirse como una ecuación vectorial de CDR transitoria con no linealidad en el término de reacción.

'''Palabras clave''': Aguas poco profundas, ablandamiento y estabilización, transporte de contaminantes, modelo depredador-presa 

==Abstract==

In this article we present the approximation of the coupled model of the equations of motion of a fluid in shallow waters with the convection-diffusion-reaction (CDR) equation of pollutant transport. This approximation is carried out using high order finite elements and using stabilised variational sub-scale methods. We write the coupled system of equations, previously discretised in time and linearised, as a transient vector equation of CDR. The stabilised finite element methods used are the known ASGS and OSS sub-scale methods, the same ones that allow us to use the same interpolation for all unknowns, as well as to deal with dominant convection and reaction flows. We consider the possibility of non-linearity in both the convective and reaction terms. We will not consider the possible development of shocks in the solution. In order to examine the accuracy and robustness of the ASGS and OSS methods, we present four test cases: mesh convergence, transport of a pollutant in a square cavity, transport of a pollutant in the Gulf of Roses and at the river Guadalquivir mouth, and the predator-prey model, which can be written as a transient CDR vector equation with non-linearity in the reaction term.

'''Keywords''': Shallow waters, softening and stabilization, transport of pollutants, predator-prey model

==1. Introducción==

Un contaminante es aquel componente que está presente en el agua a niveles perjudiciales para la vida de los seres humanos, plantas y animales (<span id='citeF-1'></span>[[#cite-1|[1]]], véase también los informes de la  U.S. Environmental Protecction Agency, https://www.epa.gov/). La simulación del transporte de contaminantes para la predicción de la concentración de contaminantes en ríos, lagos, lagunas y regiones costeras tiene importancia estratégica en el análisis y diseño de soluciones de los problemas de contaminación ambiental del planeta.

En el fenómeno físico del movimiento de contaminantes se deben distinguir tres procesos: la difusión, la convección y la reacción.  La difusión es el proceso físico debido al cual el contaminante se mueve como resultado del movimiento intermolecular de las partículas de ambas sustancias, el fluido que la transporta (agua) y el contaminante. La convección es el movimiento del soluto (contaminante) debido al movimiento del agua, por lo cual si el agua permanece en reposo no hay convección. Finalmente, la reacción tiene en cuenta el posible efecto de crecimiento o decrecimiento de un contaminante por factores externos (calor, especies con las que puede combinarse); en el caso de haber distintas especies de contaminantes, la reacción puede modelar la interacción entre ellas.

La formulación del problema del transporte de contaminantes se fundamenta, como la de todos los fenómenos físicos, en las ecuaciones de equilibrio y las ecuaciones constitutivas. Las ecuaciones del transporte de sustancias disueltas o en suspensión en el flujo se basan en el principio de conservación de la masa de dichas sustancias en caso de que no haya reacción  o en el modelo de crecimiento o decrecimiento en caso de que la haya.  En el problema que queremos abordar, la dinámica del transporte de contaminantes involucra los modelos de flujos del campo de velocidades en aguas someras y del transporte, difusión y reacción de contaminantes, los cuales tienen características de la ecuación transitoria de convección, difusión y reacción (CDR).

En este trabajo vamos a utilizar un modelo matemático que acopla las ecuaciones del movimiento de un fluido en aguas poco profundas con la ecuación de CDR del transporte de contaminantes. Las ecuaciones del movimiento de un fluido en aguas poco profundas que vamos a utilizar son las obtenidas a partir de las ecuaciones de Navier-Stokes para fluidos viscosos, y que se encuentran descritas por ejemplo en <span id='citeF-2'></span>[[#cite-2|[2]]], acopladas a la ecuación de CDR del transporte de contaminantes como en <span id='citeF-3'></span><span id='citeF-4'></span><span id='citeF-5'></span>[[#cite-3|[3,4,5]]], entre otras referencias. El sistema acoplado de ecuaciones se puede escribir en forma conservativa en términos de la profundidad del fluido, de las velocidades y de las concentraciones de contaminantes promediadas en el plano donde el flujo medio tiene lugar. La presencia del término fuente en la ecuación de momento, la consideración de la batimetría variable, la fricción en el fondo del lecho, la tensión en la superficie libre del agua debida al viento y el efecto de Coriolis son características que se pueden tomar en cuenta para obtener modelos físicos más realistas del flujo de aguas someras. El resultado es un sistema no lineal acoplado de ecuaciones diferenciales parciales. Para simplificar el modelo, supondremos que los contaminantes son pasivos, en el sentido de que no afectan la densidad del fluido y por lo tanto el flujo; en caso de no ser así, se podría introducir fácilmente algún modelo de flotación que tuviera en cuenta el flujo generado por la variabilidad en la concentración de contaminantes.

En la literatura existe mucha investigación numérica en el área del transporte de contaminantes en aguas poco profundas, aplicando métodos de diferencias finitas <span id='citeF-6'></span>[[#cite-6|[6]]], métodos de volúmenes finitos <span id='citeF-7'></span><span id='citeF-8'></span><span id='citeF-11'></span><span id='citeF-12'></span><span id='citeF-13'></span>[[#cite-7|[7,8,11,12,13]]], métodos de elementos finitos <span id='citeF-4'></span><span id='citeF-9'></span><span id='citeF-10'></span>[[#cite-4|[4,9,10]]], métodos tipo-Godunov <span id='citeF-18'></span>[[#cite-18|[18]]] y métodos de Boltzmann <span id='citeF-14'></span><span id='citeF-15'></span><span id='citeF-16'></span>[[#cite-14|[14,15,16]]].

En este trabajo vamos a realizar la aproximación de las ecuaciones del movimiento de un fluido en aguas poco profundas acopladas a la ecuación de CDR del transporte de contaminantes mediante los métodos VMS (''Variational Muli-Scale'') llamados ASGS (''Algebraic SubGrid Scale'') y OSS (''Orthogonal Subgrid Scale'') (véase <span id='citeF-17'></span>[[#cite-17|[17]]]). El propósito de este trabajo es tanto aplicar estas técnicas, ya conocidas, a problemas no lineales como el transporte de contaminantes en aguas poco profundas y el modelo depredador-presa, como experimentar con elementos finitos de orden elevado (hasta el cuarto orden).

Veremos en primer lugar que mediante un cambio de variables apropiado y linealizando el sistema acoplado de ecuaciones en aguas poco profundas con las ecuaciones del transporte de contaminantes, cada iteración del proceso no lineal puede escribirse como una ecuación vectorial de CDR transitoria lineal. Plantearemos la aproximación en el tiempo y el espacio de esta ecuación vectorial de CDR transitoria, para luego particularizarla al sistema acoplado de las ecuaciones del flujo en aguas poco profundas con el transporte de contaminantes.

En lo concerniente a la discretización, primero procederemos a discretizar en el tiempo mediante un esquema de integración temporal en diferencias finitas y luego realizaremos la aproximación en el espacio mediante elementos finitos, aunque el enfoque más común es proceder a la inversa, es decir, primero discretizar en el espacio y luego aproximar el sistema resultante de ecuaciones diferenciales en el tiempo. Sin embargo, el enfoque de discretizar primero en el tiempo nos permite desacoplar los errores provenientes de la discretización temporal de los correspondientes a la discretización espacial. Aunque ambos caminos para el método estandar de Galerkin dan el mismo problema discreto, esto no es así para todos los métodos basados en VMS. En lo referente a la discretización temporal, usaremos la regla trapezoidal generalizada, que es el método en diferencias finitas más sencillo; sin embargo, cualquier otro método en diferencias finitas o incluso en elementos finitos podría aplicarse. En particular, en los ejemplos numéricos usamos esquemas de integración temporal de mayor orden que la regla trapezoidal.

En cuanto a la discretización espacial, es bien conocido que el método estándar de Galerkin presenta inestabilidades numéricas cuando los términos convectivo o de reacción son dominantes con respecto al término difusivo, además de por la necesidad dada por la condición inf-sup de usar diferentes interpolaciones de elementos finitos para la profundidad, las velocidades promediadas y la concentración de contaminante. Estos dos inconvenientes pueden superarse utilizando una formulación estabilizada. Aquí vamos a usar las formulaciones basadas en subescalas introducidas por Hughes y otros en <span id='citeF-20'></span><span id='citeF-21'></span>[[#cite-20|[20,21]]] para la ecuación CDR escalar estacionaria y generalizado a sistemas lineales estacionarios en <span id='citeF-22'></span><span id='citeF-23'></span>[[#cite-22|[22,23]]], y a sistemas lineales transitorios en <span id='citeF-24'></span><span id='citeF-25'></span><span id='citeF-26'></span>[[#cite-24|[24,25,26]]], dando lugar a las formulaciones ASGS y OSS para resolver el problema de CDR vectorial transitorio, como en <span id='citeF-27'></span>[[#cite-27|[27]]]. La idea básica es aproximar el efecto de la componente de la solución continua que no puede ser resuelta por la malla de elementos finitos en términos de la solución en el espacio de elementos finitos (método ASGS), o en términos de una subescala ortogonal respecto del producto escalar de <math display="inline">L^{2}</math> al espacio de elementos finitos  <span id='citeF-28'></span><span id='citeF-29'></span>[[#cite-28|[28,29]]] (método OSS). Hay que hacer notar que no consideraremos el posible desarrollo de choques en la solución, es decir, discontinuidades que pueden llegar a producirse en el límite invíscido. Para tomar en cuenta esto es necesario desarrollar técnicas de captura de discontinuidades que no consideraremos en este trabajo, así como tampoco el tratamiento de subescalas dinámicas.

El artículo está organizado como sigue. En la siguiente sección se describe la física del problema a ser resuelto y la obtención de la ecuación vectorial de CDR transitoria que representa el problema acoplado de las ecuaciones del movimiento de un fluido en aguas poco profundas con las ecuaciones del transporte de contaminantes. Las aproximaciones y formulaciones de los métodos variacionales estabilizados ASGS y OSS para la ecuación vectorial de CDR transitoria están descritos en la sección 3. En la sección 4 se presentan pruebas numéricas de convergencia en malla con soluciones analíticas conocidas, un ejemplo clásico de la literatura del transporte de un contaminante en una cavidad cuadrada, dos ejemplos prácticos del transporte de un contaminante, uno en el golfo de Roses y otro en la desembocadura del río Guadalquivir, y un ejemplo del modelo depredador-presa. Finalmente, en base a los resultados obtenidos, en la Sección 5 se detallan algunas conclusiones.

==2. Planteamiento del problema del transporte de contaminantes==

===2.1 Ecuaciones de gobierno===

Las ecuaciones del movimiento de un fluido en aguas poco profundas se obtienen por integración vertical de las ecuaciones del movimiento, bien sea de fluidos viscosos (ecuaciones de Navier-Stokes) o invíscidos (ecuaciones de Euler); véase por ejemplo <span id='citeF-2'></span><span id='citeF-19'></span>[[#cite-2|[2,19]]]. En este trabajo vamos a utilizar las ecuaciones obtenidas a partir de las ecuaciones de Navier-Stokes para un fluido Newtoniano incompresible e isotérmico, incluyendo la tensión superficial provocada por los efectos del viento así como también los efectos de Coriolis, al igual que batimetría variable, fricción en el fondo del lecho y distribución hidrostática de presiones. Vamos a considerar el problema de mecánica de fluidos acoplado con el de transporte  de contaminantes. Esta será en esencia la diferencia del presente trabajo con respecto a <span id='citeF-27'></span>[[#cite-27|[27]]], por lo que muchos de los desarrollos que presentamos aquí se encuentran también en esta referencia. Omitiremos aquellos desarrollos que no afectan la compresión de este artículo, aunque para permitir la lectura auto-contenida se repetirán muchos otros.

Sea <math display="inline">\Omega </math> un dominio acotado <math display="inline">\mathbb{R}^{2}</math> con coordenadas cartesianas <math display="inline">\mathbf{x}=(x_{1},x_{2}).</math> Consideremos un fluido que se mueve durante el intervalo de tiempo <math display="inline">[0,T]</math> en el dominio de <math display="inline">\mathbb{R}^{3}</math> dado por <math display="inline">\Omega \times (-H(\mathbf{x}),\eta (\mathbf{x},t)),</math> donde la coordenada <math display="inline">x_{3}=-H(\mathbf{x})</math> representa la batimetría del terreno y <math display="inline">x_{3}=\eta (\mathbf{x},t)</math>, que será incógnita del problema, es la superficie libre del fluido con <math display="inline">\mathbf{x}\in \Omega </math>, <math display="inline">t\in [0,T]</math>. Consideraremos que <math display="inline">h(\mathbf{x},t):=H(\mathbf{x})+\eta (\mathbf{x},t)</math> es mucho menor que el diámetro de <math display="inline">\Omega </math>, por lo que supondremos válida la aproximación de las ecuaciones de movimiento en aguas poco profundas. Dichas ecuaciones acopladas a la ecuación del transporte de contaminantes pueden escribirse en forma indicial como

<span id="eq-1"></span>
<span id="eq-2"></span>
<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _{t}\left(hU_{i}\right)+\partial _{j}\left(hU_{i}U_{j}\right)+\partial _{i}\frac{1}{2}g\left(h^{2}-H^{2}\right)-\partial _{j}\left(\frac{h}{\rho }\bar{\tau }_{ji}\right)+Q_{i}=0,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|-
| style="text-align: center;" | <math>  \partial _{t}h+\partial _{i}\left(hU_{i}\right)=0,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|-
| style="text-align: center;" | <math>  \partial _{t}\left(h\varphi _{\alpha }\right)+\partial _{i}\left(hU_{i}\varphi _{\alpha }\right)-\partial _{i}\left(hk_{ij}^{\alpha }\partial _{j}\varphi _{\alpha }\right)+S_{\alpha \beta }\varphi _{\beta }=0, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}
|}

con

<span id="eq-4"></span>
<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\overline{\tau }_{ij}=\mu _{H}\left(\partial _{j}U_{i}+\partial _{i}U_{j}-\frac{2}{3}\delta _{ij}\partial _{k}U_{k}\right),</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|-
| style="text-align: center;" | <math>  Q_{i}=h\hat{f}\left(-\delta _{i1}U_{2}+\delta _{i2}U_{1}\right)-g\left(h-H\right)\partial _{i}H+\frac{h}{\rho }\partial _{i}p_{a}-\frac{1}{\rho }\tau _{3i}^{s}+gU_{i}U_{0}\mathrm{Ch^{-1}}h^{-2},\quad i=1,2, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}
|}

y

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>U_{0}=\sqrt{U_{1}^{2}+U_{2}^{2}}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

donde <math display="inline">U_{i}(\mathbf{x},t)</math> <math display="inline">\left(i=1,2\right)</math> son las componentes de la velocidad promediadas en la profundidad, <math display="inline">\delta _{ij}</math> es la delta de Kronecker, <math display="inline">\partial _{t}</math> denota la derivada temporal y <math display="inline">\partial _{i}</math> la derivada con respecto a la <math display="inline">i</math>-ésima coordenada cartesiana <math display="inline">(i=1,2)</math> en <math display="inline">\Omega </math>. En las ecuaciones anteriores y en lo que sigue, índices repetidos indican la notación de sumación de Einstein, y los subíndices <math display="inline">i,</math> <math display="inline">j,</math> <math display="inline">k</math> recorren desde 1 hasta 2. En las expresiones previas, <math display="inline">g</math> es la aceleración de la gravedad, <math display="inline">\rho </math> la densidad del fluido, <math display="inline">\bar{\tau }_{ij}</math> son las tensiones viscosas promediadas, <math display="inline">\mu _{H}</math> es el coeficiente de viscosidad dinámico, <math display="inline">\hat{f}</math> es el parámetro de Coriolis, <math display="inline">p_{a}</math> es la presión atmosférica, <math display="inline">\tau _{3i}^{s}</math> <math display="inline">\left(i=1,2\right)</math> son las tensiones en la superficie libre del agua debidas al viento, Ch es el coeficiente de Chézy para el rozamiento en el fondo y <math display="inline">U_{0}</math> es la norma euclídea de <math display="inline">\left(U_{1},U_{2}\right)</math>. Usaremos letras griegas como subíndices para caracterizar los contaminantes, cuyo número supondremos que es <math display="inline">N</math>. Así,  la concentración del contaminante <math display="inline">\alpha </math> promediada en la profundidad se representa por <math display="inline">\varphi _{\alpha }</math>. <math display="inline">S_{\alpha \beta }</math> son las componentes de la matriz de términos de reacción de la ecuación de conservación  de contaminantes, con <math display="inline">\alpha </math>, <math display="inline">\beta = 1,\dots ,N</math>; en los casos de reacción no lineal, <math display="inline">S_{\alpha \beta }</math> puede depender de las concentraciones de contaminantes, es decir, <math display="inline">S_{\alpha \beta } = S_{\alpha \beta }(\varphi _1,\dots ,\varphi _N)</math>. <math display="inline">k_{ij}^{\alpha }</math> es el coeficiente empírico de difusión (dispersión) para el contaminante <math display="inline">\alpha </math>; nótese que por simplificar hemos supuesto que la dispersión del contaminante <math display="inline">\alpha </math> no se ve afectada por el resto de contaminantes, aunque esta situación podríamos incluirla fácilmente.

Las ecuaciones [[#eq-1|(1)]], [[#eq-2|(2)]] y [[#eq-3|(3)]] representan un sistema de <math display="inline">(3+N)</math> ecuaciones, cuyo vector de incógnitas es

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{U}^{T}=\left[U_{1},U_{2},h,\varphi _1,\dots ,\varphi _N\right]. </math>
|}
|}

Dichas ecuaciones constituyen el modelo acoplado del movimiento de un fluido en aguas poco profundas y el transporte de contaminantes. Las mismas se deben completar con las condiciones iniciales y de contorno. Los valores iniciales para la elevación del agua, las velocidades promediadas y la concentración del contaminante deben ser conocidos como condiciones iniciales. Con respecto a las condiciones de contorno, es bien conocido que el problema está bien puesto si la velocidad promediada se especifica en el contorno de entrada del dominio computacional <math display="inline">\Omega </math>, mientras que en el resto del contorno se prescribe la componente normal de las tensiones viscosas y la elevación del fluido. Para las concentraciones de contaminantes, además de la condición inicial es necesario conocer como condiciones de contorno sus valores en el contorno de entrada del dominio de cálculo y los flujos en el contorno de salida; esta situación permitiría en particular considerar nulas las difusiones de contaminantes.

Debido a la complejidad del modelo descrito anteriormente, en la literatura existen muchas variantes de estas ecuaciones con distintas simplificaciones, como por ejemplo eliminar el término de difusión considerando convección pura; este es el caso de <span id='citeF-3'></span><span id='citeF-4'></span><span id='citeF-5'></span>[[#cite-3|[3,4,5]]]. En nuestro trabajo utilizaremos las ecuaciones ([[#eq-1|1]]), ([[#eq-2|2]]) y ([[#eq-3|3]]) sin eliminar o hacer simplificaciones de ninguno de sus términos. Mediante un cambio de variables transformaremos las ecuaciones para poder escribirlas como una ecuación vectorial de CDR transitoria de la forma

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _{t}\mathbf{u}-\partial _{i}\left(\mathbf{K}_{ij}\partial _{j}\mathbf{u}\right)+\partial _{i}\left(\mathbf{A}_{i}\mathbf{u}\right)+\mathbf{Su}=\mathbf{F}\quad i,j=1,2, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

para la que más adelante, en la Sección [[#3 Aproximación espacial y estabilización numérica|3]], presentaremos su aproximación numérica por elementos finitos mediante los métodos variacionales de subescalas ASGS y OSS. En las ecuaciones anteriores, <math display="inline">\mathbf{u}</math> es el vector de incógnitas, <math display="inline">\mathbf{F}</math> el vector de fuerzas conocido y <math display="inline">\mathbf{K}_{ij}</math>, <math display="inline">\mathbf{A}_{i}</math>, <math display="inline">\mathbf{S}</math> matrices dadas, que pueden depender de <math display="inline">\mathbf{u}</math> en problemas no lineales.

===2.2 Cambio de variables===

Realizando un cambio de variables, <math display="inline">\varphi _{a,\alpha }:=h\varphi _{\alpha }</math>, <math display="inline">u_{i}:=hU_{i}</math> y <math display="inline">P:=\frac{1}{2}g\left(h^{2}-H^{2}\right)</math>, con <math display="inline">\partial _{t}P=hg\partial _{t}h,</math> las ecuaciones [[#eq-1|(1)]], [[#eq-2|(2)]] y ([[#eq-3|3]]) se pueden escribir de la siguiente manera:

<span id="eq-8"></span>
<span id="eq-9"></span>
<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _{t}u_{i}+\partial _{j}\left(U_{j}u_{i}\right)+\partial _{i}P-\partial _{j}\frac{\mu _{H}}{\rho }\left[\left(\partial _{i}u_{j}+\partial _{j}u_{i}-\frac{2}{3}\delta _{ij}\partial _{k}u_{k}\right)-\left(U_{j}\partial _{i}h+U_{i}\partial _{j}h-\frac{2}{3}\delta _{ij}U_{k}\partial _{k}h\right)\right] </math>
|-
| style="text-align: center;" | <math>  \qquad +\hat{f}\left(-\delta _{i1}u_{2}+\delta _{i2}u_{1}\right)+gU_{0}\mathrm{Ch}^{-1}h^{-3}u_{i}-g\left(h-H\right)\partial _{i}H+\frac{h}{\rho }\partial _{i}p_{a}-\frac{1}{\rho }\tau _{3i}^{s}=0,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|-
| style="text-align: center;" | <math>  \frac{1}{hg}\partial _{t}P+\partial _{i}u_{i}=0,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|-
| style="text-align: center;" | <math>  \partial _{t}\varphi _{a,\alpha }+\partial _{i}\left(U_{i}\varphi _{a,\alpha }\right)-\partial _{i}\left(k_{ij}^{\alpha }\partial _{j}\varphi _{a,\alpha }-k_{ij}^{\alpha }\varphi _{\alpha }\partial _{j}h\right)+S_{\alpha \beta }h^{-1}\varphi _{a,\beta }=0. </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}
|}

Las ecuaciones ([[#eq-8|8]]), ([[#eq-9|9]]) y ([[#eq-10|10]]) se pueden escribir fácilmente en el formato de la ecuación ([[#eq-7|7]]). Sin embargo esto lo haremos en el caso del problema linealizado y discretizado en el tiempo.

===2.3 Discretización en el tiempo y linealización===

Para la discretización temporal del modelo acoplado de las ecuaciones del movimiento de un fluido en aguas poco profundas y del transporte de contaminantes, consideremos una partición uniforme del intervalo de tiempo <math display="inline">0=t^{0}<t^{1}<\dots{<}t^{M}=T</math>, con <math display="inline">\delta t=t^{n+1}-t^{n}</math>, donde <math display="inline">\delta t</math> es el tamaño del paso de tiempo, considerado constante. Aquí y en adelante, usaremos el superíndice <math display="inline">n</math> para indicar el nivel de tiempo&nbsp;<math display="inline">n</math>.

Aunque la discretización en el tiempo puede realizarse mediante cualquier aproximación de la derivada temporal en diferencias finitas o elementos finitos, aquí utilizaremos la regla trapezoidal generalizada. Para una función genérica <math display="inline">f</math>, sea

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f^{n+\theta }=\theta f^{n+1}+\left(1-\theta \right)f^{n},\qquad{0}\leq \theta \leq{1.} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

Aproximaremos

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _{t}^{n}f:=\frac{1}{\delta t}(f^{n+1}-f^{n})=\frac{1}{\theta \delta t}(f^{n+\theta }-f^{n})\approx (\partial _{t}f)^{n+\theta } </math>
|}
|}

Esta aproximación es de segundo orden en <math display="inline">\delta t</math> en el caso en que <math display="inline">\theta=1/2</math> (método de Crank-Nicolson) y de primer orden si <math display="inline">\theta \not=1/2</math>. Para <math display="inline">\theta=1</math> corresponde al método de Euler implícito o el método de diferencias hacia atrás de primer orden (BDF1).

En cuanto a la linealización, utilizaremos el método de punto fijo de Picard. En cada iteración de cada paso de tiempo, el problema a resolver es de la forma:

<span id="eq-12"></span>
<span id="eq-13"></span>
<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _{t}^{n}u_{i}+\partial _{j}\left(a_{j}u_{i}\right)+\partial _{i}p-\partial _{j}\hat{\tau }_{ji}+S_{u,ij}u_{j}+b_{i}  =0,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|-
| style="text-align: center;" | <math>  \frac{1}{h_{0}g}\delta _{t}^{n}p+\partial _{i}u_{i}  =0,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|-
| style="text-align: center;" | <math>  \delta _{t}^{n}\varphi _{a,\alpha }+\partial _{j}\left(a_{j}\varphi _{a,\alpha }\right)-\partial _{i}\left(k_{ij}^{\alpha }\partial _{j}\varphi _{a,\alpha }-k_{ij}^{\alpha }\varphi _{0}^{\alpha }\partial _{j}h_{0}\right)+S_{\alpha \beta }h^{-1}\varphi _{a,\beta }=0. </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}
|}

Usando un superíndice adicional para el contador de iteraciones, las variables introducidas en ([[#eq-12|12]]), ([[#eq-13|13]]) y ([[#eq-14|14]]) para la iteración <math display="inline">i</math>-ésima en el paso de tiempo <math display="inline">n</math> son

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{2}   a_{j}=U_{j}^{n+\theta ,i-1},\quad    u_{j}=\left(hU_{j}\right)^{n+\theta ,i},</math>
|-
| style="text-align: center;" | <math>   h_{0}=h^{n+\theta ,i-1},\quad    p=P^{n+\theta ,i},</math>
|-
| style="text-align: center;" | <math>   \delta _{t}^{n}p=\frac{1}{\theta \delta t}\left(p-P^{n}\right),\quad    \delta _{t}^{n}u_{j}=\frac{1}{\theta \delta t}\left[u_{j}-\left(hU_{j}\right)^{n}\right],</math>
|-
| style="text-align: center;" | <math>   \hat{\tau }_{ij}=\frac{\mu _{H}}{\rho }\left(\partial _{i}u_{j}+\partial _{j}u_{i}-\frac{2}{3}\delta _{ij}\partial _{k}u_{k}\right),\quad    \tau _{ij}^{*}=\frac{\mu _{H}}{\rho }\left(a_{j}\partial _{i}h_{0}+a_{i}\partial _{j}h_{0}-\frac{2}{3}\delta _{ij}a_{k}\partial _{k}h_{0}\right),</math>
|-
| style="text-align: center;" | <math>   S_{u,11}=S_{u,22}=\mathrm{Ch}^{-1}gU_{0}^{n+\theta ,i-1}h_{0}^{-3},\quad    U_{0}^{n+\theta ,i-1}=\sqrt{\left(U_{1}^{n+\theta ,i-1}\right)^{2}+\left(U_{2}^{n+\theta ,i-1}\right)^{2}},</math>
|-
| style="text-align: center;" | <math>   S_{u,12}=-S_{u,21}=-\hat{f},\quad    b_{i}=b_{i}^{*}+\partial _{j}\tau _{ji}^{*},</math>
|-
| style="text-align: center;" | <math>   b_{i}^{*}=-g\left(h_{0}-H\right)\partial _{i}H+\frac{1}{\rho }h_{0}\partial _{i}p_{a}-\frac{1}{\rho }\tau _{3i}^{S},\quad    \varphi _{0}^{\alpha }=\varphi _{\alpha }^{n+\theta ,i-1}, </math>
|-
| style="text-align: center;" | <math>  \delta _{t}^{n}\varphi _{a,\alpha }=\frac{1}{\theta \delta{t}}\left[\varphi _{a,\alpha }-\left(h\varphi _{\alpha }\right)^{n}\right],   \varphi _{a,\alpha }=\left(h\varphi _{\alpha }\right)^{n+\theta ,i}. </math>
|}
|}

El super-índice <math display="inline">i</math> es el contador de las iteraciones debidas a la no linealidad. Obsérvese que la forma linealizada de las tensiones viscosas <math display="inline">\overline{\tau }_{ij}</math> dada en la ecuación ([[#eq-4|4]]) corresponde a <math display="inline">\hat{\tau }_{ij}-\tau _{ij}^{*},</math> donde <math display="inline">\hat{\tau }_{ij}</math> está escrito en términos de las variables de la iteración actual, mientras que <math display="inline">\tau _{ij}^{*}</math> se calcula con los valores de las variables correspondientes a la iteración anterior. Nótese que hemos usado el símbolo <math display="inline">S_u</math> para indicar las componentes de la matriz de términos reactivos provenientes de la ecuación de conservación de la cantidad de movimiento, para distinguirlas de las provenientes de la ecuación de conservación de contaminantes, aunque ambas contribuirán a la matriz de términos reactivos en la forma vectorial que presentamos a continuación.

===2.4 Forma vectorial===

Para escribir las ecuaciones ([[#eq-12|12]]), ([[#eq-13|13]]) y ([[#eq-14|14]]) en la forma vectorial ([[#eq-7|7]]), aunque discretizada en el tiempo, escribimos las mismas como un sistema de <math display="inline">(3+N)</math> ecuaciones con <math display="inline">(3+N)</math> incógnitas. Las ecuaciones correspondientes a las velocidades y altura de la lámina de agua son las mismas que en <span id='citeF-27'></span>[[#cite-27|[27]]], mientras que la ecuación correspondiente a la concentración de los contaminantes es:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _{t}^{n}\varphi _{a,\alpha }-k_{11}^{\alpha }\partial _{1}^{2}\varphi _{a,\alpha }-k_{12}^{\alpha }\partial _{1}\partial _{2}\varphi _{a,\alpha }-k_{21}^{\alpha }\partial _{2}\partial _{1}\varphi _{a,\alpha }-k_{22}^{\alpha }\partial ^{2}\varphi _{a,\alpha }+a_{1}\partial _{1}\varphi _{a,\alpha }+a_{2}\partial _{2}\varphi _{a,\alpha }+\left(\partial _{1}a_{1}+\partial _{2}a_{2}\right)\varphi _{a,\alpha } </math>
|-
| style="text-align: center;" | <math>  \qquad +k_{11}^{\alpha }\left(\partial _{1}\varphi _{0}^{\alpha }\partial _{1}h_{0}+\varphi _{0}^{\alpha }\partial _{1}^{2}h_{0}\right)+k_{12}^{\alpha }\left(\partial _{1}\varphi _{0}^{\alpha }\partial _{2}h_{0}+\varphi _{0}^{\alpha }\partial _{1}\partial _{2}h_{0}\right)</math>
|-
| style="text-align: center;" | <math>  \qquad +k_{21}^{\alpha }\left(\partial _{2}\varphi _{0}^{\alpha }\partial _{1}h_{0}+\varphi _{0}^{\alpha }\partial _{2}\partial _{1}h_{0}\right)+k_{22}^{\alpha }\left(\partial _{2}\varphi _{0}^{\alpha }\partial _{2}h_{0}+\varphi _{0}^{\alpha }\partial _{2}^{2}h_{0}\right)+S_{\alpha \beta }h_0^{-1}\varphi _{a,\beta }=0. </math>
|}
|}

En esta ecuación hay suma para el índice <math display="inline">\beta </math>, pero no para <math display="inline">\alpha </math>, que es el contador de contaminante considerado.

Para simplificar la notación, consideremos un único contaminante, caracterizado por el índice genérico <math display="inline">\alpha </math>. La forma vectorial ([[#eq-7|7]]) del sistema de ecuaciones final, previamente discretizadas en el tiempo y linealizadas con el método de punto fijo de Picard, está dada por:

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[\begin{array}{c}\delta _{t}^{n}u_{1}\\ \delta _{t}^{n}u_{2}\\ \begin{array}{c}\frac{1}{h_{0}g}\delta _{t}^{n}p\\ \delta _{t}^{n}\varphi _{a,\alpha } \end{array} \end{array}\right] -\frac{\mu _{H}}{\rho }\left[\begin{array}{cccc}\frac{4}{3} & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & k_{11}^{*} \end{array}\right]\partial _{1}^{2}\left[\begin{array}{c}u_{1}\\ u_{2}\\ \begin{array}{c}p\\ \varphi _{a,\alpha } \end{array} \end{array}\right]-\frac{\mu _{H}}{\rho }\left[\begin{array}{cccc}0 & \frac{1}{6} & 0 & 0\\ \frac{1}{6} & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & k_{12}^{*} \end{array}\right]\partial _{1}\partial _{2}\left[\begin{array}{c}u_{1}\\ u_{2}\\ \begin{array}{c}p\\ \varphi _{a,\alpha } \end{array} \end{array}\right]</math>
|-
| style="text-align: center;" | <math>  \qquad  -\frac{\mu _{H}}{\rho }\left[\begin{array}{cccc}1 & 0 & 0 & 0\\ 0 & \frac{4}{3} & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & k_{22}^{*} \end{array}\right]\partial _{2}^{2}\left[\begin{array}{c}u_{1}\\ u_{2}\\ \begin{array}{c}p\\ \varphi _{a,\alpha } \end{array} \end{array}\right]-\frac{\mu _{H}}{\rho }\left[\begin{array}{cccc}0 & \frac{1}{6} & 0 & 0\\ \frac{1}{6} & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & k_{21}^{*} \end{array}\right]\partial _{2}\partial _{1}\left[\begin{array}{c}u_{1}\\ u_{2}\\ \begin{array}{c}p\\ \varphi _{a,\alpha } \end{array} \end{array}\right]</math>
|-
| style="text-align: center;" | <math>  \qquad +\left[\begin{array}{cccc}a_{1} & 0 & 1 & 0\\ 0 & a_{1} & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 0 & 0 & a_{1} \end{array}\right]\partial _{1}\left[\begin{array}{c}u_{1}\\ u_{2}\\ \begin{array}{c}p\\ \varphi _{a,\alpha } \end{array} \end{array}\,\right]+\left[\begin{array}{cccc}a_{2} & 0 & 0 & 0\\ 0 & a_{2} & 1 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 0 & a_{2} \end{array}\right]\partial _{2}\left[\begin{array}{c}u_{1}\\ u_{2}\\ \begin{array}{c}p\\ \varphi _{a,\alpha } \end{array} \end{array}\right]</math>
|-
| style="text-align: center;" | <math>  \qquad  +\left[\begin{array}{cccc}\partial _{1}a_{1}+\partial _{2}a_{2}+S_{u,11} & S_{u,12} & 0 & 0\\ S_{u,21} & \partial _{1}a_{1}+\partial _{2}a_{2}+S_{u,22} & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \partial _{1}a_{1}+\partial _{2}a_{2} + S_{\alpha \alpha } h_0^{-1} \end{array}\right]\left[\begin{array}{c}u_{1}\\ u_{2}\\ \begin{array}{c}p\\ \varphi _{a,\alpha } \end{array} \end{array}\right]=\left[\begin{array}{c}-Q_{1}\\ -Q_{2}\\ \begin{array}{c}0\\ -Q_{4} \end{array} \end{array}\right], </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

con

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k_{11}^{*}=\frac{\rho }{\mu _{H}}k_{11}^{\alpha },\quad k_{12}^{*}=\frac{\rho }{\mu _{H}}k_{12}^{\alpha },\quad k_{21}^{*}=\frac{\rho }{\mu _{H}}k_{21}^{\alpha },\quad k_{22}^{*}=\frac{\rho }{\mu _{H}}k_{22}^{\alpha },</math>
|-
| style="text-align: center;" | <math>  Q_{1}=b_{1}^{*}+\partial _{1}\tau _{11}^{*}+\partial _{2}\tau _{21}^{*}, </math>
|-
| style="text-align: center;" | <math>  Q_{2}=b_{2}^{*}+\partial _{1}\tau _{12}^{*}+\partial _{2}\tau _{22}^{*}, </math>
|-
| style="text-align: center;" | <math>  Q_{4}=k_{11}^{\alpha }\left(\partial _{1}\varphi _{0}^{\alpha }\partial _{1}h_{0}+\varphi _{0}^{\alpha }\partial _{1}^{2}h_{0}\right)+k_{12}^{\alpha }\left(\partial _{1}\varphi _{0}^{\alpha }\partial _{2}h_{0}+\varphi _{0}^{\alpha }\partial _{1}\partial _{2}h_{0}\right)</math>
|-
| style="text-align: center;" | <math>  \qquad +k_{21}^{\alpha }\left(\partial _{2}\varphi _{0}^{\alpha }\partial _{1}h_{0}+\varphi _{0}^{\alpha }\partial _{2}\partial _{1}h_{0}\right)+k_{22}^{\alpha }\left(\partial _{2}\varphi _{0}^{\alpha }\partial _{2}h_{0}+\varphi _{0}^{\alpha }\partial _{2}^{2}h_{0}\right). </math>
|}
|}

En caso de considerar <math display="inline">N</math> contaminantes, el coeficiente 4-4 de la matriz de términos reactivos sería de hecho una matriz <math display="inline">N \times N</math>, con el coeficiente <math display="inline">\alpha \beta </math> dado por <math display="inline">(\partial _{1}a_{1}+\partial _{2}a_{2}) \delta _{\alpha \beta } + S_{\alpha \beta } h_0^{-1}</math>.

La ecuación ([[#eq-15|15]]) corresponde a la ecuación vectorial de CDR transitoria ([[#eq-7|7]]) de las ecuaciones acopladas del movimiento de un fluido en aguas poco profundas y del transporte de contaminantes, discretizadas en el tiempo y linealizadas. Las matrices <math display="inline">\mathbf{K}_{ij}</math>, <math display="inline">\mathbf{A}_{i},</math> <math display="inline">\mathbf{S}</math> y <math display="inline">\mathbf{F}</math>, de nuevo para el caso de un contaminante, son

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{K}_{11}=\frac{\mu _{H}}{\rho }\left[\begin{array}{cccc}\frac{4}{3} & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & k_{11}^{*} \end{array}\right],\quad \mathbf{K}_{12}=\frac{\mu _{H}}{\rho }\left[\begin{array}{cccc}0 & \frac{1}{6} & 0 & 0\\ \frac{1}{6} & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & k_{12}^{*} \end{array}\right],\quad \mathbf{K}_{12}=\frac{\mu _{H}}{\rho }\left[\begin{array}{cccc}0 & \frac{1}{6} & 0 & 0\\ \frac{1}{6} & 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & k_{21}^{*} \end{array}\right],</math>
|-
| style="text-align: center;" | <math>   \mathbf{K}_{22}=\frac{\mu _{H}}{\rho }\left[\begin{array}{cccc}1 & 0 & 0 & 0\\ 0 & \frac{4}{3} & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & k_{22}^{*} \end{array}\right],\quad \mathbf{A}_{1}=\left[\begin{array}{cccc}a_{1} & 0 & 1 & 0\\ 0 & a_{1} & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 0 & 0 & a_{1} \end{array}\right],\quad \mathbf{A}_{2}=\left[\begin{array}{cccc}a_{2} & 0 & 0 & 0\\ 0 & a_{2} & 1 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 0 & a_{2} \end{array}\right],</math>
|-
| style="text-align: center;" | <math>   \mathbf{S}=\left[\begin{array}{cccc}\partial _{1}a_{1}+\partial _{2}a_{2}+S_{u,11} & S_{u,12} & 0 & 0\\ S_{u,21} & \partial _{1}a_{1}+\partial _{2}a_{2}+S_{u,22} & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \partial _{1}a_{1}+\partial _{2}a_{2} + S_{\alpha \alpha } h_0^{-1} \end{array}\right],\quad \mathbf{F=}\left[\begin{array}{c}-Q_{1}\\ -Q_{2}\\ \begin{array}{c}0\\ -Q_{4} \end{array} \end{array}\right]. </math>
|}
|}

Obsérvese que las matrices <math display="inline">\mathbf{K}_{ij}</math> son simétricas. Los coeficientes de las matrices <math display="inline">\mathbf{A}_{i},</math> <math display="inline">\mathbf{S}</math> y <math display="inline">\mathbf{F}</math> no todos son constantes, pues algunos dependen de las velocidades del flujo, elevación del nivel del agua, su batimetría y sus gradientes, cuyos valores corresponden a la iteración anterior resultado de la linealización. Asimismo, en el caso de reacciones no lineales entre contaminantes los coeficientes <math display="inline">S_{\alpha \beta }</math> pueden depender de las concentraciones de estos evaluadas en la iteración anterior del actual paso de tiempo.

En el caso de <math display="inline">N</math> contaminantes, el vector de incógnitas es

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{u}^{T}=\left[u_{1},u_{2},p,\varphi _{a,1},\dots ,\varphi _{a,N}\right]. </math>
|}
|}

==3. Aproximación espacial y estabilización numérica==

Una vez hemos conseguido escribir las ecuaciones del problema en la forma ([[#eq-7|7]]), el proceso es exactamente el mismo que el seguido en <span id='citeF-27'></span>[[#cite-27|[27]]] para llevar a cabo la aproximación espacial, por lo que aquí nos centraremos solamente en el problema discreto final y la forma de tratar la ecuación de transporte de contaminantes.

La aproximación de Galerkin del problema es estándar. Dicha aproximación parte de la forma variacional del problema. Una vez discretizado el dominio de cálculo en una partición de elementos finitos <math display="inline">\{ \Omega ^{e}\} </math> de <math display="inline">\Omega </math>, con <math display="inline">e=1,\dots ,n_{{el}}</math>, siendo <math display="inline">n_{{el}}</math> el número de elementos de la partición, podemos construir el espacio de elementos finitos que aproxima el espacio donde se encuentra la solución del problema continuo. El método de Galerkin consiste en tomar las funciones de test de la forma variacional en el mismo espacio.

En el problema que nos ocupa, el método de Gakerkin puede ser inestable por tres motivos. En primer lugar, la interpolación de la altura de la lámina de agua no puede ser independiente de la de la velocidad promediada: ambas interpolaciones tienen que cumplir la condición inf-sup que garantiza que el problema está bien puesto desde el punto de vista matemático. En segundo lugar, si los términos viscosos (o difusivos) son pequeños comparados con los convectivos, pueden aparecer oscilaciones en todo el dominio de cálculo, como es bien conocido. Este problema afecta tanto a las ecuaciones mecánicas, para altura de agua y velocidad promedio, como a las ecuaciones de transporte de contaminantes. Finalmente, si los términos reactivos son dominantes respecto de los difusivos, pueden aparecer también oscilaciones numéricas, esta vez localizadas cerca de los contornos y de las capas límite. Aunque este problema pudiera parecer poco importante, estas oscilaciones localizadas pueden ser muy perjudiciales en problemas no lineales como el que nos ocupa, pues pueden impedir que los esquemas iterativos converjan.

Todas las inestabilidades descritas se pueden evitar usando métodos de elementos finitos estabilizados. Como se ha dicho anteriormente, en este trabajo usaremos métodos basados en el concepto de VMS. En <span id='citeF-17'></span>[[#cite-17|[17]]] se presenta una revisión de los mismos aplicados a problemas de mecánica de fluidos, mientras que en <span id='citeF-27'></span>[[#cite-27|[27]]] se describe su aplicación a un problema formalmente análogo al que nos ocupa, aunque sin las ecuaciones de transporte de contaminantes.

En el problema que estamos considerando, la forma final discreta del problema a resolver es:

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\delta ^n_{t}\mathbf{u}_{h},\mathbf{v}_{h}\right)+B\left(\mathbf{u}_{h},\mathbf{v}_{h}\right)+\underset{e=1}{\overset{n_{{el}}}{\sum }}\int _{\Omega ^{e}}\left(-\mathcal{L}^{*}\left(\mathbf{v}_{h}\right)\right)^{T}\boldsymbol{\tau }^{e}\mathbf{P}\left[\partial _{t}\mathbf{u}_{h}+\mathcal{L}\left(\mathbf{u}_{h}\right)-\mathbf{F}\right]{\,{d}}\Omega =\left\langle \mathbf{F},\mathbf{v}_{h}\right\rangle{.} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

Para explicar los distintos términos de este problema discreto, empecemos diciendo que se entiende que esta ecuación corresponde a un cierto paso de tiempo de la discretización temporal, aunque hemos omitido el superíndice del paso de tiempo para no cargar la notación. Asimismo, hemos visto que estas ecuaciones corresponden también a una cierta iteración del proceso iterativo resultado de la linealización del problema.

El subíndice <math display="inline">h</math> hace referencia a funciones en el espacio de elementos finitos (que aquí nada tiene que ver con la altura de la lámina de agua)). Así, <math display="inline">\mathbf{u}_{h}</math> es la aproximación a la incógnita <math display="inline">\mathbf{u}</math>, mientras que <math display="inline">\mathbf{v}_{h}</math> es una función de test arbitraria, cumpliendo las adecuadas condiciones de contorno. Hemos usado el símbolo <math display="inline">\langle \cdot ,\cdot \rangle </math> para indicar la integral sobre el dominio de cálculo del producto de dos funciones, reemplazándolo por <math display="inline">(\cdot ,\cdot )</math> cuando estas funciones son de cuadrado integrable, como sucede con el término temporal.

El operador <math display="inline">\mathcal{L}</math> es el asociado al problema ([[#eq-7|7]]), es decir

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{L}\left(\mathbf{u}\right):=-\partial _{i}\left(\mathbf{K}_{ij}\partial _{j}\mathbf{u}\right)+\partial _{i}\left(\mathbf{A}_{i}\mathbf{u}\right)+\mathbf{Su}, </math>
|}
|}

mientras que <math display="inline">\mathcal{L}^{*}</math> denota su adjunto (sin condiciones de contorno), dado por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{L}^{*}\left(\mathbf{v}\right):=-\partial _{i}\left(\mathbf{K}^T_{ij}\partial _{j}\mathbf{v}\right)-\mathbf{A}^T_{i}\partial _{i}\mathbf{v}+\mathbf{S}^T\mathbf{v}. </math>
|}
|}

La forma <math display="inline">B</math> en la ecuación ([[#eq-16|16]]) es la asociada al problema escrito en forma débil, es decir, la del método de Galerkin. Suponiendo que las matrices <math display="inline">\mathbf{A}_{i}</math> se descomponen en una parte que no se integra por partes, <math display="inline">\mathbf{A}^c_{i}</math>, y otra que sí, <math display="inline">\mathbf{A}^f_{i}</math> (y que por consiguiente contribuye a los flujos de contorno), y que las condiciones de contornos son homogéneas, para simplificar, en nuestro caso tenemos que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B(\mathbf{u},\mathbf{v})  =\int _{\Omega }\partial _{i}\mathbf{v}^{T}\mathbf{K}_{ij}\partial _{j}\mathbf{u}{\,{d}}\Omega{+\int}_{\Omega }\mathbf{v}^{T}\partial _{i}\left(\mathbf{A}_{i}^{c}\mathbf{u}\right){\,{d}}\Omega{-\int}_{\Omega }\partial _{i}\mathbf{v}^{T}\mathbf{A}_{i}^{f}\mathbf{u}{\,{d}}\Omega{+\int}_{\Omega }\mathbf{v}^{T}\mathbf{Su}{\,{d}}\Omega{.} </math>
|}
|}

Solo nos queda por describir el término llamado de estabilización en la ecuación ([[#eq-16|16]]). En primer lugar, el operador <math display="inline">\mathbf{P}</math> define el método que estamos considerando. En el caso en que <math display="inline">\mathbf{P} = \mathbf{I}</math>, la identidad, nos referiremos a este método como ASGS, mientras que si <math display="inline">\mathbf{P}</math> es la proyección ortogonal al espacio de elementos finitos, el método resultante es el llamado OSS. Los detalles de cómo motivar dichos métodos y cómo se particularizan al problema que nos ocupa (aunque sin las concentraciones de contaminantes), pueden encontrarse en <span id='citeF-17'></span><span id='citeF-27'></span>[[#cite-17|[17,27]]]. Solo falta describir el cálculo de la llamada matriz de parámetros de estabilización <math display="inline">\boldsymbol{\tau }^{e}</math>, calculada en cada elemento <math display="inline">e</math> de la partición de elementos finitos. Tomaremos dicha matriz como diagonal, de la forma

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\tau }^{e} =\mathrm{diag}(\tau _{1}^{e},\tau _{1}^{e},\tau _{2}^{e},\tau _{3}^{e},\dots ,\tau _{2+N}^{e}) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

con

<span id="eq-18"></span>
<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau _{1}^{e}=\left[c_{1}\frac{\nu _{H}}{\left(\frac{h^{e}}{d^{2}}\right)^{2}}+c_{2}\frac{|\boldsymbol{a}|}{\frac{h^{e}}{d}}+c_{3}|S_{u,11}|+c_{4}|S_{u,12}|\right]^{-1},\mbox{ }\tau _{2}^{e}=\frac{\left(\frac{h^{e}}{d}\right)^{2}}{c_{1}\tau _{1}^{e}}, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|-
| style="text-align: center;" | <math>  \tau _{2+\alpha }^{e}=\left[c_{1}\frac{\min \{ k_{11}^{\alpha }, k_{22}^{\alpha }\} }{\left(\frac{h^{e}}{d^{2}}\right)^{2}}+c_{2}\frac{|\boldsymbol{a}|}{\frac{h^{e}}{d}}+c_{3}|S_{\alpha \alpha }| h_0^{-1}\right]^{-1}, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}
|}

donde <math display="inline">c_{1},</math> <math display="inline">c_{2},</math> <math display="inline">c_{3}</math> y <math display="inline">c_{4}</math> son parámetros numéricos, y que para los ejemplos numéricos hemos adoptado <math display="inline">c_{1}=12,</math> <math display="inline">c_{2}=2</math> y <math display="inline">c_{3}=c_{4}=1</math>. En estas expresiones, <math display="inline">\nu _{H}=\frac{\mu _{H}}{\rho }</math> es el coeficiente de viscosidad cinemática, <math display="inline">|\boldsymbol{a}|</math> es la norma de la velocidad y <math display="inline">d</math> es el orden polinomial de interpolación.

La justificación de los parámetros de estabilización usados, así como algunas técnicas numéricas que hemos empleado en elementos de alto orden, puede consultarse en <span id='citeF-17'></span><span id='citeF-30'></span>[[#cite-17|[17,30]]] (véase también <span id='citeF-31'></span><span id='citeF-32'></span>[[#cite-31|[31,32]]]).

==4. Resultados de los experimentos numéricos==

Con el propósito de evaluar las formulaciones de los métodos variacionales de subescalas ASGS y OSS con elementos finitos de alto orden en la ecuación vectorial de CDR transitoria, con no linealidad en los términos convectivo y de reacción, hemos realizado diferentes experimentos numéricos, tales como pruebas de convergencia en malla para el problema del transporte de un contaminante, la aproximación de un ejemplo numérico del transporte de un contaminante en una cavidad cuadrada que se encuentra en la literatura, un ejemplo práctico de distribución del transporte de un contaminante en el golfo de Roses y en la desembocadura del río Guadalquivir, y un ejemplo del modelo depredador-presa.

El método de linealización utilizado es el de Picard. El criterio de convergencia para la no linealidad se ha fijado en la norma <math display="inline">L^{2},</math> es decir, el proceso iterativo se considera convergido cuando

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \frac{\left\Vert \mathbf{u}_{h}^{i}-\mathbf{u}_{h}^{i-1}\right\Vert _{L^{2}}}{\left\Vert \mathbf{u}_{h}^{i}\right\Vert _{L^{2}}}< {tol}, </math>
|}
|}

donde <math display="inline">i</math> es el contador de iteraciones en un paso de tiempo dado y <math display="inline">tol</math> es la tolerancia de convergencia, fijada en <math display="inline">10^{-5}</math> en todos los ejemplos. Aquí y en lo que sigue, las normas <math display="inline">L^{2}</math> se extienden a todo el dominio <math display="inline">\Omega{.}</math> En los casos del transporte de contaminantes estudiaremos el transporte de un solo contaminante, excepto en el ejemplo del modelo depredador-presa.

===4.1 Pruebas de convergencia en malla===

Las pruebas de convergencia en malla para el problema del transporte de un contaminante que hemos realizado consisten en calcular el error en norma <math display="inline">L^{2}</math> entre el valor exacto de una solución manufacturada y la aproximación numérica de una de las incógnitas del problema, el cual debe ser menor o igual que el valor teórico. Sabemos que si el error del método se comporta como el error de la interpolación, situación que puede considerarse óptima, este error debe cumplir con

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>e=\left\Vert u-u_{h}\right\Vert _{L^{2}}\leqslant Ch^{d+1}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

donde <math display="inline">u_{h}</math> es la aproximación de la incógnita <math display="inline">u</math> en el espacio de elementos finitos, <math display="inline">C</math> es una constante positiva, <math display="inline">h</math> es el diámetro de la partición de elementos finitos (no confundir con la altura de la lámina de agua) y <math display="inline">d</math> el grado polinómico de la función de forma correspondiente a la interpolación de la aproximación <math display="inline">u_{h}</math>. Consideraremos pruebas de convergencia en malla para <math display="inline">1\leq d\leq{4}</math>.

Si graficamos la ecuación ([[#eq-20|20]]) (en el caso de la igualdad) en el plano, con <math display="inline">\mathrm{log}\,e</math> en las ordenadas y <math display="inline">\mathrm{log}\,h</math> en las abscisas, tenemos una línea recta, donde <math display="inline">d+1</math> es la pendiente teórica óptima de convergencia. Las pruebas de convergencia en malla consisten entonces en verificar que la pendiente calculada sea precisamente <math display="inline">d+1</math>. Para ello, seleccionamos como solución exacta tanto para las velocidades <math display="inline">U_1</math>, <math display="inline">U_2</math> como para la elevación de la superficie libre del agua <math display="inline">\eta </math> y para la concentración de un contaminante <math display="inline">\varphi </math>, la función polinómica espacio-temporal

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \eta (x,y,t)=U_{1}(x,y,t)=U_{2}(x,y,t)=\varphi (x,y,t)=x^{6}y^{6}(1-x^{6})(1-y^{6})t, </math>
|}
|}

e imponemos el término fuente <math display="inline">\mathbf{F}</math> en la ecuación vectorial de CDR transitoria ([[#eq-7|7]]) mediante las siguientes expresiones tomadas de las ecuaciones ([[#eq-1|1]]), ([[#eq-2|2]]) y ([[#eq-3|3]]):

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f_{i}=  \partial _{t}\left(hU_{i}\right)+\partial _{j}\left(hU_{i}U_{j}\right)+\partial _{i}\frac{1}{2}g\left(h^{2}-H^{2}\right)-\partial _{j}\left(\frac{h}{\rho }\bar{\tau }_{ji}\right)+Q_{i}\qquad i=1,2,</math>
|-
| style="text-align: center;" | <math>  f_{3} =  \partial _{t}h+\partial _{i}\left(hU_{i}\right),</math>
|-
| style="text-align: center;" | <math>  f_{4} =  \partial _{t}\left(h\varphi \right)+\partial _{i}\left(hU_{i}\varphi \right)-\partial _{i}\left(hk_{ij}\partial _{j}\varphi \right), </math>
|}
|}

siendo <math display="inline">\mathbf{F}=[f_{1},f_{2},f_{3},f_{4}]^{T}</math> y <math display="inline">\bar{\tau }_{ji}</math> y <math display="inline">Q_{i}</math> dados por las ecuaciones ([[#eq-4|4]]) y ([[#eq-5|5]]), respectivamente. Hemos tomado batimetría constante con profundidad de 1 m, es decir, <math display="inline">H=1</math> m, aceleración de la gravedad <math display="inline">g=10</math> <math display="inline">\mathrm{m/s^{2}}</math>, viscosidad cinemática <math display="inline">\nu _{H}=10^{-3}</math> <math display="inline">\mathrm{m^{2}/s}</math> y los coeficientes empíricos de dispersión <math display="inline">k_{11}=k_{22}=10^{-3}</math> <math display="inline">\mathrm{m^{2}/s},</math> <math display="inline">k_{12}=k_{21}=0.</math> El parámetro de Coriolis <math display="inline">\hat{f},</math> las tensiones de superficie libre debidas al viento <math display="inline">\tau _{3i}^{s},</math> las variaciones de presión atmosférica <math display="inline">p_{a}</math> y el coeficiente de fricción en el fondo Ch no han sido tomados en cuenta para los ejemplos de las pruebas de convergencia en malla, como tampoco el término reactivo de la ecuación de transporte de contaminantes.

A continuación detallamos los resultados de los experimentos numéricos de las pruebas de convergencia en malla con solución analítica conocida, tanto para el método ASGS como para el OSS. En cada gráfica se puede apreciar con línea continua la pendiente teórica de convergencia y con línea con apéndices sobre ella la pendiente calculada. Presentamos los resultados para los distintos grados polinómicos de las funciones de forma estudiados, es decir, elementos lineales, cuadráticos, cúbicos y de cuarto orden triangulares <math display="inline">P_{1},\;P_{2},\;P_{3},\;P_{4}</math> y elementos cuadrangulares <math display="inline">Q_{1},\;Q_{2},\;Q_{3},\;Q_{4}</math>, respectivamente.

El dominio computacional es el cuadrado <math display="inline">\left[0,1\right]\times \left[0,1\right]</math> y el intervalo de tiempo es <math display="inline">\left[0,1\right]</math>, con incrementos para cada paso de tiempo <math display="inline">\delta t=0.2</math>. La malla de elementos finitos consiste en triángulos o cuadrados formando una malla regular. Para cada grado polinómico de las funciones de forma <math display="inline">P_{1},\;P_{2},\;P_{3},\;P_{4}</math> o <math display="inline">Q_{1},\;Q_{2},\;Q_{3},\;Q_{4},</math> hemos calculado los errores con norma <math display="inline">L^{2}</math> para ocho mallas de tamaño <math display="inline">h</math>, siendo <math display="inline">\frac{1}{h}=15,\;20,\;25,\;30,\;35,\;40,\;45,\;50,</math> el número de partes en que se ha dividido cada lado del cuadrado. En la  [[#table-1|Tabla 1]] se muestran el número de elementos y el número de nodos para cada refinamiento, tanto para elementos triangulares como para elementos cuadrangulares. Los valores de las constantes algorítmcas <math display="inline">c_{i}</math>, <math display="inline">i=1,2,3,4,</math> se han calibrado a <math display="inline">c_{1}=12,</math> <math display="inline">c_{2}=2,</math> <math display="inline">c_{3}=1,</math> <math display="inline">c_{4}=1</math> y el integrador temporal utilizado es BDF3 (diferencias finitas hacia atrás de cuatro niveles de tiempo y tercer orden en <math display="inline">\delta t</math>).

<div class="center" style="font-size: 75%;">'''Tabla 1'''. Refinamiento para elementos triangulares y cuadrangulares</div>

<div id='table-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! colspan='2' | Triángulos !!  colspan='4' | Número de nodos
|- style="text-align:center"
| <math display="inline">1/h</math> 
| elem. 
| <math>P1</math>
| <math>P2</math>
| <math>P3</math>
| <math>P4</math>
|- style="text-align:center"
|  15 
| 450 
| 256 
| 961 
| 2116 
| 3721
|-style="text-align:center"
| 20 
| 800 
| 441 
| 1681 
| 3721 
| 6561
|-style="text-align:center"
| 25 
| 1250 
| 676 
| 2601 
| 5776 
| 10201
|-style="text-align:center"
| 30 
| 1800 
| 961 
| 3721 
| 8281 
| 14641
|-style="text-align:center"
| 35 
| 2450 
| 1296 
| 5041 
| 11236 
| 19881
|-style="text-align:center"
| 40 
| 3200 
| 1681 
| 6561 
| 14641 
| 25921
|-style="text-align:center"
| 45 
| 4050 
| 2116 
| 8281 
| 18496 
| 32761
|- style="text-align:center"
| 50 
| 5000 
| 2601 
| 10201 
| 22801 
| 40401
|}
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! colspan='2' | Cuadrados !!  colspan='4' | Número de nodos
|- style="text-align:center"
| <math display="inline">1/h</math> 
| elem. 
| <math>Q1</math>
| <math>Q2</math>
| <math>Q3</math>
| <math>Q4</math>
|- style="text-align:center"
|  15 
| 225 
| 256 
| 961 
| 2116 
| 3721
|-style="text-align:center"
| 20 
| 400 
| 441 
| 1681 
| 3721 
| 6561
|-style="text-align:center"
| 25 
| 625 
| 676 
| 2601 
| 5776 
| 10201
|-style="text-align:center"
| 30 
| 900 
| 961 
| 3721 
| 8281 
| 14641
|-style="text-align:center"
| 35 
| 1250 
| 1296 
| 5041 
| 11236 
| 19881
|-style="text-align:center"
| 40 
| 1600 
| 1681 
| 6561 
| 14641 
| 25921
|-style="text-align:center"
| 45 
| 2025 
| 2116 
| 8281 
| 18496 
| 32761
|- style="text-align:center"
| 50 
| 2500 
| 2601 
| 10201 
| 22801 
| 40401
|}


En el encabezado de las gráficas de las pruebas de convergencia en malla de las Figuras [[#img-1|1]] y [[#img-2|2]] se adjuntan dos filas con los valores de las pendientes calculadas; la primera fila corresponde a los valores de las pendientes de las rectas que pasan por los primeros 5 puntos y la segunda fila son las pendientes de las rectas que pasan por los últimos 5 puntos de los 8 puntos del refinamiento de malla evaluados <math display="inline">\left(\frac{1}{h}=15,\;20,\;25,\;30,\;35,\;40,\;45,\;50\right)</math>. De esta manera podemos visualizar numéricamente la tendencia de la convergencia de las pendientes calculadas hacia el valor teórico. Para el cálculo de dichas pendientes se ha utilizado el método de los mínimos cuadrados; la simbología <math display="inline">m(P1),</math> <math display="inline">m(P2),</math> <math display="inline">m(P3),</math> <math display="inline">m(P4),</math> corresponde a las pendientes para elementos triangulares: lineales, cuadráticos, cúbicos y de cuarto orden respectivamente, mientras que para elementos cuadrangulares lineales, cuadráticos, cúbicos y de cuarto orden la simbología es <math display="inline">m(Q1),</math> <math display="inline">m(Q2),</math> <math display="inline">m(Q3),</math> <math display="inline">m(Q4)</math>. Los valores de las pendientes teóricas para las variables de la velocidad y de concentración de un contaminante son <math display="inline">m(P1)=m(Q1)=2</math>, <math display="inline">m(P2)=m(Q2)=3</math>, <math display="inline">m(P3)=m(Q3)=4</math>, <math display="inline">m(P4)=m(Q4)=5</math>, y las pendientes teóricas para la variable <math display="inline">\eta </math>, que es la elevación de la superficie libre del agua, son <math display="inline">m(P1)=m(Q1)=1</math>, <math display="inline">m(P2)=m(Q2)=2</math>, <math display="inline">m(P3)=m(Q3)=3</math>, <math display="inline">m(P4)=m(Q4)=4</math>.

En la  [[#img-1|Figura 1]], se muestra la convergencia en malla para la primera componente de la velocidad usando elementos triangulares y en la  [[#img-2|Figura 2]] la convergencia en malla para la concentración de un contaminante. Las otras variables muestran igualmente convergencia óptima para todos los elementos considerados. Lo mismo sucede con los elementos cuadrangulares.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: auto;"
|-
|style="padding-top:10px;" | [[Image:Draft_Codina_913592092-test-ASGS_U1.png|600px|]]
|style="padding-top:10px;" | [[Image:Draft_Codina_913592092-test-OSS_U1.png|600px|Convergencia en malla ASGS-OSS, elementos triangulares, velocidad U₁.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" style="padding-bottom:10px;" |  '''Figura 1'''. Convergencia en malla ASGS-OSS, elementos triangulares, velocidad <math>U_1</math>
|}

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: auto;"
|-
|style="padding-top:10px;" | [[Image:Draft_Codina_913592092-test-ASGS_fi.png|600px|]]
|style="padding-top:10px;" | [[Image:Draft_Codina_913592092-test-OSS_fi.png|600px|Convergencia en malla ASGS-OSS, elementos triangulares, concentración φ.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" style="padding-bottom:10px;" |  '''Figura 2'''. Convergencia en malla ASGS-OSS, elementos triangulares, concentración <math>\varphi </math>
|}

===4.2 Transporte de un contaminante en una cavidad cuadrada===

Este caso de referencia ha sido ampliamente examinado en la literatura con diferentes métodos, en <span id='citeF-6'></span>[[#cite-6|[6]]] con el método de diferencias finitas, con el método de volúmenes finitos en <span id='citeF-3'></span><span id='citeF-5'></span>[[#cite-3|[3,5]]], con el método SUPG en <span id='citeF-4'></span>[[#cite-4|[4]]]; en todos los casos se ha considerado convección pura, es decir, el coeficiente empírico de dispersión del transporte de un contaminante se ha tomado nulo. En el presente trabajo vamos a reexaminar el caso de referencia con el método variacional ASGS y utilizando funciones de forma cuadráticas, cúbicas y de cuarto orden, en convección dominante, es decir, considerando el coeficiente empírico de dispersión <math display="inline">k_{11}=k_{22}=10^{-3}</math> <math display="inline">\mathrm{m^{2}/s}</math> y <math display="inline">k_{12}=k_{21}=0</math> <math display="inline">\mathrm{m^{2}/s}</math>.

El dominio computacional es un canal cuadrado de 9 km <math display="inline">\times 9</math> km. En cuanto al mallado, con el objeto de comparar los resultados entre elementos lineales y de alto orden, hemos considerado conveniente utilizar refinamientos de malla para elementos <math display="inline">Q_{1},</math> <math display="inline">Q_{2},</math> <math display="inline">Q_{3}</math> y <math display="inline">Q_{4}</math> que den el mismo número de nodos entre si. Tomamos como referencia el refinamiento de malla para <math display="inline">Q_{2}</math>, que consta de <math display="inline">90 \times 90</math> elementos cuadrados uniformes de 9 nodos, es decir, el espaciado en las direcciones <math display="inline">x</math> e <math display="inline">y</math> es <math display="inline">\delta x=\delta y=100</math> <math display="inline">\mathrm{m}</math>, con un total de 8100 elementos y 32761 nodos. De esta manera, la malla de referencia para elementos <math display="inline">Q_{1}</math> consta de 180<math display="inline">\times </math>180 elementos cuadrados uniformes de 4 nodos, es decir <math display="inline">\delta x=\delta y=50</math> <math display="inline">\mathrm{m}</math>, con un total de 32400 elementos y 32761 nodos. La malla de referencia para elementos <math display="inline">Q_{3}</math> consta de <math display="inline">60 \times 60</math> elementos cuadrados uniformes de 16 nodos, es decir <math display="inline">\delta x=\delta y=150</math> <math display="inline">\mathrm{m}</math>, con un total de 3600 elementos y 32761 nodos. Finalmente, la malla de referencia para elementos <math display="inline">Q_{4}</math> consta de 45<math display="inline">\times </math>45 elementos cuadrados uniformes de 25 nodos, es decir <math display="inline">\delta x=\delta y=200</math> <math display="inline">\mathrm{m}</math>, con un total de 2025 elementos y 32761 nodos.

El tamaño del paso de tiempo es uniforme, con <math display="inline">\delta t=20</math> <math display="inline">\mathrm{s}</math>, y el integrador temporal es BDF3 en todos los casos.

Un flujo uniforme con <math display="inline">u_1=u_2=0.5</math> <math display="inline">\mathrm{m/s}</math> y <math display="inline">H=0.2485</math> <math display="inline">\mathrm{m}</math> está impuesto en todo el dominio.  La concentración inicial de contaminante está dada por la siguiente suposición de dos distribuciones Gausianas centradas respectivamente en <math display="inline">x_{1}=y_{1}=1400</math> <math display="inline">\mathrm{m}</math> y en <math display="inline">x_{2}=y_{2}=2400</math> <math display="inline">\mathrm{m}</math>, dada por:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \varphi (x,y, t = 0)=\varphi _{1}\exp{\left(\frac{\left(x-x_{1}\right)^{2}+\left(y-y_{1}\right)^{2}}{\sigma _{1}^{2}}\right)}+\varphi _{2}\exp{\left(\frac{\left(x-x_{2}\right)^{2}+\left(y-y_{2}\right)^{2}}{\sigma _{2}^{2}}\right),} </math>
|}
|}

donde <math display="inline">\varphi _{1}=10</math>, <math display="inline">\varphi _{2}=6.5</math> y <math display="inline">\sigma _{1}=\sigma _{2}=264</math>.

Dada la baja difusión del problema, la solución exacta es aquella concentración de contaminante que se mueve diagonalmente a través del dominio con velocidad constante, preservándose su forma durante todo el tiempo. La Figura [[#img-3|3]] ilustra los resultados usando el elemento <math display="inline">Q_{2}</math> comparados con la solución teórica en diferentes tiempos, observándose en todos los casos resultados completamente satisfactorios comparados con la solución exacta.

En la [[#table-2|Tabla 2]], presentamos los valores máximos y mínimos de concentración de contaminante para las funciones de forma <math display="inline">Q_{1}</math>, <math display="inline">Q_{2}</math>, <math display="inline">Q_{3}</math> y <math display="inline">Q_{4}</math>, en donde observamos que a partir de <math display="inline">Q_{2}</math> el error de aproximación al valor exacto máximo es del orden de las centésimas, por lo cual en la [[#img-3|Figura 3]] hemos usado los resultados del elemento <math display="inline">Q_{2}</math>.

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"
|-
|style="padding:10px;" | [[Image:Draft_Codina_913592092-test-1600_3D.png|600px|]]
|style="padding:10px;" | [[Image:Draft_Codina_913592092-test-1600_corte_x.png|480px|]]
|-
|[[Image:Draft_Codina_913592092-test-4800_3D.png|600px|]]
|[[Image:Draft_Codina_913592092-test-4800_corte_x.png|480px|]]
|-
|style="padding:10px;" | [[Image:Draft_Codina_913592092-test-9600_3D.png|600px|]]
|style="padding:10px;" | [[Image:Draft_Codina_913592092-test-9600_corte_x.png|480px|Contornos 3D (izquierda) y distribución del contaminante (derecha) a lo largo de la diagonal de la cavidad cuadrada, con elementos Q₂.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" style="padding:10px;" |  '''Figura 3'''. Contornos 3D (izquierda) y distribución del contaminante (derecha) a lo largo de la diagonal de la cavidad cuadrada, con elementos <math>Q_2</math>
|}


<div class="center" style="font-size: 75%;">'''Tabla 2'''. Valores máximos y mínimos de concentración de contaminante en la cavidad cuadrada para <math>t=9600</math>&nbsp;s</div>

<div id='tab-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
! !!  Solución exacta !!  <math>Q_{1}</math> !!  <math>Q_{2}</math> !!  <math>Q_{3}</math> !!  <math>Q_{4}</math>
|- style="text-align:center"
|  style="text-align:left;" |  Número de elementos 
|  -
|  32400 
|  8100 
|  3600 
|  2025
|- style="text-align:center"
|   style="text-align:left;" |Número de nodos 
|  -
|  32761 
|  32761 
|  32761 
|  32761
|- style="text-align:center"
|   style="text-align:left;" |Máximo valor de <math display="inline">\varphi </math> 
|  10 
|  9.312 
|  9.976 
|  10.03 
|  10.04
|- style="text-align:center"
|  style="text-align:left;" | Mínimo valor de <math display="inline">\varphi </math> 
|  0 
|  <math>-0.01032</math>
|  <math>-3.702\times{10}^{-10}</math>
|  <math>-2.986\times{10}^{-5}</math>
|  <math>-3.533\times{10}^{-6}</math>
|}

===4.3 Transporte de un contaminante en el golfo de Roses y en la desembocadura del río Guadalquivir===

A continuación presentamos dos ejemplos numéricos del transporte de un contaminante. Ambos ejemplos corresponden a situaciones físicas reales, con geometrías complejas. El objetivo es mostrar que el modelo numérico propuesto puede ser utilizado en este tipo de situaciones, aunque no disponemos de valores experimentales con los cuales comparar los resultados y por consiguiente los mostraremos solamente de manera cualitativa, sin entrar en su discusión. Desde el punto de vista numérico, el valor de estos resultados es comprobar que no presentan ningún tipo de oscilación numérica, pese a utilizar igual interpolación para todas las variables y tratarse de situaciones de flujo con convección dominante.

Consideramos una concentración constante de un solo contaminante <math display="inline">\varphi=1</math> en el flujo de entrada y coeficientes de difusión constantes <math display="inline">k_{11}=k_{22}=10^{-3}</math> <math display="inline">\mathrm{m^{2}/s}</math> y <math display="inline">k_{12}=k_{21}=0</math> <math display="inline">\mathrm{m^{2}/s}</math>. La suposición de que los coeficientes de difusión son constantes es debido a que los efectos de turbulencia no han sido considerados en el modelo. También asumimos que la concentración del contaminante se anula en las fronteras sólidas, por lo que prescribimos a cero dicha concentración  tanto en el perfil de las costas como en las riveras del río, en los dos problemas tratados.  En el resto de los contornos consideramos la concentración de contaminante libre, con concentración inicial igual a cero. En ambos ejemplos numéricos hemos considerado batimetría constante con una profundidad de <math display="inline">H=20.0</math>&nbsp;m (no hemos considerado pues la disminución de la profundidad hasta <math display="inline">H=0</math>), la aceleración de la gravedad <math display="inline">g=10</math> m/s<math display="inline">^{2}</math> y la viscosidad cinemática <math display="inline">\nu _{H}=10^{-6}</math> m<math display="inline">^{2}</math>/s. El parámetro de Coriolis <math display="inline">\hat{f}</math>, las tensiones de superficie libre debidas al viento <math display="inline">\tau _{3i}^{s}</math>, las variaciones de presión atmosférica <math display="inline">p_{a}</math> y el coeficiente de fricción en el fondo <math display="inline">\mbox{Ch}</math> no han sido tomados en cuenta en ninguno de los dos ejemplos.

El primer ejemplo es la simulación del transporte de un contaminante de concentración constante <math display="inline">\varphi=1</math>, debido a una corriente marina que viene desde el norte en el golfo de Roses sobre la costa catalana, al sur del Cabo de Creus. La geometría del problema y los resultados numéricos se muestran en las Figuras [[#img-4|4]], [[#img-5|5]] y [[#img-6|6]]. El dominio de simulación tiene dimensiones de aproximadamente 80 <math display="inline">\times </math> 40 km<math display="inline">^{2}</math>. La malla de elementos finitos consiste en 6062 elementos triangulares lineales con 3217 nodos.

Una corriente marina de 0.1 m/s que llega desde el norte ha sido prescrita en el borde superior del dominio. Esto corresponde a una velocidad como flujo de entrada en la parte superior del dominio. En la costa la velocidad se prescribe a cero, mientras que en los contornos inferior y derecho del dominio se deja libre, donde la sobreelevación del agua se prescribe a cero.

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 90%;"
|-
|style="padding:10px;" | [[Image:Draft_Codina_913592092-test-malla.png|600px|Golfo de Roses. Geometría y malla de elementos finitos.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;-bottom" |  '''Figura 4'''. Golfo de Roses. Geometría y malla de elementos finitos
|}


<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|-
|  <math>t</math>=50000 s
| style="padding:10px;" | [[Image:Draft_Codina_913592092-test-velocidad_50000.png|300px|Creus/velocidad_50000]]
| style="padding:10px;" | [[Image:Draft_Codina_913592092-test-elevacion_50000.png|300px|Creus/elevacion_50000]]
|-
| <math>t</math>=100000 s
| [[Image:Draft_Codina_913592092-test-velocidad_100000.png|300px|Creus/velocidad_100000]]
| [[Image:Draft_Codina_913592092-test-elevacion_100000.png|300px|Creus/elevacion_100000]]
|-
| <math>t</math>=150000 s
| style="padding:10px;" | [[Image:Draft_Codina_913592092-test-velocidad_50000.png|300px|Creus/velocidad_50000]]
| style="padding:10px;" | [[Image:Draft_Codina_913592092-test-elevacion_150000.png|300px|Creus/elevacion_150000]]
|-
| <math>t</math>=200000 s
| [[Image:Draft_Codina_913592092-test-velocidad_200000.png|300px|Creus/velocidad_200000]]
| [[Image:Draft_Codina_913592092-test-elevacion_200000.png|300px|Creus/elevacion_200000]]
|-
| 
| style="font-size:75%;" | (a)  
| style="font-size:75%;" |(b) 
|}
|-
|- style="text-align: center; font-size: 75%;"
| colspan="3" style="padding-bottom:10px;" |  '''Figura 5'''. Flujo alrededor del golfo de Roses. (a) velocidad (máx: 0.29 m/s), (b) elevación (INCOG3), superficie libre del agua (máx: 2.52 mm, mín: -4.05 mm)
|}


<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|-
| <math>t</math>=50000 s
| style="padding:10px;" | [[Image:Draft_Codina_913592092-test-vector_velocidad_50000.png|300px|Creus/vector_velocidad_50000]]
| style="padding:10px;" | [[Image:Draft_Codina_913592092-test-contaminanate_50000.png|300px|Creus/contaminanate_50000]]
|-
| <math>t</math>=100000 s
| [[Image:Draft_Codina_913592092-test-vector_velocidad_100000.png|300px|Creus/vector_velocidad_100000]]
| [[Image:Draft_Codina_913592092-test-contaminanate_100000.png|300px|Creus/contaminanate_100000]]
|-
|  <math>t</math>=150000 s
| style="padding:10px;" |[[Image:Draft_Codina_913592092-test-vector_velocidad_150000.png|300px|Creus/vector_velocidad_150000]]
| style="padding:10px;" | [[Image:Draft_Codina_913592092-test-contaminanate_150000.png|300px|Creus/contaminanate_150000]]
|-
| <math>t</math>=200000 s
| [[Image:Draft_Codina_913592092-test-vector_velocidad_200000.png|300px|Creus/vector_velocidad_200000]]
| [[Image:Draft_Codina_913592092-test-contaminanate_200000.png|300px|Creus/contaminanate_200000]]
|-
| 
| style="font-size:75%;" |(a) 
| style="font-size:75%;" |(b) 
|}
|-
|- style="text-align: center; font-size: 75%;"
| colspan="3" style="padding-bottom:10px;" | '''Figura 6'''. Flujo alrededor del golfo de Roses. (a) vectores de velocidad, (b) concentración del transporte de un contaminante (INGOG 4).
|}

En la segunda simulación numérica consideramos el transporte de un contaminante de concentración constante <math display="inline">\varphi=1</math> en la desembocadura del río Guadalquivir, en la costa del sur de España. La geometría del problema y los resultados numéricos se muestran en las Figuras [[#img-7|7]], [[#img-8|8]] y [[#img-9|9]].

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 90%;"
|-
|style="padding:10px;" |[[Image:Draft_Codina_913592092-test-malla-dup1.png|600px|Desembocadura del río Guadalquivir. Geometría y malla de elementos finitos]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;" | '''Figura 7'''. Desembocadura del río Guadalquivir. Geometría y malla de elementos finitos
|}

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|-
| <math>t</math>=50000 s
| style="padding:10px;" |[[Image:Draft_Codina_913592092-dup1.png|300px|Guadalquivir/velocidad_50000]]
| style="padding:10px;" |[[Image:Draft_Codina_913592092-dup1.png|300px|Guadalquivir/elevacion_50000]]
|-
| <math>t</math>=100000 s
| [[Image:Draft_Codina_913592092-dup1.png|300px|Guadalquivir/velocidad_100000]]
| [[Image:Draft_Codina_913592092-dup1.png|300px|Guadalquivir/elevacion_100000]]
|-
| <math>t</math>=150000 s
| style="padding:10px;" |[[Image:Draft_Codina_913592092-test-velocidad_150000.png|300px|Guadalquivir/velocidad_150000]]
| style="padding:10px;" |[[Image:Draft_Codina_913592092-dup1.png|300px|Guadalquivir/elevacion_150000]]
|-
| <math>t</math>=200000 s
| [[Image:Draft_Codina_913592092-dup1.png|300px|Guadalquivir/velocidad_200000]]
| [[Image:Draft_Codina_913592092-dup1.png|300px|Guadalquivir/elevacion_200000]]
|-
| 
| style="font-size:75%;" |(a) 
| style="font-size:75%;" |(b) 
|}
|-
|- style="text-align: center; font-size: 75%;"
| colspan="3" style="padding-bottom:10px;" | '''Figura 8'''. Flujo en la desembocadura del Guadalquivir. (a) velocidad (máx: 0.83 m/s), (b) elevación (INCOG3), superficie libre del agua (máx: 1.97 cm, mín: -1.42 cm)
|}

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|-
| <math>t</math>=50000 s
| style="padding:5px;" |[[Image:Draft_Codina_913592092-test-vector_velocidad_50002.png|300px|Guadalquivir/vector_velocidad_50002]]
| style="padding:5px;" |[[Image:Draft_Codina_913592092-test-contaminante_50000.png|300px|Guadalquivir/contaminante_50000]]
|-
| <math>t</math>=100000 s
| [[Image:Draft_Codina_913592092-test-vector_velocidad_100002.png|300px|Guadalquivir/vector_velocidad_100002]]
| [[Image:Draft_Codina_913592092-test-contaminante_100000.png|300px|Guadalquivir/contaminante_100000]]
|-
| <math>t</math>=150000 s
| [[Image:Draft_Codina_913592092-test-vector_velocidad_150002.png|300px|Guadalquivir/vector_velocidad_150002]]
| [[Image:Draft_Codina_913592092-test-contaminante_150000.png|300px|Guadalquivir/contaminante_150000]]
|-
| <math>t</math>=200000 s
| [[Image:Draft_Codina_913592092-test-vector_velocidad_200002.png|300px|Guadalquivir/vector_velocidad_200002]]
| [[Image:Draft_Codina_913592092-test-contaminante_200000.png|300px|Guadalquivir/contaminante_200000]]
|-
| 
| style="font-size:75%;" |(a) 
| style="font-size:75%;" |(b)
|}
|-
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding-bottom:10px;" | '''Figura 9'''. Flujo en la desembocadura del Guadalquivir. (a) vector de velocidad, (b) concentración del transporte de un contaminante (INGOG 4)
|}


La malla de elementos finitos consiste en 4955 elementos triangulares lineales con 2645 nodos.

Un flujo de corriente marina de <math display="inline">0.1</math> m/s está prescrito en el contorno inferior en dirección este-oeste y en la desembocadura del río se prescribe una corriente del río de 0.3 m/s, igualmente en dirección este-oeste. En el resto de la costa y en las riveras del río la velocidad se prescribe a cero, y se deja libre en los contornos izquierdo y superior, donde la sobreelevación del agua está prescrita a cero.

Los dos ejemplos numéricos considerados presentan flujos complejos, con importantes variaciones temporales y espaciales. Esto hace que los parámetros de estabilización sean muy variables elemento a elemento. Este hecho dificulta la convergencia de los esquemas iterativos. Puesto que <math display="inline">\tau _1</math> y <math display="inline">\tau _{3}</math> (en el caso de <math display="inline">N=1</math> contaminantes) son escalas de tiempo, a menudo se toman proporcionales a <math display="inline">\delta t</math>, y esto se puede justificar con argumentos diversos <span id='citeF-24'></span>[[#cite-24|[24]]]. En los problemas en los que la convergencia del esquema iterativo es costosa, hemos tomado una cota máxima para <math display="inline">\tau _1</math> y <math display="inline">\tau _3</math> proporcional a <math display="inline">\delta t</math>, de manera que los hemos redefinido en cada elemento <math display="inline">e</math> como

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau _{1}^{e}\leftarrow \min (\tau _{1}^{e},\alpha \delta t),\quad \tau _{3}^{e}\leftarrow \min (\tau _{3}^{e},\beta \delta t), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

donde <math display="inline">\alpha </math> y <math display="inline">\beta </math> son parámetros que hay que ajustar dependiendo del tamaño de la malla y del grado del polinomio de las funciones de forma. El valor de <math display="inline">\tau _{2}^{e}</math> se sigue calculando mediante la expresión dada en ([[#eq-18|18]]).

Mediante experimentación numérica hemos calibrado los parámetros <math display="inline">\alpha </math> y <math display="inline">\beta </math>. Para el primer ejemplo hemos considerado <math display="inline">\alpha=15</math> y <math display="inline">\beta=6</math> con un paso de tiempo uniforme <math display="inline">\delta t=</math>100 s y para el segundo ejemplo <math display="inline">\alpha=3</math> y <math display="inline">\beta=5</math> con <math display="inline">\delta t=</math>2500 s. En ambos casos hemos usado el integrador temporal BDF2 y presentamos resultados numéricos en diferentes instantes de tiempos de simulación, <math display="inline">t=50000,\;100000,\;150000</math> y <math display="inline">200000</math>&nbsp;s.

===4.4 Modelo depredador-presa===

Vamos a considerar como problema a resolver el modelo depredador-presa mostrado en <span id='citeF-33'></span>[[#cite-33|[33]]], que es un sistema acoplado de ecuaciones transitorias de convección-difusión-reacción no lineales en el término de reacción y que describen espacial y temporalmente la dinámica poblacional de una presa y un depredador en términos de funciones de densidad continuas. Usaremos este problema como ejemplo de sistema con términos reactivos no lineales, con dos concentraciones que compiten entre ellas (en este caso no necesariamente contaminantes). Nos centraremos en el problema propiamente del modelo depredador-presa, suponiendo que la velocidad es conocida. Si embargo, la extensión al caso en el que la velocidad resulta de resolver las ecuaciones del modelo de aguas someras es inmediata. En este problema además es frecuente considerar la posibilidad de que la velocidad de transporte del depredador sea distinta a la de la presa, por lo que cuando introduzcamos estas velocidades de transporte mantendremos esta posible generalización.

====4.4.1 Formulación del modelo depredador-pressa====

La forma general del modelo es como sigue:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _{t}\varphi _{1}=k_{1}(\partial _{1}^{2}\varphi _{1}+\partial _{2}^{2}\varphi _{1})+P(\varphi _{1})-E(\varphi _{1},\varphi _{2}), </math>
|-
| style="text-align: center;" | <math>  \partial _{t}\varphi _{2}=k_{2}(\partial _{1}^{2}\varphi _{2}+\partial _{2}^{2}\varphi _{2})+\kappa E(\varphi _{1},\varphi _{2})-\mu (\varphi _{2}), </math>
|}
|}

donde <math display="inline">\varphi _{1}</math> y <math display="inline">\varphi _{2}</math> representan las densidades de población de la presa y del depredador, respectivamente. <math display="inline">k_{1}</math> y <math display="inline">k_{2}</math> son los coeficientes de difusión. <math display="inline">P(\varphi _{1})</math> es la función de crecimiento de la población de la presa. La función <math display="inline">E(\varphi _{1},\varphi _{2})</math> representa el acto de depredación resultante en una disminución de la población de la presa y en un incremento en la población del depredador. <math display="inline">\kappa </math> es la eficiencia de la depredación, la cual determina el incremento de la población del depredador, siendo <math display="inline">0\leq \kappa \leq{1}</math>. Finalmente, <math display="inline">\mu (\varphi _{2})</math> es la función que determina la cantidad de mortalidad del depredador.

Una de las posibilidades para modelar las funciones <math display="inline">P</math>, <math display="inline">E</math> y <math display="inline">\mu </math> es:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _{t}\varphi _{1}=k_{1}(\partial _{1}^{2}\varphi _{1}+\partial _{2}^{2}\varphi _{1})+C\varphi _{1}\left(1-\frac{\varphi _{1}}{K}\right)-B\left(\frac{\varphi _{1}\varphi _{2}}{\varphi _{1}+H}\right)</math>
|-
| style="text-align: center;" | <math>  \partial _{t}\varphi _{2}=k_{2}(\partial _{1}^{2}\varphi _{2}+\partial _{2}^{2}\varphi _{2})+\kappa B\left(\frac{\varphi _{1}\varphi _{2}}{\varphi _{1}+H}\right)-M\varphi _{2}, </math>
|}
|}

donde la función <math display="inline">E(\varphi _{1},\varphi _{2})</math> está modelada por el llamado modelo de Hollinger tipo II, en el cual <math display="inline">B</math> representa la tasa de depredación y <math display="inline">H</math> es la densidad media de saturación de la presa. Para modelar el crecimiento de la población de la presa usamos el llamado modelo logístico, donde <math display="inline">C</math> representa la tasa de crecimiento de la población de la presa y <math display="inline">K</math> es la capacidad de carga del sistema, y denota la cantidad máxima de población de la presa que es soportada por el dominio; en este trabajo asumimos el valor de la unidad para este parámetro. La función de mortalidad <math display="inline">\mu </math> del depredador está dada por un término lineal, siendo <math display="inline">M</math> su coeficiente.

El modelo anterior considera que no hay transporte convectivo ni del depredador ni de la presa. Este transporte puede deberse a que depredador y presa se encuentren en un medio fluido en movimiento, que se correspondería con el modelo general que hemos tratado anteriormente, pero también puede ser necesario introducir la convección a través de un patrón de migración estacional o migración hacia regiones de disponibilidad de recursos. El término de advección resultante modela el movimiento masivo de poblaciones. Con esto se obtiene el sistema depredador-presa que finalmente consideraremos, dado por el sistema de ecuaciones:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _{t}\varphi _{1}-k_{1}(\partial _{1}^{2}\varphi _{1}+\partial _{2}^{2}\varphi _{2})+a_{11}\partial _{1}\varphi _{1}+a_{12}\partial _{2}\varphi _{1}-s_{11}(1-\varphi _{1})\varphi _{1}+s_{12}\left(\frac{\varphi _{1}}{1+\alpha _{1}\varphi _{1}}\right)\varphi _{2}=0,</math>
|-
| style="text-align: center;" | <math>  \partial _{t}\varphi _{2}-k_{2}(\partial _{1}^{2}\varphi _{2}+\partial _{2}^{2}\varphi _{2})+a_{21}\partial _{1}\varphi _{2}+a_{22}\partial _{2}\varphi _{2}-s_{21}\left(\frac{\varphi _{2}}{1+\alpha _{2}\varphi _{1}}\right)\varphi _{1}+s_{22}\varphi _{2}=0, </math>
|}
|}

al cual hay que añadir las condiciones de contorno y condiciones iniciales. En estas ecuaciones, <math display="inline">s_{11}=C,</math> <math display="inline">K=1,</math> <math display="inline">s_{12}=\frac{B}{H},</math> <math display="inline">\alpha _{1}=\frac{1}{H},</math> <math display="inline">s_{21}=\frac{\kappa B}{H}</math>, <math display="inline">s_{22}=M.</math> El campo de velocidad para la presa es <math display="inline">\boldsymbol{a}_{1}=\left[a_{11},a_{12}\right]^{T}</math> y para el depredador <math display="inline">\boldsymbol{a}_{2}=\left[a_{21},a_{22}\right]^{T}</math>.

Escribiendo el sistema depredador-presa en el formato de la ecuación vectorial de CDR transitoria ([[#eq-7|7]]) tenemos:

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\begin{array}{c}\left[\begin{array}{c}\partial _{t}\varphi _{1}\\ \partial _{t}\varphi _{2} \end{array}\right]-\left[\begin{array}{cc}k_{1} & 0\\ 0 & k_{2} \end{array}\right]\partial _{1}^{2}\left[\begin{array}{c}\varphi _{1}\\ \varphi _{2} \end{array}\right]-\left[\begin{array}{cc}k_{1} & 0\\ 0 & k_{2} \end{array}\right]\partial _{2}^{2}\left[\begin{array}{c}\varphi _{1}\\ \varphi _{2} \end{array}\right]+\left[\begin{array}{cc}a_{11} & 0\\ 0 & a_{21} \end{array}\right]\partial _{1}\left[\begin{array}{c}\varphi _{1}\\ \varphi _{2} \end{array}\right]\end{array}\\ \\ +\left[\begin{array}{cc}a_{12} & 0\\ 0 & a_{22} \end{array}\right]\partial _{2}\left[\begin{array}{c}\varphi _{1}\\ \varphi _{2} \end{array}\right]+\left[\begin{array}{cc}-s_{11}\left(1-\varphi _{1}\right)& s\left(\frac{\varphi _{1}}{1+\alpha _{1}\varphi _{1}}\right)\\ -s_{21}\left(\frac{\varphi _{2}}{1+\alpha _{2}\varphi _{1}}\right)& s_{22} \end{array}\right]\left[\begin{array}{c}\varphi _{1}\\ \varphi _{2} \end{array}\right]=\left[\begin{array}{c}0\\ 0 \end{array}\right]. \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

La matriz de los parámetros de estabilización utilizada para este problema, que está descrita en <span id='citeF-33'></span>[[#cite-33|[33]]], es <math display="inline">\boldsymbol{\tau }=\mathrm{diag}(\tau _{1},\tau _{2})</math>, con:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau _{1}=\left[\frac{4k_{1}}{\left(\frac{h}{d^{2}}\right)^{2}}+\frac{2|\boldsymbol{a}_{1}|}{\frac{h}{d}}+s_{11}\left(1-\max \left(\varphi _{1}\right)\right)+s_{12}\left(1-\max \left(\varphi _{2}\right)\right)\right]^{-1},</math>
|-
| style="text-align: center;" | <math> \tau _{2}=\left[\frac{4k_{2}}{\left(\frac{h}{d^{2}}\right)^{2}}+\frac{2|\boldsymbol{a}_{2}|}{\frac{h}{d}}+s_{21}\left(1-\max \left(\varphi _{1}\right)\right)+s_{22}\left(1-\max \left(\varphi _{2}\right)\right)\right]^{-1}, </math>
|}
|}

donde recordemos que <math display="inline">h</math> corresponde al diámetro del elemento y <math display="inline">d</math> es el orden polinomial de interpolación. Las variables <math display="inline">\varphi _{1}</math> y <math display="inline">\varphi _{2}</math> las tomamos calculadas en la iteración anterior del paso de tiempo actual.

====4.4.2 Pruebas numéricas y resultados====

Consideremos la ecuación ([[#eq-22|22]]) del modelo depredador-presa dentro de un cuadrado unitario <math display="inline">\Omega =\left[0,1\right]\times \left[0,1\right]</math>, con condiciones de contorno de Dirichlet homogéneas y condiciones iniciales de las densidades de población de la presa y del depredador dadas por las distribuciones normales ([[#eq-23|23]]) y ([[#eq-24|24]]), respectivamente, escritas a continuación y mostradas en la Figura [[#img-10|10]]:

<span id="eq-23"></span>
<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\varphi _{1}(x,y,t=0)=\exp \left[{-50\left(\left(x-0.25\right)^{2}+\left(y-0.25\right)^{2}\right)}\right],</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|-
| style="text-align: center;" | <math> \varphi _{2}(x,y,t=0)=\exp \left[{-50\left(\left(x-0.75\right)^{2}+\left(y-0.75\right)^{2}\right)}\right]. </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}
|}

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Codina_913592092-test-Presa_inicial.png|424px|]]
|[[Image:Draft_Codina_913592092-test-Depredador_inicial.png|424px|Condición inicial de la densidad de población de la presa (izquierda) y del depredador (derecha).]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 10:''' Condición inicial de la densidad de población de la presa (izquierda) y del depredador (derecha).
|}

En la Tabla [[#table-3|3]] se encuentran los valores de los coeficientes de reacción <math display="inline">s_{11},</math> <math display="inline">s_{12},</math> <math display="inline">s_{21}</math> y <math display="inline">s_{22}</math> para varios casos de prueba que hemos considerado. Para todos estos casos de prueba los coeficientes de difusión se tomaron como <math display="inline">k_{11}=k_{22}=10^{-4}</math>&nbsp;m<math display="inline">^2</math>/s. El campo de velocidades se mantiene constante tanto para el depredador como para la presa. Las componentes de velocidad de la presa son <math display="inline">a_{11}=0.5</math>&nbsp;m/s, <math display="inline">a_{12}=0.5</math>&nbsp;m/s y para el depredador son <math display="inline">a_{21}=-0.5</math>&nbsp;m/s, <math display="inline">a_{22}=-0.5</math>&nbsp;m/s, con lo cual las poblaciones del depredador y de la presa son conducidas en direcciones opuestas para encontrarse frente a frente una con la otra. No consideramos ningún término fuente  y las constantes <math display="inline">\alpha _{1}</math> y <math display="inline">\alpha _{2}</math> se tomaron a la unidad. Con esto el sistema de ecuaciones del modelo depredador-presa se escribe:

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\begin{array}{c}\left[\begin{array}{c}\partial _{t}\varphi _{1}\\ \partial _{t}\varphi _{2} \end{array}\right]-\left[\begin{array}{cc}10^{-4} & 0\\ 0 & 10^{-4} \end{array}\right]\partial _{1}^{2}\left[\begin{array}{c}\varphi _{1}\\ \varphi _{2} \end{array}\right]-\left[\begin{array}{cc}10^{-4} & 0\\ 0 & 10^{-4} \end{array}\right]\partial _{2}^{2}\left[\begin{array}{c}\varphi _{1}\\ \varphi _{2} \end{array}\right]+\left[\begin{array}{cc}0.5 & 0\\ 0 & -0.5 \end{array}\right]\partial _{1}\left[\begin{array}{c}\varphi _{1}\\ \varphi _{2} \end{array}\right]\end{array}\\ \\ +\left[\begin{array}{cc}0.5 & 0\\ 0 & -0.5 \end{array}\right]\partial _{2}\left[\begin{array}{c}\varphi _{1}\\ \varphi _{2} \end{array}\right]+\left[\begin{array}{cc}-s_{11}\left(1-\varphi _{1}\right)& \quad s_{12}\left(\frac{\varphi _{1}}{1+\varphi _{1}}\right)\\ -s_{21}\left(\frac{\varphi _{2}}{1+\varphi _{1}}\right)& s_{22} \end{array}\right]\left[\begin{array}{c}\varphi _{1}\\ \varphi _{2} \end{array}\right]=\left[\begin{array}{c}0\\ 0 \end{array}\right]. \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-3'></span>Tabla. 3 Casos de prueba del modelo depredador-presa para diferentes coeficientes de reacción
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   
| style="border-left: 2px solid;border-right: 2px solid;" | <math>s_{11}</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>s_{12}</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>s_{21}</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>s_{22}</math>
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   Caso 1 
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
| style="border-left: 2px solid;border-right: 2px solid;" | 0
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  Caso 2 
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
| style="border-left: 2px solid;border-right: 2px solid;" | 1
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  Caso 3 
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
| style="border-left: 2px solid;border-right: 2px solid;" | 3 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.1
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  Caso 4 
| style="border-left: 2px solid;border-right: 2px solid;" | 1 
| style="border-left: 2px solid;border-right: 2px solid;" | 0 
| style="border-left: 2px solid;border-right: 2px solid;" | 3 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.1
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  Caso 5 
| style="border-left: 2px solid;border-right: 2px solid;" | 1 
| style="border-left: 2px solid;border-right: 2px solid;" | 2 
| style="border-left: 2px solid;border-right: 2px solid;" | 3 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.1

|}

La malla de elementos finitos usada es regular y consta de <math display="inline">50 \times 50</math>  elementos cuadrados de 4 nodos, es decir <math display="inline">\delta x=\delta y=0.02</math>, con un total de 2500 elementos y 2601 nodos. El intervalo de tiempo es <math display="inline">\left[0,1\right]</math> y el tamaño del paso de tiempo tomado es uniforme, con <math display="inline">\delta t=0.2</math>. Hemos usado el integrador temporal BDF2 y el método variacional multiescala ASGS descrito anteriormente. Presentamos resultados numéricos en diferentes instantes de tiempos de simulación, <math display="inline">t=0.2,</math> <math display="inline">0.4</math>, <math display="inline">0.6,</math> <math display="inline">0.8</math>, <math display="inline">1.0</math>, para todos los casos de prueba.

En las Figuras [[#img-11|11]], [[#img-12|12]], [[#img-13|13]], [[#img-14|14]] y [[#img-15|15]] se encuentran los resultados de los casos 1, 2, 3, 4 y 5, respectivamente, considerados en la Tabla [[#table-3|3]]. Los resultados numéricos en todos los casos son los esperados y corresponden con la realidad física de cada caso, observando además que son coincidentes con los resultados mostrados en <span id='citeF-33'></span>[[#cite-33|[33]]].

Las diferencias de los resultados de las densidades de población del depredador y de la presa con elementos de alto orden <math display="inline">Q_{2}</math> (con 10201 nodos), <math display="inline">Q_{3}</math> (con 22801 nodos) y <math display="inline">Q_{4}</math> (con 40401 nodos) están en el orden de las milésimas con respecto a los resultados con elementos lineales para <math display="inline">t=1.0</math>, tal como observamos en las Tablas [[#table-4|4]] y [[#table-5|5]], por lo que dependiendo de la escala del problema es suficiente la precisión con elementos lineales.

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|-
| <math>t</math>=0.2 s
| [[Image:Draft_Codina_913592092-test-Presa_0_2.png|300px|"Caso-1/Presa_0_2".png]]
| [[Image:Draft_Codina_913592092-test-Depredador_0_2.png|300px|"Caso-1/Depredador_0_2".png]]
|-
| <math>t</math>=0.4 s
| [[Image:Draft_Codina_913592092-test-Presa_0_4.png|300px|"Caso-1/Presa_0_4".png]]
| [[Image:Draft_Codina_913592092-test-Depredador_0_4.png|300px|"Caso-1/Depredador_0_4".png]]
|-
| <math>t</math>=0.6 s
| [[Image:Draft_Codina_913592092-test-Presa_0_6.png|300px|"Caso-1/Presa_0_6".png]]
| [[Image:Draft_Codina_913592092-test-Depredador_0_6.png|300px|"Caso-1/Depredador_0_6".png]]
|-
| <math>t</math>=0.8 s
| [[Image:Draft_Codina_913592092-test-Presa_0_8.png|300px|"Caso-1/Presa_0_8".png]]
| [[Image:Draft_Codina_913592092-test-Depredador_0_8.png|300px|"Caso-1/Depredador_0_8".png]]
|-
| <math>t</math>=1.0 s
| [[Image:Draft_Codina_913592092-test-Presa_1_0.png|300px|"Caso-1/Presa_1_0".png]]
| [[Image:Draft_Codina_913592092-test-Depredador_1_0.png|300px|"Caso-1/Depredador_1_0".png]]

|}

|-

|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Caso 1. Densidad de población de la presa (izquierda). Densidad de población del depredador (derecha). 
|}

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|-
| <math>t</math>=0.2 s
| [[Image:Draft_Codina_913592092-test-Presa_0_2-dup1.png|300px|"Caso-2/Presa_0_2".png]]
| [[Image:Draft_Codina_913592092-test-Depredador_0_2-dup1.png|300px|"Caso-2/Depredador_0_2".png]]
|-
| <math>t</math>=0.4 s
| [[Image:Draft_Codina_913592092-test-Presa_0_4-dup1.png|300px|"Caso-2/Presa_0_4".png]]
| [[Image:Draft_Codina_913592092-test-Depredador_0_4-dup1.png|300px|"Caso-2/Depredador_0_4".png]]
|-
| <math>t</math>=0.6 s
| [[Image:Draft_Codina_913592092-test-Presa_0_6-dup1.png|300px|"Caso-2/Presa_0_6".png]]
| [[Image:Draft_Codina_913592092-test-Depredador_0_6-dup1.png|300px|"Caso-2/Depredador_0_6".png]]
|-
| <math>t</math>=0.8 s
| [[Image:Draft_Codina_913592092-test-Presa_0_8-dup1.png|300px|"Caso-2/Presa_0_8".png]]
| [[Image:Draft_Codina_913592092-test-Depredador_0_8-dup1.png|300px|"Caso-2/Depredador_0_8".png]]
|-
| <math>t</math>=1.0 s
| [[Image:Draft_Codina_913592092-test-Presa_1_0-dup1.png|300px|"Caso-2/Presa_1_0".png]]
| [[Image:Draft_Codina_913592092-test-Depredador_1_0-dup1.png|300px|"Caso-2/Depredador_1_0".png]]

|}

|-

|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Caso 2. Densidad de población de la presa (izquierda). Densidad de población del depredador (derecha). 
|}

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|-
| <math>t</math>=0.2 s
| [[Image:Draft_Codina_913592092-test-Caso3_presa_0_2.png|300px|"Caso-3/Caso3_presa_0_2".png]]
| [[Image:Draft_Codina_913592092-test-Caso3_depredador_0_2.png|300px|"Caso-3/Caso3_depredador_0_2".png]]
|-
| <math>t</math>=0.4 s
| [[Image:Draft_Codina_913592092-test-Caso3_presa_0_4.png|300px|"Caso-3/Caso3_presa_0_4".png]]
| [[Image:Draft_Codina_913592092-test-Caso3_depredador_0_4.png|300px|"Caso-3/Caso3_depredador_0_4".png]]
|-
| <math>t</math>=0.6 s
| [[Image:Draft_Codina_913592092-test-Caso3_presa_0_6.png|300px|"Caso-3/Caso3_presa_0_6".png]]
| [[Image:Draft_Codina_913592092-test-Caso3_depredador_0_6.png|300px|"Caso-3/Caso3_depredador_0_6".png]]
|-
| <math>t</math>=0.8 s
| [[Image:Draft_Codina_913592092-test-Caso3_presa_0_8.png|300px|"Caso-3/Caso3_presa_0_8".png]]
| [[Image:Draft_Codina_913592092-test-Caso3_depredador_0_8.png|300px|"Caso-3/Caso3_depredador_0_8".png]]
|-
| <math>t</math>=1.0 s
| [[Image:Draft_Codina_913592092-test-Caso3_presa_1_0.png|300px|"Caso-3/Caso3_presa_1_0".png]]
| [[Image:Draft_Codina_913592092-test-Caso3_depredador_1_0.png|300px|"Caso-3/Caso3_depredador_1_0".png]]

|}

|-

|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Caso 3. Densidad de población de la presa (izquierda). Densidad de población del depredador (derecha). 
|}

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|-
| <math>t</math>=0.2 s
| [[Image:Draft_Codina_913592092-test-Caso4_presa_0_2.png|300px|"Caso-4/Caso4_presa_0_2".png]]
| [[Image:Draft_Codina_913592092-test-Caso4_depredador_0_2.png|300px|"Caso-4/Caso4_depredador_0_2".png]]
|-
| <math>t</math>=0.4 s
| [[Image:Draft_Codina_913592092-test-Caso4_presa_0_4.png|300px|"Caso-4/Caso4_presa_0_4".png]]
| [[Image:Draft_Codina_913592092-test-Caso4_depredador_0_4.png|300px|"Caso-4/Caso4_depredador_0_4".png]]
|-
| <math>t</math>=0.6 s
| [[Image:Draft_Codina_913592092-test-Caso4_presa_0_6.png|300px|"Caso-4/Caso4_presa_0_6".png]]
| [[Image:Draft_Codina_913592092-test-Caso4_depredador_0_6.png|300px|"Caso-4/Caso4_depredador_0_6".png]]
|-
| <math>t</math>=0.8 s
| [[Image:Draft_Codina_913592092-test-Caso4_presa_0_8.png|300px|"Caso-4/Caso4_presa_0_8".png]]
| [[Image:Draft_Codina_913592092-test-Caso4_depredador_0_8.png|300px|"Caso-4/Caso4_depredador_0_8".png]]
|-
| <math>t</math>=1.0 s
| [[Image:Draft_Codina_913592092-test-Caso4_presa_1_0.png|300px|"Caso-4/Caso4_presa_1_0".png]]
| [[Image:Draft_Codina_913592092-test-Caso4_depredador_1_0.png|300px|"Caso-4/Caso4_depredador_1_0".png]]

|}

|-

|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' Caso 4. Densidad de población de la presa (izquierda). Densidad de población del depredador (derecha). 
|}

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|
{|  style="text-align: center; margin: 1em auto;min-width:50%;width:100%;"
|-
| <math>t</math>=0.2 s
| [[Image:Draft_Codina_913592092-test-Caso5_presa_0_2.png|300px|"Caso-5/Caso5_presa_0_2".png]]
| [[Image:Draft_Codina_913592092-test-Caso5_depredador_0_2.png|300px|"Caso-5/Caso5_depredador_0_2".png]]
|-
| <math>t</math>=0.4 s
| [[Image:Draft_Codina_913592092-test-Caso5_presa_0_4.png|300px|"Caso-5/Caso5_presa_0_4".png]]
| [[Image:Draft_Codina_913592092-test-Caso5_depredador_0_4.png|300px|"Caso-5/Caso5_depredador_0_4".png]]
|-
| <math>t</math>=0.6 s
| [[Image:Draft_Codina_913592092-test-Caso5_presa_0_6.png|300px|"Caso-5/Caso5_presa_0_6".png]]
| [[Image:Draft_Codina_913592092-test-Caso5_depredador_0_6.png|300px|"Caso-5/Caso5_depredador_0_6".png]]
|-
| <math>t</math>=0.8 s
| [[Image:Draft_Codina_913592092-test-Caso5_presa_0_8.png|300px|"Caso-5/Caso5_presa_0_8".png]]
| [[Image:Draft_Codina_913592092-test-Caso5_depredador_0_8.png|300px|"Caso-5/Caso5_depredador_0_8".png]]
|-
| <math>t</math>=1.0 s
| [[Image:Draft_Codina_913592092-test-Caso5_presa_1_0.png|300px|"Caso-5/Caso5_presa_1_0".png]]
| [[Image:Draft_Codina_913592092-test-Caso5_depredador_1_0.png|300px|"Caso-5/Caso5_depredador_1_0".png]]

|}

|-

|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Caso 5. Densidad de población de la presa (izquierda). Densidad de población del depredador (derecha). 
|}


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-4'></span>Tabla. 4 Densidad de población de la presa (<math>\varphi _{1}</math>)
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   
| style="border-left: 2px solid;border-right: 2px solid;" | Caso 1 
| style="border-left: 2px solid;border-right: 2px solid;" | Caso 2 
| style="border-left: 2px solid;border-right: 2px solid;" | Caso 3 
| style="border-left: 2px solid;border-right: 2px solid;" | Caso 4 
| style="border-left: 2px solid;border-right: 2px solid;" | Caso 5
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">Q_{1}</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48614 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48614 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48614 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.69991 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.56545
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">Q_{2}</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48847 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48847 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48847 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.70560 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.56853
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">Q_{3}</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48847 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48847 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48847 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.70709 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.56928
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">Q_{4}</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48848 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48848 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48848 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.70786 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.56982

|}


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-5'></span>Tabla. 5 Densidad de población del depredador (<math>\varphi _{2}</math>)
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   
| style="border-left: 2px solid;border-right: 2px solid;" | Caso 1 
| style="border-left: 2px solid;border-right: 2px solid;" | Caso 2 
| style="border-left: 2px solid;border-right: 2px solid;" | Caso 3 
| style="border-left: 2px solid;border-right: 2px solid;" | Caso 4 
| style="border-left: 2px solid;border-right: 2px solid;" | Caso 5
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">Q_{1}</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48614 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.22615 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.63921 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.70337 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.65541
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">Q_{2}</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48847 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.22685 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.63872 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.71223 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.66228
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">Q_{3}</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48847 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.22684 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.63954 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.71348 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.66352
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | <math display="inline">Q_{4}</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.48848 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.22684 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.64002 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.71412 
| style="border-left: 2px solid;border-right: 2px solid;" | 0.66422

|}

==5 Conclusiones==

En este artículo se ha presentado la aproximación del modelo acoplado de las ecuaciones del movimiento de un fluido en aguas poco profundas junto con las ecuaciones de convección-difusión-reacción del transporte de contaminantes mediante formulaciones estabilizadas de elementos finitos de alto orden (hasta el cuarto orden).

Las formulaciones estabilizadas de elementos finitos para la ecuación vectorial de CDR transitoria, ASGS y OSS, nos permiten utilizar igual interpolación en todas las variables, que para el modelo del transporte de contaminantes en aguas poco profundas son las dos componentes de velocidad <math display="inline">U_{1}</math>, <math display="inline">U_{2}</math>, la elevación <math display="inline">\eta </math> de la superficie libre del agua y las concentraciones de contaminantes. Asimismo, las formulaciones ASGS y OSS para la ecuación vectorial de CDR transitoria nos permiten tratar con flujos de convección y reacción dominantes, incluyendo problemas no lineales en los términos convectivos y reactivos. Como ejemplo de problema de reacción no lineal hemos presentado el modelo depredador-presa, con dos concentraciones acopladas a través del término reactivo como incógnitas.

Para utilizar elementos finitos de alto orden, en la matriz de los parámetros de estabilización se ha incluido la dependencia con el orden polinomial de interpolación.

Los resultados de las gráficas de convergencia muestran claramente que las formulaciones para los métodos estabilizados ASGS y OSS son capaces de aproximar correctamente el problema vectorial transitorio de CDR con convección dominante, particularizado al sistema acoplado de cuatro ecuaciones correspondientes al flujo en aguas poco profundas y transporte de un contaminante.

Los resultados calculados del transporte de un contaminante en una cavidad cuadrada no muestran errores significativos de dispersión o disipación, por lo que concuerdan muy bien con los valores teóricos, y son mejores que los resultados presentados en la literatura, como en <span id='citeF-3'></span><span id='citeF-4'></span><span id='citeF-5'></span>[[#cite-3|[3,4,5]]].

Asimismo, se observa que los resultados obtenidos del transporte de un contaminante en el golfo de Roses y en la desembocadura del río Guadalquivir  son completamente satisfactorios, tanto para el campo de velocidades y el nivel de la superficie libre del agua como para la distribución de la concentración de un contaminante.

Finalmente los resultados del modelo depredador-presa son coherentes con la realidad física de cada caso, y muestran la robustez de las formulaciones estabilizadas ASGS y OSS para la ecuación vectorial de CDR transitoria en problemas de reacción dominante y de no linealidad en el término de reacción.

==Agradecimientos==

R. Codina agradece la ayuda recibida a través del programa ICREA Academia, del Gobierno de Catalunya.

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  Z. Zlatev. The Air Pollution Problem. In: Computer Treatment of Large Air Pollution Models. Environmental Science and Technology Library, vol 2. Springer, Dordrecht, 1995.

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]'''  O.C. Zienkiewicz and R.L. Taylor. The Finite Element Method, volume 3, Fluid Mechanics. Butterworth- Heinemann, 5 edition, 2000.

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]'''  S. Li and C. J. Duffy. Fully-coupled modeling of shallow water flow and pollutant transport on unstructured grids. Procedia Environmental Sciences, 13:20982121, 2012.

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''  F. Behzadi and J. C. Newman III. A Semi-Discrete SUPG Method for Contaminant Transport in Shallow Water Models. Procedia Computer Science, 80:1313-1323, 2016.

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  F. Benkhaldoun, I. Elmahi, M. Seaïd. Well-balanced finite volume schemes for pollutant transport by shallow water equations on unstructured meshes. J. Comput. Phys. 1:180-203, 2007.

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  T. Komatsu, A. K. Ohgushi, and K. Asai. Refined numerical scheme for advective transport in diffusion simulation. J. Hydraul. Eng., 123:41-50, 1997.

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]'''  L. Begnudelli, B. F. Sanders, Unstructured grid finite-volume algorithm for shallow-water flow and scalar transport with wetting and drying. Journal of Hydraulic Engineering, 132(4):371-384, 2006.

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]'''  L. Cea, M. E. Vázquez-Cedón, Unstructured finite volume discretisation of bed friction and convective flux in solute transport models linked to the sallow water equations. Journal of Computational Physics. 231(8):3317-3330, 2012.

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''  V. Caleffi, A. Valiani. A 2D local discontinuous Galerkin method for contaminant transport in channel bends. Computers & Fluids. 88:629-642, 2013.

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]'''  S. Pavan, J. M. Hervouet, M. Ricchiuto, and R. Ata. A second order residual based predictorcorrector approach for time dependent pollutant transport. Journal of Computational Physics. 318:122141, 2016.

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''  L. Cai, W.-X. Xie, J.-H. Feng, and J. Zhou.  Computations of transport of pollutant in shallow water. Applied Mathematical Modelling. 31:490498, 2007.

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]'''  L. Postma, J.-M. Hervouet, Compatibility between finite volumes and finite elements using solutions of shallow water equations for substance transport, Int. J. Numer. Methods Fluids. 53(9):14951507, 2007.

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]'''  S. Pavan, R. Ata, J.-M. Hervouet, Finite volume schemes and residual distribution schemes for pollutant transport on unstructured grids, Environ. Earth Sci. 74:73377356, 2015.

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  H. Liu, J. G. Zhou, M. Li, Y. Zhao, Multi-block Lattice Boltzmann Simulations of Solute Transport in Shallow Water Flows, Advances in Water Resources (2013).

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''  H. Hammou, I. Ginzburg and M. Boulerhcha, Two-relaxation-times lattice Boltzmann schemes for solute transport in unsaturated water flow with a focus on stability. Adv Wat Res. 4(6):779-793, 2011.

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]'''  J. G. Zhou, A lattice Boltzmann method for solute transport. Int J Numer Methods Fluids. 61:848-63, 2009.

<div id="cite-17"></div>
'''[[#citeF-17|[17]]]'''   R. Codina, S. Badia, J. Baiges and J. Principe, Variational Multiscale Methods in Computational Fluid Dynamics,  In:  ''Encyclopedia of Computational Mechanics'' (eds E. Stein, R. Borst and T.J.R. Hughes), doi: 10.1002/9781119176817.ecm21172017, 2017.

<div id="cite-18"></div>
'''[[#citeF-18|[18]]]'''  B. Rogers, M. Fujihara and A. G. L. Borthwick, Adaptive Q-tree Godunov-type scheme for shallow water equations. Int J Numer Methods Fluids. 35:247-80, 2001.

<div id="cite-19"></div>
'''[[#citeF-19|[19]]]'''  H. Yamamoto. Lecture 8: The Shallow-Water Equations, 

https://www.coursehero.com/file/35480149/othersonline-HW4pdf/

<div id="cite-20"></div>
'''[[#citeF-20|[20]]]'''  T. J. R. Hughes. Multiscale phenomena: Green's functions, the Drichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Computer Methods in Applied Mechanics and Engineering'', ''Vol. 127, No. 1-4, pp. 387-401, 1995.

<div id="cite-21"></div>
'''[[#citeF-21|[21]]]'''  T. J. R Hughes, G. O. Feijóo, L. Mazzei, and J. B. Quincy. The variational multiscale method-a paradign for computational mechanics. Computer Methods in Applied Mechanics and Egineering, 166:3-24, 1998.

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]'''  R. Codina, On stabilized finite element methods for linear systems of convection-diffusion-reaction equations, Computer Methods in Applied Mechanics and Engineering 188:61-82, 2000.

<div id="cite-23"></div>
'''[[#citeF-23|[23]]]'''  R. Codina, A stabilized finite element method for generalized stationary incompressible flows, Computer Methods in Applied Mechanics and Enginineering 190:2681-2706, 2001.

<div id="cite-24"></div>
'''[[#citeF-24|[24]]]'''  R. Codina, Comparison of some finite element methods for solving the diffusion-convection-reaction equation, Computer Methods in Applied Mechanics and Engineering 156:185-210, 1998.

<div id="cite-25"></div>
'''[[#citeF-25|[25]]]'''  R. Codina and J. Blasco, Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales, Computing and Visualization in Science 4 (2002) 167-174.

<div id="cite-26"></div>
'''[[#citeF-26|[26]]]'''  R. Codina, Stabilized finite element approximation of transient incompressible flows using orthogonal subscales, Computer Methods in Applied Mechanics and Engineering, 191:4295-4321, 2002.

<div id="cite-27"></div>
'''[[#citeF-27|[27]]]'''  A. Villota and R. Codina, Approximation of the shallow water equations with higher order finite elements and variational multiscale methods, Rev. int. métodos numér. cálc. diseño ing. (2018). Vol. 34, 1, 28 URL https://www.scipedia.com/public/Villota_Codina_a.

<div id="cite-28"></div>
'''[[#citeF-28|[28]]]'''  R. Codina, Finite Element Approximation of the Convection-Diffusion Equation: Subgrid-Scale Spaces, Local Instabilities and Anisotropic Space-Time Discretizations, in Bail 2010 - Boundary and Interior Layers, Computational and Asymptotic Methods, Lecture Notes in Computational Science and Engineering, Eds. C. Clavero, J.L. Gracia and F.J. Lisbona, volume 81, pages 85-97. Springer, 2011.

<div id="cite-29"></div>
'''[[#citeF-29|[29]]]'''  R. Codina, J. M. González-Ondina, G. Díaz-Hernández and J. Principe, Finite element approximation of the modified Boussinesq equations using a stabilized formulation, International Journal for Numerical Methods in Fluids, 57:1249-1268, 2008.

<div id="cite-30"></div>
'''[[#citeF-30|[30]]]'''   A. Villota and R. Codina, Approximation of the scalar convection-diffusion-reaction equation with stabilized finite element formulations of high order, Rev. int. métodos numér. cálc. diseño ing. (2019). Vol. 35, (1), 6 URL https://www.scipedia.com/public/Villota_Codina_a

<div id="cite-31"></div>
'''[[#citeF-31|[31]]]'''  J. Principe and R. Codina. On the stabilization parameter in the subgrid scale approximation of scalar convection-diffusion-reaction equations on distorted meshes. Computer Methods in Applied Mechanics and Engineering, 199:1386-1402, 2010.

<div id="cite-32"></div>
'''[[#citeF-32|[32]]]'''  R. Codina,  On hp convergence of stabilized finite element approximations of the convection-diffusion equation, SeMA Journal, 75:591-606, 2018.

<div id="cite-33"></div>
'''[[#citeF-33|[33]]]'''  S. K. Sheshachala and R. Codina, Finite element modeling of nonlinear reaction-diffusion-advection systems of equations, International Journal of Numerical Methods for Heat & Fluid Flow, https://doi.org/10.1108/HFF-02-2018-0077, 2018.