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_smag- Constante de Smagorinsky (adim)

${\textstyle C_{\rho w}}$- Calor específico del agua (m² /°C s² )

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

g- Aceleración de la gravedad (m/s² )

H- Profundidad total (m)

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

k_p- Constante de ajuste en la ecuación de estado (°C)

K_TH- Coeficiente de dispersión en el plano horizontal (m² C)

K_TV- Coeficiente de dispersión en el plano vertical (m² /s)

ℓ- Longitud de mezclado (m)

P- Presión hidrostática (Pa)

P_atm- Presión atmosférica (Pa)

P_rT- Numero de Prandtl turbulento (adim)

Q- Flujo de calor neto (kg/s³ )

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_b- Distancia vertical desde el fondo hasta un punto cualquiera (m)

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

ν_TH- Coeficiente de viscosidad turbulenta horizontal (m² /s)

ν_TV- Coeficiente de viscosidad turbulenta vertical (m² /s)

ΔA- Área de control (m² )

Δ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³ )

Δ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³ )

ρ₀- Densidad de referencia (kg/m³ )

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 [1]  and [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 [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 [4]  and [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 [6] , [7]  and [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 [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 [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 [1]  and [2] , así como los efectos producidos por el oleaje incidente debido a las obras de protección [11]  and [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:

${\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}}$ $+{\nu }_{TH}\left({\frac {{\partial }^{2}U}{\partial x^{2}}}+{\frac {{\partial }^{2}U}{\partial y^{2}}}\right)+{\frac {\partial }{\partial z}}\left({\nu }_{TV}{\frac {\partial U}{\partial z}}\right)+f_{x},$

( 1)

${\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}}$ $+{\nu }_{TH}\left({\frac {{\partial }^{2}V}{\partial x^{2}}}+{\frac {{\partial }^{2}V}{\partial y^{2}}}\right)+{\frac {\partial }{\partial z}}\left({\nu }_{TV}{\frac {\partial V}{\partial z}}\right)-f_{y}.$

( 2)

Ecuación de continuidad:

${\frac {\partial U}{\partial x}}+{\frac {\partial V}{\partial y}}+$${\frac {\partial W}{\partial z}}=0.$

( 3)

Ecuación para la superficie libre:

${\frac {\partial U}{\partial x}}+{\frac {\partial }{\partial x}}\left({\int }_{z}^{\eta }Udz\right)+$${\frac {\partial }{\partial y}}\left({\int }_{z}^{\eta }Vdz\right)=$$0.$

( 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 [13] :

$\rho (S,T,P)={\frac {{\rho }_{0}}{1-{\frac {P}{k_{p}}}}}.$

( 5)

donde k_p es una constante de ajuste que se determina por módulos [13] , y es función de las variaciones de P , T y S .

Las ecuaciones para la temperatura y la salinidad son:

${\frac {\partial T}{\partial t}}+U{\frac {\partial T}{\partial x}}+V{\frac {\partial T}{\partial y}}+W{\frac {\partial T}{\partial z}}=$ $+K_{TH}\left({\frac {{\partial }^{2}T}{\partial x^{2}}}+{\frac {{\partial }^{2}T}{\partial y^{2}}}\right)+{\frac {\partial }{\partial z}}\left(K_{TV}{\frac {\partial T}{\partial z}}\right)\pm Q,$

( 6)

${\frac {\partial S}{\partial t}}+U{\frac {\partial S}{\partial x}}+V{\frac {\partial S}{\partial y}}+W{\frac {\partial S}{\partial z}}=$ $+K_{TH}\left({\frac {{\partial }^{2}S}{\partial x^{2}}}+{\frac {{\partial }^{2}S}{\partial y^{2}}}\right)+{\frac {\partial }{\partial z}}\left(K_{TV}{\frac {\partial S}{\partial z}}\right).$

( 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 (3) como:

${\frac {\partial W}{\partial z}}=-\left({\frac {\partial U}{\partial x}}+\right.$$\left.{\frac {\partial V}{\partial y}}\right).$

( 8)

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

$-{\frac {1}{{\rho }_{0}}}{\frac {\partial P}{\partial x}}=$$-{\frac {\rho g}{{\rho }_{0}}}{\frac {\partial \eta }{\partial x}}-$${\frac {g}{{\rho }_{0}}}{\int }_{z}^{\eta }{\frac {\partial {\rho }^{'}}{\partial z}}dz-$${\frac {1}{{\rho }_{0}}}{\frac {\partial P_{atm}}{\partial x}}.$

( 9)

El primer término de la derecha de la ecuación (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 (2) .

2.1. Los coeficientes de viscosidad turbulenta

El coeficiente de viscosidad turbulenta horizontal ν_TH que aparece en las ecuaciones (1) y (2) se estima mediante el modelo de Smagorinsky como sigue:

${\nu }_{TH}=C_{smag}\Delta x\Delta y$ $\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],$

( 10)

donde C_smag 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 [14] . Para el cálculo del coeficiente de viscosidad turbulenta vertical ν_TV se utiliza el modelo siguiente [15] :

${\nu }_{TV}=l{\sqrt {{\left({\frac {\partial U}{\partial z}}\right)}^{2}+{\left({\frac {\partial V}{\partial z}}\right)}^{2}}},$

( 11)

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

$l={\begin{array}{lll}k\cdot z_{b}&si&z_{b}<0,20H\\0,20H\cdot k&si&z_{b}>0,20H,\end{array}}$

( 12)

donde k es la constante de Von Kármán igual a 0,40 y z_b es la distancia vertical desde el fondo hasta el punto donde se calcula ν_TV .

Los coeficientes de difusión de temperatura y salinidad de las ecuaciones (6) y (7) son evaluados como ${\textstyle K_{TH}=}$${\frac {{\nu }_{TH}}{{Pr}_{T}}}$ y ${\textstyle K_{TV}=}$${\frac {{\nu }_{TV}}{{Pr}_{T}}}$ , donde Pr_T es el número de Prandtl turbulento igual a 0,9 [15]  and [16] .

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

En la ecuación (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 [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:

$Q={\frac {\Delta Q(\Delta A)}{\rho C_{\rho w}(\Delta V)}}.$

( 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 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 .

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

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 figura 2 .

Figura 2. Discretización en la vertical.

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 (1) , (2) , (6) y (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 (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 [15] , [16] , [18] , [19]  and [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 [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 [18]  and [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 [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 [22] determinada por la siguiente expresión:

$R^{2}=1-{\frac {\Sigma {\left({\varphi }_{obs}-{\varphi }_{sim}\right)}^{2}}{\Sigma {\left({\varphi }_{obs}-{\overline {\varphi }}\right)}^{2}}},$

( 14)

donde φ_obs es un dato observado o medido y φ_sim es un dato obtenido de la simulación numérica, en el mismo punto y en el mismo instante de tiempo. El parámetro ${\textstyle {\overline {\varphi }}}$ 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³ /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 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 [1] y Botello y Rendón [2] .

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).

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 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² , con una longitud de 4.200 m en la dirección x y 3.500 m en la dirección y .

Figura 4. Batimetría de la zona de estudio.

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 [1] y en Botello y Rendón [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 fig. 5 ). El incremento del paso de tiempo en las simulaciones fue de Δt  = 1 s.

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

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 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 [1] . En la figura 6 se muestra la pluma térmica obtenida [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 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.

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 m3 /s
Temperatura de la descarga	38 °C

Figura 6. Pluma térmica medida en marzo de 1996 [1] .

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.

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 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 (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 (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 (6) y definido en la ecuación (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 (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.

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	33ups
Gasto de descarga	63 m3 /s
Temperatura de la descarga	44 °C
Gasto del río El Viejón	6,09 m3 /s
Temperatura del río El Viejón	25 °C

Figura 8. Marea anual 2005.

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

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

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 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 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.

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.

En la 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.

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

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 figura 13 se presenta la comparación entre la temperatura medida en la obra de toma (boya 1 en la 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 (14) para los datos de la 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 [22] .

Figura 13. Comparación de datos medidos (Tobs ) contra los datos simulados (Tsim ). La línea continua representa los datos históricos de la temperatura medida Tobs y la línea punteada representa los datos predichos por el modelo.

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 [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

1. [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.
2. [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.
3. [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.
4. [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.
5. [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
6. [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
7. [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
8. [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
9. [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
10. [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
11. [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)
12. [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
13. [13] Background papers and supporting data on the international equations of state of seawater. Technical papers in marine science. Technical report, UNESCO, 1981.
14. [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.
15. [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)
16. [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
17. [17] Departamento de Ingeniería Civil Proyecto Laguna Verde. Technical report, Comisión Federal de Electricidad, 2006.
18. [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
19. [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
20. [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
21. [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)
22. [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)