Numerical Modeling of Water Thermal Plumes Emitted by Thermal Power Plants

This work focuses on the study of thermal dispersion of plumes emitted by power plants into the sea. Wastewater discharge from power stations causes impacts that require investigation or monitoring. A study to characterize the physical effects of thermal plumes into the sea is carried out here by numerical modeling and field measurements. The case study is the thermal discharges of the Presidente Adolfo López Mateos Power Plant, located in Veracruz, on the coast of the Gulf of Mexico. This plant is managed by the Federal Electricity Commission of Mexico. The physical effects of such plumes are related to the increase of seawater temperature caused by the hot water discharge of the plant. We focus on the implementation, calibration, and validation of the Delft3D-FLOW model, which solves the shallow-water equations. The numerical simulations consider a critical scenario where meteorological and oceanographic parameters are taken into account to reproduce the proper physical conditions of the environment. The results show a local physical effect of the thermal plumes within the study zone, given the predominant strong winds conditions of the scenario under study.


Introduction
In recent years, coastal water bodies have experienced extensive environmental impact.Despite regulations, sea pollution in coastal areas has continued to alter physical, chemical, and biological water properties [1,2].The main source of pollution at the coast is industrial activity and specifically the discharge from power plants.
According to the U.S. Geological Survey, power plants use more water than any other industry [3].The main source of electricity in Mexico is thermal power plants that use water for cooling and then discharge that water to a heat sink at an increased temperature.The heat sink could be the ocean.When discharged, the heated water mixes into the water body to generate a thermal plume, which can cause environmental impacts [4].This anomaly in temperature can significantly disrupt the aquatic ecosystem [5].On the other hand, the efficiency of the power plant's cooling systems can be compromised if the thermal plume recirculates through the intake.For these reasons, the thermal dispersion of such plumes must be evaluated to identify the main driving parameters and the influence area within the water body.
The monitoring of the thermal plume of power plants has been carried out by several techniques, e.g., measurements with oceanographic buoys and cruises, satellites, and numerical modeling.The use of satellites is a well-developed technique for monitoring and recording surface water temperature and provides a synoptic view of the temperature in large water bodies [6,7].The use of numerical modeling for water quality studies has increased, especially as a complement to the other techniques.The viability of numerical modeling is mainly due to both the progress in multidisciplinary modeling in recent decades and the progress in computing technologies and programming performance, which makes possible the computation of complex phenomena with high resolution and reduced processing times.Some recent works that study thermal plume dispersion, specifically from power plants, through numerical modeling can be consulted in [4] and [8,9].
This paper presents a study of a thermal plume dispersion into the sea caused by a power plant.The methodology is based on field measurements and numerical simulations.The case study is the thermal plume of the Presidente Adolfo López Mateos Power Plant (CTPALM), operated by the Federal Electricity Commission of Mexico (CFE).For the CFE, it is important to monitor the plume dispersion to better manage the possible physical effects on the coast caused by the temperature increase of the seawater.Predictions from numerical simulations are compared against measurements for calibration and validation.The results show the influence area and the directionality of the thermal plume under a specific extreme weather condition.

Case Study
The CTPALM is one of the most important power plants of Mexico, producing 8000 GWh per year approximately [10].The plant uses about 90 m 3 /s of seawater of the Gulf of Mexico (GM) for cooling, then the water is returned to the GM with a temperature increase of about 5 • C through four discharge channels.Figure 1 shows the CTPALM, the location of its four discharges (marked as D1, D2, D3, and D4), and the intake structure (IS).The power plant is situated 18 km from Tuxpan City, Veracruz, Mexico, in a region where a river-lagoon-sea system interacts and comprises the discharge of the Tuxpan River, the Tampamachoco Lagoon, and the coast [11].
water temperature and provides a synoptic view of the temperature in large water bodies [6,7].The use of numerical modeling for water quality studies has increased, especially as a complement to the other techniques.The viability of numerical modeling is mainly due to both the progress in multidisciplinary modeling in recent decades and the progress in computing technologies and programming performance, which makes possible the computation of complex phenomena with high resolution and reduced processing times.Some recent works that study thermal plume dispersion, specifically from power plants, through numerical modeling can be consulted in [4] and [8,9].
This paper presents a study of a thermal plume dispersion into the sea caused by a power plant.The methodology is based on field measurements and numerical simulations.The case study is the thermal plume of the Presidente Adolfo López Mateos Power Plant (CTPALM), operated by the Federal Electricity Commission of Mexico (CFE).For the CFE, it is important to monitor the plume dispersion to better manage the possible physical effects on the coast caused by the temperature increase of the seawater.Predictions from numerical simulations are compared against measurements for calibration and validation.The results show the influence area and the directionality of the thermal plume under a specific extreme weather condition.

Case Study
The CTPALM is one of the most important power plants of Mexico, producing 8000 GWh per year approximately [10].The plant uses about 90 m 3 /s of seawater of the Gulf of Mexico (GM) for cooling, then the water is returned to the GM with a temperature increase of about 5 °C through four discharge channels.Figure 1 shows the CTPALM, the location of its four discharges (marked as D1, D2, D3, and D4), and the intake structure (IS).The power plant is situated 18 km from Tuxpan City, Veracruz, Mexico, in a region where a river-lagoon-sea system interacts and comprises the discharge of the Tuxpan River, the Tampamachoco Lagoon, and the coast [11].

Field Measurements
Information about the physical parameters prevailing in the study zone was obtained by field investigation and bibliographic compilation.Some information was provided by the CFE and the

Field Measurements
Information about the physical parameters prevailing in the study zone was obtained by field investigation and bibliographic compilation.Some information was provided by the CFE and the National Water Commission of Mexico (CONAGUA).The field investigation was conducted in October 2011, when the prevailing weather conditions correspond to a season known as Nortes in Mexico.In these latitudes, the Nortes season coincides with winter, approximately from October to February.This season is characterized by a series of cold fronts that generate strong winds and sudden changes in weather, especially along the coast of the GM.In this region, the force of a cold front is able to drive the mass of water in the direction of its movement, which is generally from northwest to southeast.A decrease in the air temperature and the seawater temperature is also an effect of the passage of a cold front.
The field investigation was planned to provide real-time and long-term measurements.Various instruments were installed inshore just in front of the power plant at many observation points.The parameters measured were currents, tides, seawater temperature, and bathymetry.The data were used to identify the driving forces of the study area and to implement the numerical model.
To measure currents and sea level variations, an Acoustic Doppler Profiler (ADCP) was used.The instrument is a SonTek ADCP, high-performance, 3D water current profiler, with a frequency of 500 kHz, equipped with a pressure sensor.The ADCP was mounted on the seafloor at a depth of 12.15 m, about 1 km seaward, directly in front of Discharge 3 (see Figure 2).
National Water Commission of Mexico (CONAGUA).The field investigation was conducted in October 2011, when the prevailing weather conditions correspond to a season known as Nortes in Mexico.In these latitudes, the Nortes season coincides with winter, approximately from October to February.This season is characterized by a series of cold fronts that generate strong winds and sudden changes in weather, especially along the coast of the GM.In this region, the force of a cold front is able to drive the mass of water in the direction of its movement, which is generally from northwest to southeast.A decrease in the air temperature and the seawater temperature is also an effect of the passage of a cold front.
The field investigation was planned to provide real-time and long-term measurements.Various instruments were installed inshore just in front of the power plant at many observation points.The parameters measured were currents, tides, seawater temperature, and bathymetry.The data were used to identify the driving forces of the study area and to implement the numerical model.
To measure currents and sea level variations, an Acoustic Doppler Profiler (ADCP) was used.The instrument is a SonTek ADCP, high-performance, 3D water current profiler, with a frequency of 500 kHz, equipped with a pressure sensor.The ADCP was mounted on the seafloor at a depth of 12.15 m, about 1 km seaward, directly in front of Discharge 3 (see Figure 2).The bathymetry along the coast in front of the power plant was measured using an echo sounder with GPS, for points distributed in a rectangular mesh.The depth measured by the sounder was corrected for the astronomical tidal conditions.Figure 3 shows the bathymetry of the study zone, which is characterized by contour lines parallel to the shoreline, with a relatively smooth gradient perpendicular to the shore.The bathymetry along the coast in front of the power plant was measured using an echo sounder with GPS, for points distributed in a rectangular mesh.The depth measured by the sounder was corrected for the astronomical tidal conditions.Figure 3 shows the bathymetry of the study zone, which is characterized by contour lines parallel to the shoreline, with a relatively smooth gradient perpendicular to the shore.Time series of seawater temperature were measured at two different points located in front of the power plant.A mooring with two thermistors was installed at each point, one thermistor at 3 m depth (surface) and the other at 7 m depth (bottom), approximately.The location of each mooring is shown in Figure 2. To observe the influence of the four thermal discharges of the CTPALM in the sea at the sea surface, four thermal traces were performed using boat-based CTD on 14 October 2011, during the passage of a cold front.The first trace was oriented parallel to the shoreline from the Tuxpan River discharge to a point near Discharge 3; the second trace was also oriented parallel to the shoreline starting at a point near Discharge 3 and extending approximately 750 m northward; the third trace was oriented perpendicular to the shoreline starting at a point near Discharge 3 and extending seaward approximately 700 m; and the fourth trace was also oriented perpendicular to the shoreline starting further out in the sea and returning coastward to a point near Discharge 1 (see Figure 4).These measurements were taken up to 0.1 m depth at most, in the course of the traces.

Numerical Model
The numerical model used here is the Delft3D-FLOW model.This model was developed by Deltares Institute of Netherlands and has been applied around the world for many studies.The numerical model simulates unsteady flow in water bodies and can be applied in a two-dimensional Time series of seawater temperature were measured at two different points located in front of the power plant.A mooring with two thermistors was installed at each point, one thermistor at 3 m depth (surface) and the other at 7 m depth (bottom), approximately.The location of each mooring is shown in Figure 2. To observe the influence of the four thermal discharges of the CTPALM in the sea at the sea surface, four thermal traces were performed using boat-based CTD on 14 October 2011, during the passage of a cold front.The first trace was oriented parallel to the shoreline from the Tuxpan River discharge to a point near Discharge 3; the second trace was also oriented parallel to the shoreline starting at a point near Discharge 3 and extending approximately 750 m northward; the third trace was oriented perpendicular to the shoreline starting at a point near Discharge 3 and extending seaward approximately 700 m; and the fourth trace was also oriented perpendicular to the shoreline starting further out in the sea and returning coastward to a point near Discharge 1 (see Figure 4).These measurements were taken up to 0.1 m depth at most, in the course of the traces.Time series of seawater temperature were measured at two different points located in front of the power plant.A mooring with two thermistors was installed at each point, one thermistor at 3 m depth (surface) and the other at 7 m depth (bottom), approximately.The location of each mooring is shown in Figure 2. To observe the influence of the four thermal discharges of the CTPALM in the sea at the sea surface, four thermal traces were performed using boat-based CTD on 14 October 2011, during the passage of a cold front.The first trace was oriented parallel to the shoreline from the Tuxpan River discharge to a point near Discharge 3; the second trace was also oriented parallel to the shoreline starting at a point near Discharge 3 and extending approximately 750 m northward; the third trace was oriented perpendicular to the shoreline starting at a point near Discharge 3 and extending seaward approximately 700 m; and the fourth trace was also oriented perpendicular to the shoreline starting further out in the sea and returning coastward to a point near Discharge 1 (see Figure 4).These measurements were taken up to 0.1 m depth at most, in the course of the traces.

Numerical Model
The numerical model used here is the Delft3D-FLOW model.This model was developed by Deltares Institute of Netherlands and has been applied around the world for many studies.The numerical model simulates unsteady flow in water bodies and can be applied in a two-dimensional

Numerical Model
The numerical model used here is the Delft3D-FLOW model.This model was developed by Deltares Institute of Netherlands and has been applied around the world for many studies.The numerical model simulates unsteady flow in water bodies and can be applied in a two-dimensional or three-dimensional mode, depending on the vertical structure of the flow under study.The system of equations consists of the horizontal equations of motion, the continuity equation, and the transport equations for conservative constituents.In the horizontal direction, Delft3D-FLOW uses orthogonal curvilinear co-ordinates whereas in the vertical direction uses a σ or a Z co-ordinate system.In combination with an appropriate set of initial and boundary conditions, the system of equations is solved on a finite difference staggered grid.Specific details about the time integration method, the particular spatial discretization of the advection terms, turbulence modeling, and more model scopes can be consulted in [12][13][14][15][16][17][18][19][20], respectively.Moreover, some recent similar works to the present one can be consulted in [21][22][23].

Implementation of the Numerical Model
The computational domain implemented for the numerical simulations is shown in Figure 2. The surface covers a coastal front of 9 km and extends 3.7 km seaward.The intake structure of the CTPALM, together with its four discharges, is located on the coastal line boundary.The northern, southern, and eastern boundaries are considered sea open boundaries.The discrete domain consists of a two-dimensional curvilinear grid of 290 × 128 cell elements (see Figure 2).The use of a two-dimensional grid is justified by a barotropic vertical flow structure, observed in the study zone with the aid of the ADCP measurements (see Section 3.1).The elements are quasi-rectangular with an approximate size of 30 m in each horizontal dimension.A maximum aspect ratio of 1.5 is obtained in some regions near the coastal line.
For the numerical simulations the sea level variation measured with the ADCP is imposed at the eastern open boundary and free flow is considered at the northern and southern open boundaries.The initial velocities are set to zero and the temperature is imposed uniform according with the data of the thermistors.This same initial temperature is also imposed at the sea open boundaries.The constant flow rates and constant temperatures considered in the simulation for both the intake structure and the four discharges were provided by the CFE.The physical parameters used for the numerical simulations are shown in Table 1.Winds are assumed to be uniformly distributed over the study area and vary in accordance with data reported for the study period by the weather station 766400 (owned by CONAGUA), located in Tuxpan City at Lat. 20.95 • N and Long.97.4 • W. Other climatological parameters such as air temperature and relative humidity were also obtained from the weather station and are used in the Delft3F-FLOW heat flux model.
The time period for the numerical simulations was defined with the aid of the ADCP velocity measurements and existing climatological data.The simulation time period began on 13 October and ended on 15 October 2011.During this period, the cold front Number 6 of the 2011 Nortes season was reported by the CFE.This meteorological phenomenon affected the nearshore current direction, which completely changed from northwest to southeast during the occurrence of this event, as discussed in Section 3.1.This scenario was selected for this study because it is an extraordinary event characterized by strong winds and occurred within a relatively short time window, causing sudden weather changes.The discrete time step is set to 12 s for all the numerical simulations.

Results of the Field Measurements
The ADCP measurements of the flow velocity time series at different layers and sea level variations are shown in Figure 5. Data of the ADCP reveal two significant features: the similarity of the currents through the water column, and the directionality of flow.Regarding the similarity, vertical homogeneity is evident and can be attributed to both the shallow depth of the zone and the absence or low impact of water mass flow with different physical or chemical properties.Regarding the directionality, two naturally preferred directions were observed, northwest and southeast, presumably attributable to the shoreline alignment.After a correlation analysis between the velocity of the upper layer and the velocity of the lower layer, a barotropic structure can be deduced (see Figure 6).Regarding the variability of cycle rates of the currents, a signal analysis shows that the velocity series are dominated by low frequencies, much lower than the frequency of the astronomical tide, where the higher one is around 0.016 cycles/day of a period of 25.6 days (see Figure 7).Accordingly, during the period of measurements the evolution of the currents is not dependent on tides; it must be attributable to processes having lower frequencies such as regional mass density flow or the recurrence of cold fronts.event characterized by strong winds and occurred within a relatively short time window, causing sudden weather changes.The discrete time step is set to 12 s for all the numerical simulations.

Results of the Field Measurements
The ADCP measurements of the flow velocity time series at different layers and sea level variations are shown in Figure 5. Data of the ADCP reveal two significant features: the similarity of the currents through the water column, and the directionality of flow.Regarding the similarity, vertical homogeneity is evident and can be attributed to both the shallow depth of the zone and the absence or low impact of water mass flow with different physical or chemical properties.Regarding the directionality, two naturally preferred directions were observed, northwest and southeast, presumably attributable to the shoreline alignment.After a correlation analysis between the velocity of the upper layer and the velocity of the lower layer, a barotropic structure can be deduced (see Figure 6).Regarding the variability of cycle rates of the currents, a signal analysis shows that the velocity series are dominated by low frequencies, much lower than the frequency of the astronomical tide, where the higher one is around 0.016 cycles/day of a period of 25.6 days (see Figure 7).Accordingly, during the period of measurements the evolution of the currents is not dependent on tides; it must be attributable to processes having lower frequencies such as regional mass density flow or the recurrence of cold fronts.In Figure 8, the time series of temperature measured with thermistors are shown.The time series show some deviation between the surface and the bottom temperatures the first days of measurements, which can be attributed to the cold fronts.However, the temperature measurements confirm that the flow structure is adequately mixed, because there are no significant differences between the surface and the bottom temperatures.A gradual decrease of the seawater temperature, characteristic of the transition between autumn and winter, was observed in the time series data.In Figure 8, the time series of temperature measured with thermistors are shown.The time series show some deviation between the surface and the bottom temperatures the first days of measurements, which can be attributed to the cold fronts.However, the temperature measurements confirm that the flow structure is adequately mixed, because there are no significant differences between the surface and the bottom temperatures.A gradual decrease of the seawater temperature, characteristic of the transition between autumn and winter, was observed in the time series data.In Figure 8, the time series of temperature measured with thermistors are shown.The time series show some deviation between the surface and the bottom temperatures the first days of measurements, which can be attributed to the cold fronts.However, the temperature measurements confirm that the flow structure is adequately mixed, because there are no significant differences between the surface and the bottom temperatures.A gradual decrease of the seawater temperature, characteristic of the transition between autumn and winter, was observed in the time series data.In Figure 8, the time series of temperature measured with thermistors are shown.The time series show some deviation between the surface and the bottom temperatures the first days of measurements, which can be attributed to the cold fronts.However, the temperature measurements confirm that the flow structure is adequately mixed, because there are no significant differences between the surface and the bottom temperatures.A gradual decrease of the seawater temperature, characteristic of the transition between autumn and winter, was observed in the time series data.The sea surface temperature measured with CTD along each of the four traces of Figure 4 is shown in Figure 9.In each of the four traces, a temperature increase can be observed as the boat approaches the point of discharge.The first trace shows a minimum temperature of about 28.4 • C and a peak of 31.3 • C, which yields a maximum observed temperature increase of 2.9 • C (see Figure 9a).The peak observed in first trace corresponds to the heat of Discharge 1, which is the most important source of heat of the power plant.The second peak of about 30.2 • C is evidence of the heat emitted by Discharge 3, the second most important source of heat.The other traces show lower temperature increases of about 1.8 • C (see Figure 9b-d), reaching temperatures up to 30.9 • C.These last three traces hardly shown the heat emitted by Discharges 3. The traces perpendicular to the shoreline have a quasi-linear temperature decay seaward for the first 1000 m (see Figure 9c,d).Based on the temperature decay lengths observed along the traces with CDT, an influence area of the thermal plumes can be approximated.By considering a reference seawater temperature of 28.66 • C, which is the average of the temperature measured with thermistors TM 1 and TM 2 during the simulation time period, the first and second traces show a decay length of about 2840 m, parallel to the shoreline, while the third and fourth traces show a decay length of about 1000 m, perpendicular to the shoreline.This yields a rectangular surface of 2.84 km 2 , which is assumed to be the approximate sea surface area of the thermal plumes.
The sea surface temperature measured with CTD along each of the four traces of Figure 4 is shown in Figure 9.In each of the four traces, a temperature increase can be observed as the boat approaches the point of discharge.The first trace shows a minimum temperature of about 28.4 °C and a peak of 31.3 °C, which yields a maximum observed temperature increase of 2.9 °C (see Figure 9a).The peak observed in first trace corresponds to the heat of Discharge 1, which is the most important source of heat of the power plant.The second peak of about 30.2 °C is evidence of the heat emitted by Discharge 3, the second most important source of heat.The other traces show lower temperature increases of about 1.8 °C (see Figure 9b-d), reaching temperatures up to 30.9 °C.These last three traces hardly shown the heat emitted by Discharges 1 and 3.The traces perpendicular to the shoreline have a quasi-linear temperature decay seaward for the first 1000 m (see Figure 9c,d).
Based on the temperature decay lengths observed along the traces with CDT, an influence area of the thermal plumes can be approximated.By considering a reference seawater temperature of 28.66 °C, which is the average of the temperature measured with thermistors TM 1 and TM 2 during the simulation time period, the first and second traces show a decay length of about 2840 m, parallel to the shoreline, while the third and fourth traces show a decay length of about 1000 m, perpendicular to the shoreline.This yields a rectangular surface of 2.84 km 2 , which is assumed to be the approximate sea surface area of the thermal plumes.

Calibration and Validation of the Numerical Simulations
Numerical simulations were performed to calibrate and validate the numerical model.The calibration of the dynamic field consisted of a trial wind drag coefficient C d and Manning coefficient n, adjusted in subsequent simulations until the pattern of the predicted depth-averaged horizontal velocities best matched that measured with ADCP over the simulation time period.There is not a specific strategy to define the values of C d ; it is a trial and error process in which each result was compared against measurements and then validated.The Manning coefficient has been prescribed according to the literature for different boundary types.Both coefficients were considered to be uniform over the study area and the wind drag coefficient was imposed as a function of wind speed.The parameters used in the calibration tests are presented in Table 2.The predictions of the velocity components were stored at the discrete point shown in Figure 2, which corresponds to the location of the ADCP.In Figure 10, the evolution of measured horizontal velocities and the corresponding predictions are graphed.For the case of V-velocity, some tests over-predict or under-predict the peak, whereas others yield matching velocities at the end of the period.For the case of U-velocity, most of the tests accurately predicted the pattern, though the velocity is under-predicted in some tests.
Water 2016, 8, 482 9 of 16 is not a specific strategy to define the values of d C ; it is a trial and error process in which each result was compared against measurements and then validated.The Manning coefficient has been prescribed according to the literature for different boundary types.Both coefficients were considered to be uniform over the study area and the wind drag coefficient was imposed as a function of wind speed.The parameters used in the calibration tests are presented in Table 2.The predictions of the velocity components were stored at the discrete point shown in Figure 2, which corresponds to the location of the ADCP.In Figure 10, the evolution of measured horizontal velocities and the corresponding predictions are graphed.For the case of V -velocity, some tests over-predict or under-predict the peak, whereas others yield matching velocities at the end of the period.For the case of U -velocity, most of the tests accurately predicted the pattern, though the velocity is under-predicted in some tests.The validation of the dynamic field was effectuated over the depth-averaged velocity vector using the following efficiency criteria: The validation of the dynamic field was effectuated over the depth-averaged velocity vector V = √ U 2 + V 2 using the following efficiency criteria: where RMSE is the root-mean-square error, m is the mass balance error, the indexes M and P denote measured and predicted values, respectively, and k is the total number of data points.The different coefficients obtained in each calibration test are shown in Table 3, where the range of possible values and the perfect match value are specified for both criteria in the two last columns.According to the efficiency criteria, some predictions are poorly consistent with measurements.However, despite the complexity of the case study where diverse parameters interact, the flow directionality is well predicted in all calibration tests, as Figure 10 shows.Based on both efficiency criteria, test 4 better approximates the measurements, with RMSE = 0.087 and m = 0.360.Thus, the wind drag coefficients C d and the Manning coefficient n of test 4 were retained for simulations (see Table 2).
Both efficiency criteria, the root-mean-square error and the mass balance error, were also applied to compare the temperature measured with the thermistors and the depth-averaged predicted temperature.The time series between the surface and the bottom temperatures measured at each point TM 1 and TM 2 (see Figure 8) were averaged to be able to compare with the depth-averaged predicted temperature.Figure 11 shows the evolution of depth-averaged measured and predicted temperatures and Table 4 summarizes the calculated coefficients over the simulation time period.Although both efficiency criteria provide good results, the evolution of measured and predicted temperatures shows that the numerical simulation tends to overestimate the temperature in both points TM 1 and TM 2. The point where the highest differences were obtained is TM 1, but these differences do not exceed 1 • C. At point TM 2, these differences are slightly lower.It is difficult to determine whether the predictions of temperature are good or bad given the different sources of uncertainty involved in the study case, such as the precision of the thermistors or the quality of the meteorological data considered as boundary conditions for the heat flux model.Thus, a maximum range of 1 • C in the differences between measured and predicted temperatures may be considered acceptable.
where RMSE is the root-mean-square error, m is the mass balance error, the indexes M and P denote measured and predicted values, respectively, and k is the total number of data points.
The different coefficients obtained in each calibration test are shown in Table 3, where the range of possible values and the perfect match value are specified for both criteria in the two last columns.According to the efficiency criteria, some predictions are poorly consistent with measurements.However, despite the complexity of the case study where diverse parameters interact, the flow directionality is well predicted in all calibration tests, as Figure 10 shows.Based on both efficiency criteria, test 4 better approximates the measurements, with RMSE = 0.087 and m = 0.360.Thus, the wind drag coefficients d C and the Manning coefficient n of test 4 were retained for simulations (see Table 2).
Both efficiency criteria, the root-mean-square error and the mass balance error, were also applied to compare the temperature measured with the thermistors and the depth-averaged predicted temperature.The time series between the surface and the bottom temperatures measured at each point TM 1 and TM 2 (see Figure 8) were averaged to be able to compare with the depth-averaged predicted temperature.Figure 11 shows the evolution of depth-averaged measured and predicted temperatures and Table 4 summarizes the calculated coefficients over the simulation time period.Although both efficiency criteria provide good results, the evolution of measured and predicted temperatures shows that the numerical simulation tends to overestimate the temperature in both points TM 1 and TM 2. The point where the highest differences were obtained is TM 1, but these differences do not exceed 1 °C.At point TM 2, these differences are slightly lower.It is difficult to determine whether the predictions of temperature are good or bad given the different sources of uncertainty involved in the study case, such as the precision of the thermistors or the quality of the meteorological data considered as boundary conditions for the heat flux model.Thus, a maximum range of 1 °C in the differences between measured and predicted temperatures may be considered acceptable.

Thermal Plume Dispersion
Figure 12 shows the current velocity vectors and magnitudes at time intervals of 12 h.As can be observed, the currents maintain a predominantly southward direction, which means that the driven force is the wind acting on the free surface during the passage of the cold front.During the simulation time period, the current velocities reach magnitudes up to 0.4 m/s in the offshore region, whereas the maximum velocities occur near the discharges of the power plant, reaching values greater than 1.0 m/s.In this case, the effects of tides or typical regional mass density flow on the hydrodynamics are not significant.As the winds weaken, these mechanisms have more and more influence on the currents.Other existing studies [9] have shown the significant effects of the tide, which differs from the results observed for the present specific case study.
The thermal plume dispersion is highly influenced by the magnitude and direction of the currents.Figure 13 illustrates the evolution of the predicted thermal field for the same time intervals in Figure 12.The four thermal discharges of the power plant can be clearly observed with the aid of the isotherms.Because of the predominant flow direction toward the southeast, the thermal plumes are dispersed in the same general direction.Moreover, the plumes are dispersed along the coast, at least during the cold front, for the first 24 h of the simulation time period (see Figure 13a-c).When the cold front weakens at the end of the simulation time period, a seaward thermal dispersion is observed.It is evident that the influence area of the thermal plumes of Discharges 1 and 3 is larger than the area covered by the thermal plumes of Discharges 2 and 4.This is likely because the discharge rate at Discharge 1 and 3 is significantly greater than for Discharge 2 and 4. The directionality developed by the thermal plume dispersion is undesirable from the point of view of the efficiency of the power plant, because large amounts of hot water emitted from Discharges 1 and 3 could recirculate toward the intake.As mentioned above, the efficiency of the power plant is compromised while hot water is taken for cooling.

Thermal Plume Dispersion
Figure 12 shows the current velocity vectors and magnitudes at time intervals of 12 h.As can be observed, the currents maintain a predominantly southward direction, which means that the driven force is the wind acting on the free surface during the passage of the cold front.During the simulation time period, the current velocities reach magnitudes up to 0.4 m/s in the offshore region, whereas the maximum velocities occur near the discharges of the power plant, reaching values greater than 1.0 m/s.In this case, the effects of tides or typical regional mass density flow on the hydrodynamics are not significant.As the winds weaken, these mechanisms have more and more influence on the currents.Other existing studies [9] have shown the significant effects of the tide, which differs from the results observed for the present specific case study.
The thermal plume dispersion is highly influenced by the magnitude and direction of the currents.Figure 13 illustrates the evolution of the predicted thermal field for the same time intervals used in Figure 12.The four thermal discharges of the power plant can be clearly observed with the aid of the isotherms.Because of the predominant flow direction toward the southeast, the thermal plumes are dispersed in the same general direction.Moreover, the plumes are dispersed along the coast, at least during the cold front, for the first 24 h of the simulation time period (see Figure 13a-c).When the cold front weakens at the end of the simulation time period, a seaward thermal dispersion is observed.It is evident that the influence area of the thermal plumes of Discharges 1 and 3 is larger than the area covered by the thermal plumes of Discharges 2 and 4.This is likely because the discharge rate at Discharge 1 and 3 is significantly greater than for Discharge 2 and 4. The directionality developed by the thermal plume dispersion is undesirable from the point of view of the efficiency of the power plant, because large amounts of hot water emitted from Discharges 1 and 3 could recirculate toward the intake.As mentioned above, the efficiency of the power plant is compromised while hot water is taken for cooling.Table 5 shows the calculated influence areas per range of temperature at time intervals of 12 h.The temperature differences refer to the temperature of Discharge 1, 34.89 °C, which is the highest.The percentages of heat dissipated are calculated considering the maximum temperature difference of 6.23 °C, which corresponds to the difference between the temperature of Discharge 1 and the average of the thermistor measurements of 28.66 °C.The isotherm 29.06 °C covers an area of almost 12 km 2 during the first 12 h (13 October at 12:00), which is the area of maximum influence of the simulation time period.With this temperature value, the plume has dissipated 94% of heat.After the first 24 h of the simulation time period (14 October at 00:00), the influence area of the isotherm 29.06 °C decreases to reach approximately 0.82 km 2 .This occurs when the current velocities are more intense over the whole study domain, as Figure 12 shows for 14 October at 00:00.After 36 h of the simulation time period (14 October at 12:00), the influence area of the isotherm 29.06 °C covers a surface of approximately 2.42 km 2 .At this particular time, the calculated influence area agrees well with that estimated at the sea surface with the aid of the CTD measurements the same day (see Section 3.1), which is approximately 2.84 km 2 .At all time intervals, the isotherm 29.72 °C, which corresponds to 83% of heat dissipated, significantly decreases the influence area with respect to the area covered by the isotherm 29.06 °C.The influence areas of the rest of the isotherms gradually decrease as the temperature increases to reach low heat dissipation regions, covering relatively small influence areas.Table 5 shows the calculated influence areas per range of temperature at time intervals of 12 h.The temperature differences refer to the temperature of Discharge 1, 34.89 • C, which is the highest.The percentages of heat dissipated are calculated considering the maximum temperature difference of 6.23 • C, which corresponds to the difference between the temperature of Discharge 1 and the average of the thermistor measurements of 28.66 • C. The isotherm 29.06 • C covers an area of almost 12 km 2 during the first 12 h (13 October at 12:00), which is the area of maximum influence of the simulation time period.With this temperature value, the plume has dissipated 94% of heat.After the first 24 h of the simulation time period (14 October at 00:00), the influence area of the isotherm 29.06 • C decreases to reach approximately 0.82 km 2 .This occurs when the current velocities are more intense over the whole study domain, as Figure 12 shows for 14 October at 00:00.After 36 h of the simulation time period (14 October at 12:00), the influence area of the isotherm 29.06 • C covers a surface of approximately 2.42 km 2 .At this particular time, the calculated influence area agrees well with that estimated at the sea surface with the aid of the CTD measurements the same day (see Section 3.1), which is approximately 2.84 km 2 .At all time intervals, the isotherm 29.72 • C, which corresponds to 83% of heat dissipated, significantly decreases the influence area with respect to the area covered by the isotherm 29.06 • C. The influence areas of the rest of the isotherms gradually decrease as the temperature increases to reach low heat dissipation regions, covering relatively small influence areas.
In order to observe the decay length of the thermal plumes for Discharges 1 and 3, the decay curves of temperature along the axis (seaward) were graphed at time intervals of 6 h (see Figure 14).For the case of Discharge 1 (see Figure 14a), 90% of the temperature is dissipated within about 1500 m from the outfall for all time intervals.On the other hand, the temperature of Discharge 3 (see Figure 14b) is dissipated in a shorter distance, 90% within about 750 m.Thus, the plume effects are local, which is confirmed by the CTD measurements of Figure 8.In order to observe the decay length of the thermal plumes for Discharges 1 and 3, the decay curves of temperature along the axis (seaward) were graphed at time intervals of 6 h (see Figure 14).For the case of Discharge 1 (see Figure 14a), 90% of the temperature is dissipated within about 1500 m from the outfall for all time intervals.On the other hand, the temperature of Discharge 3 (see Figure 14b) is dissipated in a shorter distance, 90% within about 750 m.Thus, the plume effects are local, which is confirmed by the CTD measurements of Figure 8.The temperature history was recorded at the intake, just between the two breakwaters of the plant, to observe thermal recirculation to the power plant (see Figure 15).Compared against 28.66 °C, the average of the thermistor measurements during the simulation time period, the maximum temperature increase observed during the simulation time period is about 1.1 °C.This increase is about 17.66% of 6.23 °C, the maximum temperature difference defined above; thus, it can be considered as not significant.Given the directionality shown by the thermal plumes dispersion, the position of the breakwater north of the plant protects the intake from recirculating thermal flow coming from Discharges 1, 2, and 3. Thus, the plume dispersion that occurred during the cold front passage did not significantly affect the efficiency of the cooling systems of the power plant, regarding thermal plume recirculating.

Conclusions
The thermal plume dispersion from the CTPALM discharges has been successfully characterized by numerical modeling.In order to analyze the relevant features of the study zone and to implement the numerical model, a field investigation was also performed.The scenario of a cold front, particularly studied here, revealed that the strong winds occurring during the passage of this phenomenon are capable of causing significant changes in the local coastal hydrodynamics and in thermal plume dispersion as a consequence.Under this scenario, northward winds along the coast result in thermal plume dispersion southward parallel to the shoreline where the power plant intake The temperature history was recorded at the intake, just between the two breakwaters of the plant, to observe thermal recirculation to the power plant (see Figure 15).Compared against 28.66 • C, the average of the thermistor measurements during the simulation time period, the maximum temperature increase observed during the simulation time period is about 1.1 • C.This increase is about 17.66% of 6.23 • C, the maximum temperature difference defined above; thus, it can be considered as not significant.Given the directionality shown by the thermal plumes dispersion, the position of the breakwater north of the plant protects the intake from recirculating thermal flow coming from Discharges 1, 2, and 3. Thus, the plume dispersion that occurred during the cold front passage did not significantly affect the efficiency of the cooling systems of the power plant, regarding thermal plume recirculating.
Water 2016, 8, 482 14 of 16 In order to observe the decay length of the thermal plumes for Discharges 1 and 3, the decay curves of temperature along the axis (seaward) were graphed at time intervals of 6 h (see Figure 14).For the case of Discharge 1 (see Figure 14a), 90% of the temperature is dissipated within about 1500 m from the outfall for all time intervals.On the other hand, the temperature of Discharge 3 (see Figure 14b) is dissipated in a shorter distance, 90% within about 750 m.Thus, the plume effects are local, which is confirmed by the CTD measurements of Figure 8.The temperature history was recorded at the intake, just between the two breakwaters of the plant, to observe thermal recirculation to the power plant (see Figure 15).Compared against 28.66 °C, the average of the thermistor measurements during the simulation time period, the maximum temperature increase observed during the simulation time period is about 1.1 °C.This increase is about 17.66% of 6.23 °C, the maximum temperature difference defined above; thus, it can be considered as not significant.Given the directionality shown by the thermal plumes dispersion, the position of the breakwater north of the plant protects the intake from recirculating thermal flow coming from Discharges 1, 2, and 3. Thus, the plume dispersion that occurred during the cold front passage did not significantly affect the efficiency of the cooling systems of the power plant, regarding thermal plume recirculating.

Conclusions
The thermal plume dispersion from the CTPALM discharges has been successfully characterized by numerical modeling.In order to analyze the relevant features of the study zone and to implement the numerical model, a field investigation was also performed.The scenario of a cold front, particularly studied here, revealed that the strong winds occurring during the passage of this phenomenon are capable of causing significant changes in the local coastal hydrodynamics and in thermal plume dispersion as a consequence.Under this scenario, northward winds along the coast result in thermal plume dispersion southward parallel to the shoreline where the power plant intake

Conclusions
The thermal plume dispersion from the CTPALM discharges has been successfully characterized by numerical modeling.In order to analyze the relevant features of the study zone and to implement the numerical model, a field investigation was also performed.The scenario of a cold front, particularly studied here, revealed that the strong winds occurring during the passage of this phenomenon are capable of causing significant changes in the local coastal hydrodynamics and in thermal plume dispersion as a consequence.Under this scenario, northward winds along the coast result in thermal plume dispersion southward parallel to the shoreline where the power plant intake is located; however, the north breakwater protects it from recirculating hot water.The main hot water source of the power plant comes from Discharges 1 and 3, with a contribution of 93% of the total discharge.Results also showed that the influence area of the thermal plume is at most 12 km 2 , dissipating almost 94% of heat within the first 1500 m seaward.
Finally, in order to implement an effective monitoring program of the thermal plume and to establish appropriate prevention and mitigation measures in terms of the environmental effects, this study has to be reproduced under different physical conditions to identify the features of other typical seasonal events.Data compilation and additional field measurements are necessary to complement numerical simulations and gather more reliable results.

Figure 1 .
Figure 1.Location of the CTPALM power plant.

Figure 1 .
Figure 1.Location of the CTPALM power plant.

Figure 2 .
Figure 2. Study domain and location of the ADCP and thermistor moorings.

Figure 2 .
Figure 2. Study domain and location of the ADCP and thermistor moorings.

Figure 3 .
Figure 3. Bathymetry of the study zone.

Figure 4 .
Figure 4. Location of the four traces performed using boat-based CTD on 14 October 2011.

Figure 3 .
Figure 3. Bathymetry of the study zone.

Figure 3 .
Figure 3. Bathymetry of the study zone.

Figure 4 .
Figure 4. Location of the four traces performed using boat-based CTD on 14 October 2011.

Figure 4 .
Figure 4. Location of the four traces performed using boat-based CTD on 14 October 2011.

Figure 5 .
Figure 5. Tidal variation and current velocity measured with ADCP.

Figure 5 .
Figure 5. Tidal variation and current velocity measured with ADCP.

Figure 6 .
Figure 6.Correlation analysis between the surface and bottom currents, which shows a match, characteristic of a barotropic structure.

Figure 7 .
Figure 7. Frequency spectrum of the surface current parallel to the coastline.

Figure 8 .
Figure 8.Time series of seawater temperature measured in front of the CTPALM.(a) Mooring TM 1; (b) Mooring TM 2.

Figure 6 .
Figure 6.Correlation analysis between the surface and bottom currents, which shows a match, characteristic of a barotropic structure.

Figure 6 .
Figure 6.Correlation analysis between the surface and bottom currents, which shows a match, characteristic of a barotropic structure.

Figure 7 .
Figure 7. Frequency spectrum of the surface current parallel to the coastline.

Figure 8 .
Figure 8.Time series of seawater temperature measured in front of the CTPALM.(a) Mooring TM 1; (b) Mooring TM 2.

Figure 7 .
Figure 7. Frequency spectrum of the surface current parallel to the coastline.

Figure 6 .
Figure 6.Correlation analysis between the surface and bottom currents, which shows a match, characteristic of a barotropic structure.

Figure 7 .
Figure 7. Frequency spectrum of the surface current parallel to the coastline.

Figure 8 .
Figure 8.Time series of seawater temperature measured in front of the CTPALM.(a) Mooring TM 1; (b) Mooring TM 2.

Figure 8 .
Figure 8.Time series of seawater temperature measured in front of the CTPALM.(a) Mooring TM 1; (b) Mooring TM 2.

Figure 9 .
Figure 9. Sea surface CTD temperature measurements in front of the CTPALM.Solid lines show the location of the discharges.(a) First trace; (b) Second trace; (c) Third trace; (d) Fourth trace.
Numerical simulations were performed to calibrate and validate the numerical model.The calibration of the dynamic field consisted of a trial wind drag coefficient d C and Manning coefficient n , adjusted in subsequent simulations until the pattern of the predicted depth-averaged horizontal velocities best matched that measured with ADCP over the simulation time period.There

Figure 9 .
Figure 9. Sea surface CTD temperature measurements in front of the CTPALM.Solid lines show the location of the discharges.(a) First trace; (b) Second trace; (c) Third trace; (d) Fourth trace.

Figure 10 .
Figure 10.Measured and predicted depth-averaged horizontal velocities of the calibration tests.(a)U -velocity; (b) V -velocity.

Figure 10 .
Figure 10.Measured and predicted depth-averaged horizontal velocities of the calibration tests.(a) U-velocity; (b) V-velocity.

Figure 15 .
Figure 15.Temperature signal at the intake of the CTPALM.

Figure 15 .
Figure 15.Temperature signal at the intake of the CTPALM.

Figure 15 .
Figure 15.Temperature signal at the intake of the CTPALM.

Table 1 .
Initial conditions and physical parameters used for the numerical simulations.

Table 2 .
Wind drag and Manning coefficients tested during calibration.

Table 2 .
Wind drag and Manning coefficients tested during calibration.

Table 3 .
Calculated coefficients for each calibration test.

Table 3 .
Calculated coefficients for each calibration test.

Table 4 .
Calculated efficiency coefficients for each thermistor.

Table 4 .
Calculated efficiency coefficients for each thermistor.

Table 5 .
Calculated influence areas of the thermal plumes per range of temperature.

Table 5 .
Calculated influence areas of the thermal plumes per range of temperature.