Numerical modeling of the thermal discharge of the Laguna Verde power station

==Resumen==

El presente trabajo contribuye al conocimiento del comportamiento de las descargas térmicas en aguas marinas costeras mediante la implementación y la aplicación de un modelo numérico que resuelve las ecuaciones de Navier-Stokes-Reynolds para aguas someras y la ecuación de la energía para la temperatura. El modelo numérico toma en cuenta el flujo de calor en la capa superficial, donde interactúan la superficie libre del mar y la atmósfera. Como caso de estudio, se analiza la dispersión de la pluma térmica de la Central Nucleoeléctrica Laguna Verde (CNLV), ubicada en el Estado de Veracruz, México. Las simulaciones numéricas se llevan a cabo teniendo en cuenta la información ambiental, que incluye información batimétrica, oceanográfica, meteorológica e hidrológica, y las condiciones de operación de la descarga de la CNLV. Con base en comparaciones con datos medidos en campo y en el criterio de eficiencia de Nash-Suffle para verificar la calidad de la solución numérica, se considera que los resultados obtenidos son bastante satisfactorios.

==Abstract==

This work contributes to the study of nuclear plant thermal discharges in coastal areas by using a numerical model which solves the Navier-Stokes-Reynolds equations for shallow waters and the energy equation for computing temperature variations. The numerical model takes into account the heat flux given in the upper layer, where the free surface and the atmosphere interact. In this study, the thermal plume dispersion from the nuclear power plant Laguna Verde, Veracruz, Mexico, is analyzed. Bathymetry, oceanographic, meteorological, hydrologic and plant operating data are used to run numerical simulations. The results are compared against observed data showing good agreement. The Nash-Suffles criterion is also applied to verify the quality of the numerical solution obtaining suitable results.

==Palabras clave==

Central nucleoeléctrica ; Dinámica costera ; Recirculación de plumas térmicas ; Modelación numérica

==Keywords==

Thermal plume dispersion ; Coastal hydrodynamics ; Nuclear power station ; Multilayer numerical modeling

==Nomenclatura==

''C''<sub>''smag''</sub>- Constante de Smagorinsky (adim)

<math display="inline">C_{\rho w}</math>- Calor específico del agua (m<sup>2</sup> /°C s<sup>2</sup> )                          

''f''- Fuerza de Coriolis por unidad de masa (1/s)

''g''- Aceleración de la gravedad (m/s<sup>2</sup> )                          

''H''- Profundidad total (m)

''k''- Constante de von Kármán (adim)

''k''<sub>''p''</sub>- Constante de ajuste en la ecuación de estado (°C)

''K''<sub>''TH''</sub>- Coeficiente de dispersión en el plano horizontal (m<sup>2</sup> C)                          

''K''<sub>''TV''</sub>- Coeficiente de dispersión en el plano vertical (m<sup>2</sup> /s)                          

ℓ- Longitud de mezclado (m)

''P''- Presión hidrostática (Pa)

''P''<sub>''atm''</sub>- Presión atmosférica (Pa)

''P''<sub>''rT''</sub>- Numero de Prandtl turbulento (adim)

''Q''- Flujo de calor neto (kg/s<sup>3</sup> )                          

''S''- Salinidad (ups)

''t''- Tiempo (s)

''T''- Temperatura del agua (°C)

''U''- Magnitud de la componente de velocidad en la dirección ''x''  (m/s)                          

''V''- Magnitud de la componente de velocidad en la dirección ''y''  (m/s)                          

''W''- Magnitud de la componente de velocidad en la dirección ''z''  (m/s)                          

''x'' , ''y'' , ''z''- Plano cartesiano de referencia (m)

''z''<sub>''b''</sub>- Distancia vertical desde el fondo hasta un punto cualquiera (m)

''η''- Variación de la superficie libre con respecto al nivel medio (m)

''ν''<sub>''TH''</sub>- Coeficiente de viscosidad turbulenta horizontal (m<sup>2</sup> /s)                          

''ν''<sub>''TV''</sub>- Coeficiente de viscosidad turbulenta vertical (m<sup>2</sup> /s)                          

Δ''A''- Área de control (m<sup>2</sup> )                          

Δ''Q''- Aporte de temperatura por flujo de calor en la superficie (°C/s)

Δ''t''- Incremento del paso de tiempo (s)

Δ''V''- Volumen de control (m<sup>3</sup> )                          

Δ''x''- Ancho de celda en la dirección ''x''  (m)                          

Δ''y''- Ancho de celda en la dirección ''y''  (m)                          

''ρ''- Densidad del agua (kg/m<sup>3</sup> )                          

''ρ''<sub>0</sub>- Densidad de referencia (kg/m<sup>3</sup> )                  

==1. Introducción==

En los últimos años, los cuerpos de agua costeros han sido objeto de contaminación ambiental por obras que ha construido el hombre. Aunque las leyes son cada vez más estrictas, la contaminación de los mares en la zona costera es cada vez más significativa, lo que ocasiona una alteración de los parámetros físicos, químicos y biológicos del agua que afecta directamente al ecosistema [[#bib0005|[1]]]  and [[#bib0010|[2]]] . Algunas de las principales fuentes de contaminación en la zona costera son las que se derivan de las actividades industriales, agrícolas y domésticas. De manera particular, en lo que se refiere a las actividades de tipo industrial, cabe mencionar las descargas térmicas de plantas de generación de energía que, al verter aguas residuales al mar a temperaturas mayores a las del cuerpo de agua receptor, provocan una contaminación térmica del agua. Desde una perspectiva ambiental, esta variación local de la temperatura del agua puede perturbar considerablemente el ecosistema acuático [[#bib0015|[3]]] , causando alteraciones en los organismos que ahí se desarrollan. Se ha demostrado que un ligero cambio en la temperatura del agua (1 o 2 °C) puede tener un impacto considerable en el ambiente [[#bib0020|[4]]]  and [[#bib0025|[5]]] . Por otro lado, desde un punto de vista técnico-operacional, en ocasiones se presenta en las plantas de generación de energía el problema de la recirculación de agua caliente por la mala ubicación de la obra de toma o de una o varias descargas, lo que se traduce en una disminución de la capacidad generadora de energía de la planta. De aquí la necesidad de estudiar y monitorizar la dispersión de la pluma térmica en el mar y determinar así su orientación y su área de influencia.      

El seguimiento de la pluma térmica de plantas de generación de energía se ha llevado a cabo mediante el uso de diferentes técnicas para medir la temperatura del mar, entre las que destacan las mediciones en sitio con boyas oceanográficas, las mediciones en cruceros oceanográficos, las imágenes de satélites radiométricos y la modelación numérica. El uso de boyas oceanográficas permite obtener información a lo largo de la columna de agua, aunque se limita a mediciones puntuales. Las mediciones realizadas en cruceros oceanográficos, además de que proporcionan información de la columna de agua, permiten medir a lo largo de diferentes transectos, cubriendo una zona determinada. Sin embargo, estas mediciones suelen tomar tiempo y su costo es generalmente elevado. Otra técnica utilizada en el estudio y seguimiento de plumas térmicas es la medición de la temperatura superficial con satélites radiométricos, que proporcionan una visión sinóptica de la temperatura en grandes extensiones de agua [[#bib0030|[6]]] , [[#bib0035|[7]]]  and [[#bib0040|[8]]] . Pese a los avances significativos de la tecnología satelital en los últimos años, esta técnica se limita a la capa superficial del mar y a la resolución espacial con que cuentan los satélites actualmente en órbita, 120 m en el caso de los Landsat-4 y 5 [[#bib0045|[9]]] , por citar un ejemplo de los satélites que cuentan con mayor resolución.      

La dinámica de fluidos computacional se utiliza cada vez más en estudios de calidad del agua y de otros temas de tipo ambiental. Las ecuaciones que expresan los principios de conservación de cantidad de movimiento, de masa y de energía, sirven para definir las propiedades macroscópicas del movimiento de los fluidos, y ayudan así al estudio de la difusión y el transporte de materia a través de un fluido en movimiento. La solución de dichas ecuaciones a través del desarrollo de modelos numéricos parece ser la técnica más viable en el estudio de contaminantes en general, entre los que se incluyen desde luego las plumas térmicas. La viabilidad en el uso de los modelos numéricos obedece principalmente a 2 razones: la primera es el progreso en la modelación multidisciplinaria en las últimas décadas, en la que se agrupan la mayoría de los factores que interactúan en el medio (hidrodinámicos, meteorológicos y de calidad del agua) en un modelo operacional capaz de reproducir un determinado fenómeno físico; la segunda es el avance en la tecnología computacional, incluyendo el uso de procesos en paralelo, que ha permitido que por primera vez sea factible utilizar modelos físicos y ecológicos de manera acoplada y con la misma resolución [[#bib0050|[10]]] .      

En el presente trabajo se lleva a cabo un estudio de dispersión térmica en el mar mediante la implementación de un modelo numérico. El caso de estudio es la pluma térmica de la Central Nucleoeléctrica Laguna Verde (CNLV), localizada sobre el litoral del Golfo de México (GDM), en el Estado de Veracruz, México (coordenadas 19°43′15″ N 96°24′23″ O). La CNLV es la única planta nuclear en el país y aporta el 3% de la generación total de energía eléctrica. Desde el punto de vista ecológico, los efectos generados por la operación de la CNLV en el ambiente ya han sido estudiados y analizados por varios autores [[#bib0005|[1]]]  and [[#bib0010|[2]]] , así como los efectos producidos por el oleaje incidente debido a las obras de protección [[#bib0055|[11]]]  and [[#bib0060|[12]]] . Por otro lado, no existe un estudio de la dispersión de la pluma térmica enfocado en determinar un posible problema de recirculación de agua caliente que pudiera afectar a las condiciones de operación de la CNLV. Así pues, el objetivo de este trabajo se plantea en 2 vertientes: validar el código numérico implementado de acuerdo con los datos medidos en sitio, y determinar mediante la modelación numérica si la CNLV se ve afectada por el problema de la recirculación de agua caliente.      

En la siguiente sección se muestra el sistema de ecuaciones que resuelve el modelo numérico; en la sección 3 se detalla el esquema de solución numérica; en la sección 4 se hace un planteamiento sobre los aspectos que se consideran importantes para la validación del modelo numérico implementado; en la sección 5 se describen los escenarios de simulación para el caso de estudio de la pluma térmica de la CNLV y se muestran y discuten los resultados. El trabajo termina con las conclusiones finales.

==2. Ecuaciones gobernantes==

Las ecuaciones que se resuelven a continuación son las ecuaciones para aguas someras deducidas a partir de las ecuaciones de Navier-Stokes-Reynolds para flujos a superficie libre. Se considera la aproximación de aguas someras, donde las velocidades en el plano horizontal son mayores que la velocidad vertical, y las escalas de las longitudes horizontales son mucho mayores que la escala de longitud vertical, es decir Δ''x''  ≈ Δ''y''  ≫ Δ''z'' .      

Así, se resuelven las ecuaciones que a continuación se enumeran.

Ecuaciones para las velocidades en el plano horizontal:

<span id='eq0005'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial x}+V\frac{\partial U}{\partial y}+W\frac{\partial U}{\partial z}=-\frac{1}{{\rho }_0}\frac{\partial P}{\partial x}</math>
|-
|<math>+{\nu }_{TH}\left(\frac{{\partial }^2U}{\partial x^2}+\frac{{\partial }^2U}{\partial y^2}\right)+\frac{\partial }{\partial z}\left({\nu }_{TV}\frac{\partial U}{\partial z}\right)+f_x,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)
|}

<span id='eq0010'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial V}{\partial t}+U\frac{\partial V}{\partial x}+V\frac{\partial V}{\partial y}+W\frac{\partial V}{\partial z}=-\frac{1}{{\rho }_0}\frac{\partial P}{\partial y}</math>
|-
|<math>+{\nu }_{TH}\left(\frac{{\partial }^2V}{\partial x^2}+\frac{{\partial }^2V}{\partial y^2}\right)+\frac{\partial }{\partial z}\left({\nu }_{TV}\frac{\partial V}{\partial z}\right)-f_y.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)
|}

Ecuación de continuidad:

<span id='eq0015'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}+</math><math>\frac{\partial W}{\partial z}=0.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)
|}

Ecuación para la superficie libre:

<span id='eq0020'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial U}{\partial x}+\frac{\partial }{\partial x}\left({\int }_z^{\eta }Udz\right)+</math><math>\frac{\partial }{\partial y}\left({\int }_z^{\eta }Vdz\right)=</math><math>0.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4)
|}

Las variaciones de densidad (''ρ'' ) son función de las variaciones de presión (''P'' ), temperatura (''T'' ) y salinidad (''S'' ), y son estimadas a través de la siguiente ecuación de estado definida por la UNESCO  [[#bib0065|[13]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\rho (S,T,P)=\frac{{\rho }_0}{1-\frac{P}{k_p}}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)
|}

donde ''k''<sub>''p''</sub>  es una constante de ajuste que se determina por módulos [[#bib0065|[13]]] , y es función de las variaciones de ''P'' , ''T''  y ''S'' .      

Las ecuaciones para la temperatura y la salinidad son:

<span id='eq0030'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial T}{\partial t}+U\frac{\partial T}{\partial x}+V\frac{\partial T}{\partial y}+W\frac{\partial T}{\partial z}=</math>
|-
|<math>+K_{TH}\left(\frac{{\partial }^2T}{\partial x^2}+\frac{{\partial }^2T}{\partial y^2}\right)+\frac{\partial }{\partial z}\left(K_{TV}\frac{\partial T}{\partial z}\right)\pm Q,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)
|}

<span id='eq0035'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial S}{\partial t}+U\frac{\partial S}{\partial x}+V\frac{\partial S}{\partial y}+W\frac{\partial S}{\partial z}=</math>
|-
|<math>+K_{TH}\left(\frac{{\partial }^2S}{\partial x^2}+\frac{{\partial }^2S}{\partial y^2}\right)+\frac{\partial }{\partial z}\left(K_{TV}\frac{\partial S}{\partial z}\right).</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)
|}

En el sistema de ecuaciones que gobiernan las aguas someras se evita la ecuación del movimiento para calcular la velocidad vertical ''W''  se aproxima a partir de la ecuación  [[#eq0015|(3)]]  como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial W}{\partial z}=-\left(\frac{\partial U}{\partial x}+\right. </math><math>\left. \frac{\partial V}{\partial y}\right).</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)
|}

Por otro lado, el término de presión en la ecuación [[#eq0005|(1)]]  es sustituido, considerando la aproximación hidrostática:

<span id='eq0045'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>-\frac{1}{{\rho }_0}\frac{\partial P}{\partial x}=</math><math>-\frac{\rho g}{{\rho }_0}\frac{\partial \eta }{\partial x}-</math><math>\frac{g}{{\rho }_0}{\int }_z^{\eta }\frac{\partial {\rho }^{{'}}}{\partial z}dz-</math><math>\frac{1}{{\rho }_0}\frac{\partial P_{atm}}{\partial x}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 9)
|}

El primer término de la derecha de la ecuación [[#eq0045|(9)]]  es el término barotrópico, que representa al gradiente de presión horizontal definido por la variación de la superficie libre. El segundo término es el término baroclínico, y describe la diferencia de densidades en la columna vertical de agua. El tercer término considera la contribución de la presión atmosférica. Análogamente, se plantea la misma definición para la variación del gradiente de presión de la ecuación [[#eq0010|(2)]] .      

===2.1. Los coeficientes de viscosidad turbulenta===

El coeficiente de viscosidad turbulenta horizontal ''ν''<sub>''TH''</sub>  que aparece en las ecuaciones [[#eq0005|(1)]]  y [[#eq0010|(2)]]  se estima mediante el modelo de Smagorinsky como sigue:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\nu }_{TH}=C_{smag}\Delta x\Delta y</math>
|-
|<math>\left[{\left(\frac{\partial U}{\partial x}\right)}^2+\frac{1}{2}{\left(\frac{\partial V}{\partial x}+\frac{\partial U}{\partial y}\right)}^2+{\left(\frac{\partial V}{\partial y}\right)}^2\right],</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10)
|}

donde ''C''<sub>''smag''</sub>  es la constante de Smagorinsky, cuyos valores fluctúan entre 0,02 y 0,5. Un valor típico que es utilizado es de 0,1 [[#bib0070|[14]]] . Para el cálculo del coeficiente de viscosidad turbulenta vertical ''ν''<sub>''TV''</sub>  se utiliza el modelo siguiente [[#bib0075|[15]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\nu }_{TV}=l\sqrt{{\left(\frac{\partial U}{\partial z}\right)}^2+{\left(\frac{\partial V}{\partial z}\right)}^2},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11)
|}

donde ℓ es la longitud de mezclado y se define como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>l=\begin{array}{lll}
k\cdot z_b & si & z_b<0,20H\\
0,20H\cdot k & si & z_b>0,20H,
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12)
|}

donde ''k''  es la constante de Von Kármán igual a 0,40 y ''z''<sub>''b''</sub>  es la distancia vertical desde el fondo hasta el punto donde se calcula ''ν''<sub>''TV''</sub> .      

Los coeficientes de difusión de temperatura y salinidad de las ecuaciones [[#eq0030|(6)]]  y [[#eq0035|(7)]]  son evaluados como <math display="inline">K_{TH}=</math><math>\frac{{\nu }_{TH}}{{Pr}_T}</math>  y <math display="inline">K_{TV}=</math><math>\frac{{\nu }_{TV}}{{Pr}_T}</math> , donde ''Pr''<sub>''T''</sub>  es el número de Prandtl turbulento igual a 0,9  [[#bib0075|[15]]]  and [[#bib0080|[16]]] .      

===2.2. El intercambio de calor océano-atmósfera===

En la ecuación [[#eq0030|(6)]]  aparece el término del flujo de calor neto ''Q''  en la superficie libre del agua, que es la suma de la radiación de longitud de onda larga, de la radiación de longitud de onda corta, de la evaporación y del calor sensible. Para cuantificar el flujo de calor neto, el modelo numérico toma en cuenta datos climatológicos y meteorológicos, tales como la magnitud y la dirección del viento, la temperatura ambiente, la nubosidad y la radiación, y la temperatura superficial y del fondo del agua  [[#bib0085|[17]]] . Para el cálculo de ''Q''  se utilizó la forma linearizada de Haney, que permite estimar el flujo de calor neto en la interfase agua-aire como sigue:

<span id='eq0065'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Q=\frac{\Delta Q(\Delta A)}{\rho C_{\rho w}(\Delta V)}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)
|}

==3. Esquema de solución numérica==

El modelo numérico desarrollado tiene la capacidad de trabajar en 3 dimensiones en multicapas, y en el caso particular en 2 dimensiones se realiza una integración en la vertical. El mallado que se utiliza es de tipo ''staggered cell''  o MAC. En la  [[#fig0005|figura 1]]  se muestran los elementos de la malla, donde el punto negro ubicado al centro de cada celda representa cualquier magnitud escalar como la temperatura y la salinidad, y los círculos en las caras indican la posición de las componentes de los vectores de velocidad ''U'' , ''V''  y ''W'' .

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr1.jpg|center|355px|Discretización en el volumen de control en 2D y 3D.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 1.

Discretización en el volumen de control en 2D y 3D.

</span>
|}

En lo que se refiere a la discretización vertical, se pueden considerar espesores constantes o variables. La distribución de las variables se muestra en la [[#fig0010|figura 2]] .

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr2.jpg|center|355px|Discretización en la vertical.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 2.

Discretización en la vertical.

</span>
|}

El método de solución numérica se basa en una discretización espacial en diferencias finitas en orden 2, y se considera una formulación temporal euleriana-lagrangiana donde los términos advectivos altamente no lineales de las ecuaciones de movimiento y transporte [[#eq0005|(1)]] , [[#eq0010|(2)]] , [[#eq0030|(6)]]  y [[#eq0035|(7)]]  se resuelven por el método de las características (parte lagrangiana) en orden 1, el cual evita oscilaciones numéricas introducidas por fuertes gradientes. Por otro lado, los términos de difusión se resuelven mediante un esquema implícito (parte euleriana) en orden 2. Para la solución discreta de la ecuación de la superficie libre (''η'' ), ecuación  [[#eq0020|(4)]] , es necesario resolver un sistema lineal pentadiagonal (o heptadiagonal en 3D), donde se utiliza el método de Gradiente Conjugado. Los detalles sobre el algoritmo de solución numérica pueden ser consultados en la bibliografía  [[#bib0075|[15]]] , [[#bib0080|[16]]] , [[#bib0090|[18]]] , [[#bib0095|[19]]]  and [[#bib0100|[20]]] .      

==4. Sobre la implementación y la validación del modelo numérico==

El modelo utilizado es validado en 4 aspectos importantes:
* Conceptualmente.
* Algorítmicamente.
* Informáticamente.
* Funcionalmente.

El primero se refiere a que el formulismo matemático que conforman las ecuaciones utilizadas debe incluir los principios físicos fundamentales que garanticen resultados aceptables. El segundo aspecto se refiere al algoritmo de solución utilizado para resolver las ecuaciones; en este caso se utiliza un algoritmo parabólico en el tiempo y elíptico en el espacio, con la implementación de un esquema euleriano-lagrangiano. La parte lagrangiana aproxima los términos advectivos de las ecuaciones de movimiento y transporte, y contempla un algoritmo con delimitador de flujo que ha demostrado eliminar eficazmente las oscilaciones numéricas generadas por estos términos [[#bib0105|[21]]] . La parte euleriana del esquema numérico resuelve con eficiencia los términos difusivos. La conjunción de estos esquemas ha sido validada con soluciones de referencia y casos controlados, y los esquemas han mostrado ser bastante precisos y exactos [[#bib0090|[18]]]  and [[#bib0095|[19]]] . El tercer aspecto hace referencia a la codificación, que es la traducción del algoritmo al lenguaje de programación FORTRAN. Se han utilizado prácticas modernas de programación, como la modularización y la creación de diagramas de flujos de datos. Se han probado por separado los distintos componentes del código en el proceso de su desarrollo, para garantizar la consistencia, la estabilidad y la convergencia de las soluciones, y se han optimizado además los tiempos de solución numérica [[#bib0100|[20]]] .      

Finalmente, el cuarto aspecto se refiere a la validación funcional, la cual consiste en la verificación del modelo con mediciones obtenidas en sitio. La interpretación de un fenómeno físico costero conlleva mediciones con instrumentación que involucran diferentes escalas temporales y espaciales. La medición de la dinámica de estos ecosistemas implica un análisis de estas escalas con métodos particulares, tanto estadísticos como determinísticos. Por ello, la calibración, y por tanto la validación de un modelo en aplicaciones costeras, es bastante compleja ya que, a diferencia de las aplicaciones controladas en laboratorio, los aparatos y mecanismos de medición no son consistentes y traen consigo un cierto grado de aproximación intrínseco. A esto hay que añadirle el factor humano, ya que la interpretación y la apreciación introducen una incertidumbre extra. Para verificar la calidad de la solución numérica con respecto a los datos observados en campo se utiliza la eficiencia de Nash-Sutcliffle [[#bib0110|[22]]]  determinada por la siguiente expresión:

<span id='eq0070'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R^2=1-\frac{\Sigma {\left({\varphi }_{obs}-{\varphi }_{sim}\right)}^2}{\Sigma {\left({\varphi }_{obs}-\overline{\varphi }\right)}^2},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 14)
|}

donde ''φ''<sub>''obs''</sub>  es un dato observado o medido y ''φ''<sub>''sim''</sub>  es un dato obtenido de la simulación numérica, en el mismo punto y en el mismo instante de tiempo. El parámetro <math display="inline">\overline{\varphi }</math>  es la varianza de los datos observados. Evidentemente, si el valor de ''R''  es cercano a la unidad, es posible considerar que los datos de la simulación son bastante aproximados a los datos medidos en campo.      

==5. Caso de estudio de la descarga térmica de la Central Nucleoeléctrica Laguna Verde==

El sistema de enfriamiento de la CNLV es de tipo abierto, el agua destinada para el enfriamiento de los condensadores no se vuelve a utilizar. Aunque es difícil definir los límites de una descarga, para la CNLV se inicia con el aumento de la temperatura que se produce en el condensador y termina cuando el efluente se mezcla y se enfría en el mar. El gasto necesario para enfriar los condensadores de la CNLV es de 63 m<sup>3</sup> /s, lo mismo que se descarga al mar a través de un canal con un incremento de temperatura. La descarga se localiza a 1.800 m de la obra de toma, aproximadamente. A 800 m al sur del canal de descarga se ubica la desembocadura del río El Viejón, como se ilustra en la [[#fig0015|figura 3]] . En esta misma figura se muestra la ubicación geográfica de 6 boyas de medición que se instalaron en la zona para obtener información oceanográfica y meteorológica, que se reportaron en Silva y Botello [[#bib0005|[1]]]  y Botello y Rendón [[#bib0010|[2]]] .

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr3.jpg|center|px|Ubicación geográfica de la obra de toma y descarga de la CNLV (coordenadas ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 3.

Ubicación geográfica de la obra de toma y descarga de la CNLV (coordenadas 19°43′15″''N'' 96°24′23″ O).                  

</span>
|}

El sitio es topográficamente plano, limitado al oeste por montañas de baja elevación (de 500 m aproximadamente) en los primeros 5 km de la costa. Batimétricamente, la zona se caracteriza por una pendiente suave que presenta una isóbata máxima en la parte noreste cuyo valor es de -28 m respecto al nivel medio del mar. En la [[#fig0020|figura 4]]  se muestra el dominio considerado para las simulaciones numéricas y su batimetría, que abarca un área aproximada de 11,7 km<sup>2</sup> , con una longitud de 4.200 m en la dirección ''x''  y 3.500 m en la dirección ''y'' .

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr4.jpg|center|355px|Batimetría de la zona de estudio.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 4.

Batimetría de la zona de estudio.

</span>
|}

Para analizar la dispersión de la pluma térmica de la CNLV se simularon 2 escenarios, los cuales corresponden a los periodos de secas (febrero-mayo) y lluvias (junio-septiembre). Dado que se contó con muy poca información histórica medida en campo en la zona de interés, tanto oceanográfica como meteorológica, se tomaron los datos disponibles de zonas aledañas para diferentes periodos anuales, particularmente de la zona del puerto de Veracruz que se ubica al sur de la CNLV. La información se procesó conjuntamente con los datos medidos en la zona de estudio reportados en Silva y Botello [[#bib0005|[1]]]  y en Botello y Rendón [[#bib0010|[2]]]  para generar las condiciones iniciales y forzantes de cada escenario. La información sobre gastos de succión y gastos y temperaturas de descarga de la CNLV fue proporcionada por la Comisión Federal de Electricidad del Gobierno Mexicano.      

El dominio de estudio se discretizó utilizando una malla rectangular de 88 × 67 celdas, conformando un total de 5.896 elementos, con espaciamientos constantes de Δ''x''  = Δ''y''  = 50 m (véase [[#fig0025|fig. 5]] ). El incremento del paso de tiempo en las simulaciones fue de Δ''t''  = 1 s.

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr5.gif|center|347px|Malla implementada para las simulaciones numéricas, 88×67 celdas.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 5.

Malla implementada para las simulaciones numéricas, 88 × 67 celdas.

</span>
|}

===5.1. Simulaciones en periodo de secas===

Los parámetros de inicialización y forzantes que se impusieron en las simulaciones numéricas se indican en la [[#tbl0005|tabla 1]] . En este caso, no se consideraron los efectos de la marea, la acción del viento sobre la superficie libre, ni el aporte hidrológico del río El Viejón. Los resultados obtenidos con el modelo numérico se validaron cualitativamente con las mediciones en sitio de Silva y Botello [[#bib0005|[1]]] . En la [[#fig0030|figura 6]]  se muestra la pluma térmica obtenida [[#bib0005|[1]]] , donde se observa que esta se dispersa hacia el sur, lejos de la obra de toma de la CNLV. De la misma forma, la secuencia de imágenes de la [[#fig0035|figura 7]] , que ilustra las isotermas y las líneas de corriente obtenidas con el modelo numérico para 120 h de simulación, muestra que la pluma térmica se dispersa con dirección sur, lo que permite concluir que el modelo reproduce el mismo comportamiento que el que se reporta en la referencia con datos medidos en sitio. Además, según estos resultados, la CNLV no se ve afectada por el problema de la recirculación de la pluma térmica hacia la obra de toma.

<span id='tbl0005'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 1.

Datos utilizados como parámetros iniciales para la época de secas

|-

! Parámetro
! Valor
|-

| Corriente litoral
| 0,10 m/s
|-

| Dirección de la corriente litoral
| Norte-Sur
|-

| Temperatura media del mar
| 25 °C
|-

| Gasto de descarga
| 63 m<sup>3</sup> /s                                                    
|-

| Temperatura de la descarga
| 38 °C
|}

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr6.jpg|center|352px|Pluma térmica medida en marzo de 1996 [1].]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 6.

Pluma térmica medida en marzo de 1996 [[#bib0005|[1]]] .                  

</span>
|}

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr7.jpg|center|px|Simulaciones correspondientes al periodo de secas: isotermas (a) 40h de ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 7.

Simulaciones correspondientes al periodo de secas: isotermas (a) 40 h de simulación, (b) 80 h de simulación, (c) 120 h de simulación; líneas de corriente (d) 40 h de simulación, (e) 80 h de simulación y (f) 120 h de simulación.

</span>
|}

===5.2. Simulaciones en periodo de lluvia===

Los valores medios de los parámetros de inicialización y forzantes que se impusieron en las simulaciones numéricas se indican en la [[#tbl0010|tabla 2]] . En este caso se tuvieron en cuenta los efectos de la marea, gracias a los datos proporcionados por las bases de datos del Centro de Investigación Científica y de Educación Superior de Ensenada (CICESE), Baja California, México ([[#fig0040|fig. 8]] ). Se consideraron también los efectos del viento sobre la superficie libre, cuyas magnitudes y direcciones se tomaron de datos medidos en la CNLV en el año 2004 ([[#fig0045|fig. 9]] ). Asimismo, se tuvo en cuenta el flujo de calor en la superficie, caracterizado por el último término del lado derecho de la ecuación [[#eq0030|(6)]]  y definido en la ecuación [[#eq0065|(13)]] . Los datos sobre las corrientes litorales se generaron a partir de información aportada por la Dirección General de Puertos y Marina Mercante, México, para el periodo 1982-1984 ([[#fig0050|fig. 10]] ). Finalmente, se valoró también el aporte hidrológico del río El Viejón, que se estimó de acuerdo con los datos medidos en sitio. El mallado que se utilizó para la simulación numérica cuenta con el mismo número de elementos que en el caso de secas.

<span id='tbl0010'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 2.

Datos utilizados para la simulacion de lluvias

|-

! Parámetro
! Valor
|-

| Corriente litoral
| 0,10 m/s
|-

| Dirección de la corriente litoral
| Sur-Norte
|-

| Temperatura media del mar
| 28 °C
|-

| Salinidad media del mar
| 33''ups''
|-

| Gasto de descarga
| 63 m<sup>3</sup> /s                                                    
|-

| Temperatura de la descarga
| 44 °C
|-

| Gasto del río El Viejón
| 6,09 m<sup>3</sup> /s                                                    
|-

| Temperatura del río El Viejón
| 25 °C
|}

<span id='fig0040'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr8.jpg|center|353px|Marea anual 2005.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 8.

Marea anual 2005.

</span>
|}

<span id='fig0045'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr9.jpg|center|361px|Rosa de los vientos en frecuencias relativas (%), medida en la CNLV, ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 9.

Rosa de los vientos en frecuencias relativas (%), medida en la CNLV, correspondiente al 2004.

</span>
|}

<span id='fig0050'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr10.jpg|center|361px|Rosa de corrientes superficiales en frecuencias relativas (%), del frente ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 10.

Rosa de corrientes superficiales en frecuencias relativas (%), del frente costero de la CNLV, para el periodo de mediciones de 1982 a 1984.

</span>
|}

Los resultados obtenidos con el modelo numérico indican que el escenario de lluvias presenta un fenómeno de recirculación de la pluma térmica en la CNLV. En este caso, debido principalmente a la intensidad y a la dirección de las corrientes, la pluma térmica se dispersa con dirección hacia la obra de toma. Los resultados se ilustran en la [[#fig0055|figura 11]] , en la que puede observarse la evolución de la pluma térmica por medio de las isotermas (''a'' ), (''b'' ) y (''c'' ). La pluma térmica alcanza a la obra de toma y, consecuentemente, afecta al sistema de enfriamiento de la CNLV al presentarse un incremento de temperatura de 2 °C respecto de la temperatura media del mar, justo en la succión. Las líneas de corriente (d), (e) y (f) de la [[#fig0055|figura 11]]  ilustran los campos hidrodinámicos correspondientes a los patrones de dispersión antes mencionados, donde se observa que la orientación de la pluma térmica obedece a la dirección predominante de las corrientes.

<span id='fig0055'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr11.jpg|center|px|Simulaciones correspondientes al periodo de lluvias: isotermas (a) 40h de ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 11.

Simulaciones correspondientes al periodo de lluvias: isotermas (a) 40 h de simulación, (b) 80 h de simulación, (c) 240 h de simulación; líneas de corriente (d) 40 h de simulación, (e) 80 h de simulación, (f) 240 h de simulación.

</span>
|}

En la [[#fig0060|figura 12]]  se muestra la evolución de la pluma térmica obtenida con el modelo numérico. Para ello, se graficaron los valores de temperatura en diferentes puntos ubicados sobre la trayectoria que se esquematiza con la línea negra gruesa, desde la descarga hasta la obra de toma. En esta misma figura, el cero de las abscisas indica el inicio de la trayectoria en la descarga, y la última coordenada, 2,4 km, indica el fin de la trayectoria en la obra de toma. El eje de las ordenadas corresponde a la temperatura. La línea negra fina permite observar que la temperatura de 44 °C impuesta en la descarga decae hasta los 30 °C en la obra de toma.

<span id='fig0060'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr12.jpg|center|362px|Evolución de la temperatura según la trayectoria.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 12.

Evolución de la temperatura según la trayectoria.

</span>
|}

Con la información medida en la obra de toma de la CNLV se pudieron verificar los resultados arrojados por el modelo numérico para el periodo de lluvias. En la [[#fig0065|figura 13]]  se presenta la comparación entre la temperatura medida en la obra de toma (boya 1 en la [[#fig0015|fig. 3]] ) y la temperatura estimada por la simulación en el mismo punto. La señal de temperatura mantiene un comportamiento similar en ambos casos, pues la media de los datos medidos en sitio es de 30,02 °C, mientras que la media de los valores simulados es de 30,26 °C. Con el objetivo de validar y establecer un parámetro cuantitativo de la calidad de la solución numérica se aplicó la ecuación [[#eq0070|(14)]]  para los datos de la [[#fig0065|figura 13]] , y se obtuvo un valor de ''R''  cercano a la unidad igual con ''R''  = 0, 998642, por lo que se considera que la solución numérica tiene una excelente calidad, según el criterio de Nash-Sutcliffle [[#bib0110|[22]]] .

<span id='fig0065'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_740780706-1-s2.0-S0213131513000163-gr13.jpg|center|355px|Comparación de datos medidos (Tobs) contra los datos simulados (Tsim). La línea ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 13.

Comparación de datos medidos (''T''<sub>''obs''</sub> ) contra los datos simulados (''T''<sub>''sim''</sub> ). La línea continua representa los datos históricos de la temperatura medida ''T''<sub>''obs''</sub>  y la línea punteada representa los datos predichos por el modelo.                  

</span>
|}

==6. Conclusiones==

El presente trabajo contribuye al conocimiento del comportamiento de las descargas térmicas en aguas someras mediante el uso de la modelación numérica. El modelo que se implementó para el estudio de la dispersión de la pluma térmica de la Central Nucleoeléctrica Laguna Verde, ubicada en el Estado de Veracruz, México, reprodujo con una excelente calidad de solución 2 escenarios: secas y lluvias. La calidad de los resultados se pudo constatar cualitativamente mediante comparaciones con información de campo y cuantitativamente según el criterio de Nash-Sutcliffle [[#bib0110|[22]]] . Así, se demuestra que la modelación numérica es una herramienta útil en el análisis y en el diseño del funcionamiento hidráulico y ambiental de algunas obras y que, si se utiliza con una correcta apreciación de sus limitaciones intrínsecas, permite reproducir diferentes escenarios en menor tiempo y a bajo costo, con excelentes resultados en comparación con otras técnicas o formas de monitoización y análisis.      

De acuerdo con los resultados obtenidos, el sistema de enfriamiento de la CNLV se ve afectado por la recirculación térmica para el periodo de lluvias ya que, debido principalmente a las corrientes marinas, la pluma térmica descargada llega a la obra de toma con valores medios cercanos a los 30 °C, que son 2 °C por encima de los 28 °C de temperatura media del mar.

==Agradecimientos==

Los autores agradecen a la Comisión Federal de Electricidad, México, y en particular al M. en I. Lázaro Aguilar por las facilidades otorgadas para la realización de este trabajo y al M. en I. Iván Campos Pérez por plantear y participar en el análisis del problema de la recirculación de la pluma térmica de la CNLV.

==References==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[1]]] H.A. Silva, A.V. Botello, editors. Golfo de México. Contaminación e impacto ambiental: diagnóstico y tendencias, chapter Evaluación del Impacto ambiental de la central Nucleoléctrica Laguna Verde. EPOMEX, SEP, 1996.</li>
<li><span id='bib0010'></span>
[[#bib0010|[2]]] A.V. Botello, J. Rendón, editores. Golfo de México. Contaminación e impacto ambiental: diagnóstico y tendencias, chapter Evaluación del Impacto ambiental de la central Nucleoléctrica Laguna Verde a 15 años de operación. EPOMEX, SEP, 2005.</li>
<li><span id='bib0015'></span>
[[#bib0015|[3]]] X. Huang, Z. Zhu, M. Xu, Q. Jing. Variation of water temperature in the southwestern daya bay before and after the operation of Daya Bay nuclear power plant. Annual research reports, Marine biology research station at Daya Bay, Beijing, China, 1998, pp. 102-112.</li>
<li><span id='bib0020'></span>
[[#bib0020|[4]]] Korea Electric Power Corporation A study on the reduction of thermal discharge effects around nuclear power plants. Technical Report No. 92-802, Korea Electric Power Corporation Ltd (KEPCL), 1994.</li>
<li><span id='bib0025'></span>
[[#bib0025|[5]]] E. Sanford, D. Bermudez, M. Bertness, S. Gaines; Flow, food supply and accord barnacle population dynamics; Marine Ecology Progress Series, 104 (1998), pp. 49–62</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] E. Sanford, D. Bermudez, M. Bertness, S. Gaines; Satellite remote sensing of marine pollution; International Journal of Remote Sensing, 14 (1993), pp. 2985–3004</li>
<li><span id='bib0035'></span>
[[#bib0035|[7]]] C.J. Donlon, P.J. Minnett, C. Gentemann, T.J. Nightingale, I.J. Barton, B. Ward,  ''et al.''; Toward improved validation of satellite sea surface skin temperature measurements for climate research; Journal of Climate, 15 (2002), pp. 353–369</li>
<li><span id='bib0040'></span>
[[#bib0040|[8]]] A. Yu-Hwan, S. Palanisamy, L. Jae-Hak, Q.K. Yong; Application of satellite infrared data for mapping of thermal plume contamination in coastal ecosystem of korea; Marine Environmental Research, 61 (2006), pp. 186–201</li>
<li><span id='bib0045'></span>
[[#bib0045|[9]]] D.E. Gibbons, G.E. Wukelic; Application of landsat thematic mapper data for coastal thermal plume analysis at diablo canyon; Photogrammetric Engineering and Remote Sensing, 55 (1989), pp. 903–909</li>
<li><span id='bib0050'></span>
[[#bib0050|[10]]] I.D. James; Modelling pollution dispersion, the ecosystem and water quality in coastal waters: a review; Environmental Modelling and Software, 17 (2002), pp. 363–385</li>
<li><span id='bib0055'></span>
[[#bib0055|[11]]] I. Herrera-Díaz, C. Couder-Castañeda, H. Ramírez-León; A case study of nearshore wave transformation processes along the coast of Mexico near the Laguna Verde nuclear power plant using a fast simulation method; Modelling and Simulation in Engineering (2010)</li>
<li><span id='bib0060'></span>
[[#bib0060|[12]]] I.E. Herrera-Diaz, C. Couder-Castaneda, H. Ramirez Leon; A Fast Simulation Method for Wave Transformation Processes in Coastal Zones; Journal of Applied Research and Technology, 9 (2) (2011), pp. 249–272</li>
<li><span id='bib0065'></span>
[[#bib0065|[13]]] Background papers and supporting data on the international equations of state of seawater. Technical papers in marine science. Technical report, UNESCO, 1981.</li>
<li><span id='bib0070'></span>
[[#bib0070|[14]]] P. Martin, C.N. Barron, L.F. Smedstad, A.J. Wallcraft, R.C. Rhodes, T.J. Campbell, et al. Software design description for the navy coastal ocean model (NCOM) version 4.0. Technical report, Naval Research Lab. Stennis Space Center MS Oceanography Div., 2008.</li>
<li><span id='bib0075'></span>
[[#bib0075|[15]]] C. Rodríguez, E. Serre, C. Rey, H. Ramírez; A numerical model for shallow-water flows: dynamics of the eddy shedding; Wseas Transactions on Environment and Development, 1 (2) (2005)</li>
<li><span id='bib0080'></span>
[[#bib0080|[16]]] C. Rodríguez, E. Serre, C. Rey, H. Ramírez; Numerical investigation of wake formation around a conical island in shallow-water flows; Proceedings of the WSEAS International Conference on Environment, Ecosystems and Development, Venice (Nov 2005), pp. 2–4</li>
<li><span id='bib0085'></span>
[[#bib0085|[17]]] Departamento de Ingeniería Civil Proyecto Laguna Verde. Technical report, Comisión Federal de Electricidad, 2006.</li>
<li><span id='bib0090'></span>
[[#bib0090|[18]]] H. Ramírez, H. Barrios, C. Rodríguez, C. Couder; Baroclinic mathematical modeling of fresh water plumes in the interaction river-sea; International Journal of Numerical Analysis & Modeling, 2 (2005), pp. 1–14</li>
<li><span id='bib0095'></span>
[[#bib0095|[19]]] H. Ramírez, C. Rodríguez, E. Herrera; Multilayer hydrodinamic model and their application to sediment transport in estuaries. Special Issue Shangai Conference, Current Trends in High Performance Computing and Its Applications; Springer-Verlag (2005) pp. 59-70</li>
<li><span id='bib0100'></span>
[[#bib0100|[20]]] C. Couder, H. Ramírez, E. Herrera; Numerical Modeling of Coupled Phenomena in Science and Engineering: Practical Use and Examples, volume 1 of Multiphysics modeling; Taylor & Francis (2008) chapter 34. Numeric optimization of the hydrodynamic model YAXUM/3D</li>
<li><span id='bib0105'></span>
[[#bib0105|[21]]] R. Martin, C. Couder; An improved unsplit and convolutional perfectly matched layer absorbing technique for the navier-stokes equations using cut-off frequency shift; CMES: Computer Modeling in Engineering & Sciences, 63 (1) (2010)</li>
<li><span id='bib0110'></span>
[[#bib0110|[22]]] M.S. Horritt; Computational Fluid Dynamics: Applications in Environmental Hydraulics, chapter 9. Parameterisation, validation and uncertainty analysis of CFD models of fluvial and flood hydraulics in the natural environment; John Wiley and Sons (2005)</li>
</ol>