Thermal Recirculation Modeling for Power Plants in an Estuarine Environment

Many power plants require large quantities of water for cooling purposes. The water taken from the source water body (e.g., lakes, estuaries, bays and rivers) circulates through the plant and returns to the source through outfall with a higher temperature. For optimal performance of the power plant, the intake inlet and discharge outlet should be meticulously placed so that the heated water will not recirculate back into the power plant. In this study, the Flow module of the Delft3D software is employed to simulate the temperature transport within the study area in three-dimensional and nested format. Model results are used to optimize the location of intake inlets, outfall outlets and diffuser port orientations. The physical processes used in the study are tidal fluctuations, winds, river discharges, salinity and temperature. The subject power plant (power plant parameters presented in this paper are realistic; however, they do not target any specific power plant within the study area) has a nominal capacity of 2600 MW and is planned to be located in Delaware Bay, USA. Existing field measurements are used to calibrate the model in a coupled two-staged fashion for main tidal constituents, currents and water temperature. The sensitivity of the model against various input parameters is tested, and conservative values are selected. The location of the intake is fixed, and the location of the outfall is changed until the thermal impact to the intake is less than 1 ◦ C. Analysis of the results shows that there is a linear logarithmic relation between the excess temperatures at the intake inlet and horizontal eddy diffusivity. The k − e turbulence closure results in higher excess temperature and a more conservative design. Extending the outfall location to the deeper portion of the estuary combined with port orientations reduces the impact by keeping the thermal plume away from the intake inlet and meeting the established criteria. It is concluded that an approximate distance of 1300 m is the optimal location for the power plant outfall outlets. In addition, the diffuser ports should not discharge the heated water toward the intake and have to be oriented away from the line connecting outfall to the intake.


Introduction
A once-through cooling power plant draws hundreds of millions, and in some cases billions, of gallons of water each day from nearby lakes, rivers or oceans [1].The withdrawn water is boiled and turns into steam, which spins the turbines in order to generate electricity.In the power industry, intake and discharge systems are known as the bloodlines of every power plant, as they deliver cold water to the plant and discharge the hot water back to the source.Their locations should be analyzed meticulously to reduce the rise in water temperature at the intake to a value established by the project (e.g., less than 1 • C).If the intake is placed improperly with respect to the plant discharge, or to the discharge from another plant, or if sea currents are very strong, a rise in temperature might occur.This will negatively affect the plant production and should be avoided [2,3].Furthermore, environmental agencies establish a regulatory mixing zone where the excess temperature is limited to permissible values.Thermal recirculation studies are frequently used to demonstrate that the thermal limits are met.However, environmental applications were not the focus of this study and, therefore, are not discussed.
The heated water released from the plant outlet transports and spreads (i.e., disperses) in the source water body by a combination of all available hydrodynamic forcing.Understanding the shape, size and timing of the thermal plume is a key factor in choosing the right location for any power plant intake.Dye and drogue studies are traditional and relatively inexpensive field methods for estimating thermal plume characteristics [4].However, these types of studies are mainly applicable to existing installations or the expansion of facilities.In addition, they do not cover all of the possible ambient conditions, as they are limited to the time of the field study.On the other hand, numerical model simulations are valuable in evaluating the impact of a wide range of ambient conditions at reasonable cost and in a manageable time frame.
Engineers, researchers and scientists have been continuously developing new methodologies and techniques to perform thermal plume studies using a variety of models in different site settings.The majority of such studies is focused on power plants located in coastal areas and not necessarily estuaries [5][6][7][8][9][10][11].The work in [5] studied the outfall excess temperature for a 6000-MW nuclear power plant in coastal nearshore waters.They used a coarse grid modified finite element version of the hydrodynamic model (CAFÉ) developed by [12].Different outfall settings including the construction of short and long jetties were tested.The excess temperature at the intake inlet ranged from 2.8 • C (long jetty) to 3.5 • C (short jetty).Considering the minimal difference (only 0.7 • C) between the two scenarios, the short jetty alternative was selected that significantly reduced project cost.
The work in [6] applied a 3D finite difference model (ASL-COCIRM) to study the recirculation of cooling waters in the Burrard Generating Station (BGS), Canada.They demonstrated that the model is capable of capturing heat recirculation into the intake with a very high grid resolution (2.5 m × 2.5 m).
The work in [7] used Delft3D-FLOW in regional format (i.e., large scale) to study the thermal recirculation between the intake and outfall of a nuclear power plant in a coastal setting.They were able to show that the excess temperate at the intake inlet is less than 0.5 • C under severe condition.Tide, current and wind were the driving forces in all of the above studies.River discharge has not been considered in any of the analysis.
The work in [13] studied effluent aspects of the Liquid Natural Gas (LNG) plant outfalls (similar to thermal recirculation) in a coastal ecosystem.Their goal was to design the plant outfall diffusers by meeting the environmental mixing zone criteria for various chemical effluents.They used the combination of CORMIX and Delft3D-FLOW to capture both near-field and far-field aspects of the mixing processes.In the near-field, mixing is dominated by the discharge condition, while in the far-field, the mixing is dominated by ambient conditions.They suggested that the tidal range and surface wind govern the diffuser design.In all of the above studies, considering the application, the combination of near-field and far-field or far-field-alone models are investigated.
As will be described in the following section, Delaware Estuary is one of the most important water bodies on the east coast of the United States.Numerous studies have been performed for this water body to address tidal/sub-tidal variations, salinity intrusion, buoyant flow and storm surge [14][15][16][17][18][19].However, the literature review did not result in any study that addresses thermal recirculation within this water body.
The objective of this study is to estimate the thermal plume behavior using the Delft3D-FLOW model software under various ambient conditions.The results of this study would be used to recommend the most economical location for the plant discharge system (i.e., outfall) relative to the intake inlet.This study is unique as it provides a new cost-effective and scientific technique that combines two-Dimensional (2D) and three-Dimensional (3D) models in a complex estuarine environment.This methodology can be used to help Engineering Procurement and Construction (EPC) projects to locate the power plant's intake and outfall systems with minimal information at the early stages of the projects.

Study Area
The subject power plant is located in Delaware Estuary (Estuary) at approximately 10 km south of its confluence with the Chesapeake and Delaware (C & D) canal.The Estuary is 210 km long and comprises tidal reaches of the Delaware River (starting in Trenton, New Jersey) and all of the Delaware Bay [18].The width of the estuary increases toward the mouth to 18 km, which extends between Cape May and Cape Henlopen.The maximum width of the estuary is 45 km (see Figure 1) at about 11 km north of Cape May.The bathymetry of the estuary varies significantly with an average value of 7 m [20].A deep channel (exceeding 28 m) runs through the estuary that connects the Atlantic Ocean to the ports and facilities located upstream of the area.The Delaware is the longest undammed river in the United States east of the Mississippi, extending 530 km from the confluence of its east and west branches at Hancock, N.Y., to the mouth of the Estuary, where it meets the Atlantic Ocean.The hydrodynamics of the Estuary are complex, as they are influenced by various forcing, such as wind, wind-generated waves, swells, river and tributary discharges, Atlantic Ocean tides and Chesapeake Bay tides through the C & D canal.The estuary belongs to a weakly stratified or well-mixed class of estuaries [19].The world's largest population of horseshoe crabs spawns in Delaware Estuary.The Delaware Estuary was designated as a Ramsar site in May of 1992.Furthermore, the estuary is one of the most important navigational channels in the United States, its second busiest waterway after the Mississippi River [21].

Plant Cooling System
A once-through cooling system is considered for the proposed 2600-MW power plant with submerged intakes and outlets.The intake heads are located at a water depth of approximately 7 m.After circulating the water through the plant, the heated cooling water downstream of the condensers would return to a seal well before entering the discharge pipes.For this study, a multi-port diffuser with 20 risers and 40 ports each 1 m in diameter is considered.The goal is to discharge the water as close to the shore as possible to minimize the required pipe material.This study evaluates various options to propose the most economically-feasible discharge point, as well as the orientation of outfall ports.It is common to consider 0.05 m 3 /s of cooling water for each 1-MW capacity of the power plant.Therefore, the proposed power plant will withdraw approximately 130 m 3 /s (i.e., 2600 MW × 0.05 m 3 /(s•MW)).Intake heads are considered to be 1 m below the Mean Low Low Eater (MLLW) level.The MLLW is approximately 0.93 m below the Mean Sea Level (MSL), which puts the inlet head 1.93 m below the MSL.The outlet of the discharge pipe will be placed 1 m above the bottom of the estuary.

Model Configurations
The hydrodynamic and transport simulations are performed with the open source three-Dimensional (3D) Delft3D-FLOW (Version 3.56.29165,Deltares, The Netherlands).This model calculates non-steady flow phenomena under various tidal and meteorological forcing on a curvilinear horizontal grid.The sigma coordinate system is used in the vertical to resolve the complex bathymetry of the region.Three-dimensional non-linear shallow water equations derived from 3D Navier-Stokes equations for an incompressible fluid, under Boussinesq assumptions along with continuity equations, are solved by an Alternating Direction Implicit (ADI)-type of formulation for the barotropic pressure [22].Both implicit (central difference) and explicit upwind (finite volume) schemes can be implemented [23].Intake and discharge are modeled as sink and source terms with optional momentum or normal factors.If the momentum option is included, the momentum of the discharge is taken into account by specifying the discharge velocity and direction.The normal-type discharge does not consider any specific aspects of the outfall jet.Both types of discharges are tested in this study to better understand their impacts.All governing partial differential equations and their associated numerical techniques are given in the Delft3D-FLOW manual [23].

Model Domain
The simulation domain is divided into overall and nested domains (see Figure 1).The overall domain is set up in a two-Dimensional (2D) fashion to generate the hydrodynamic boundary condition for its nested domain that is run in 3D format.The numerical grids are curvilinear and generated using the RGFGRID module of Delft3D.Primary grid properties, such as orthogonality, edge length, resolution, smoothness and curvature, are tested against standard values and modified as deemed necessary.

Overall
The overall domain includes 130 km of the estuary extending from the mouth to the Commodore Barry Bridge on the Delaware River (Figure 1).The offshore portion of the domain extends 90 km toward the Atlantic Ocean with a width of approximately 195 km.The grid contains 10,400 cells covering an area of 17,000 km 2 .The highest horizontal resolution (14-40 m) is required in the C & D canal.Larger grid cells (up to 5.5 km wide) are considered far offshore (Figure 2).The overall grid has a resolution of 400 m along the estuary and 200 m perpendicular to the shore near the power plant region.No vertical cell is considered for the overall grid, as it is run only in 2D format.

Nested
The nested domain starts 5.5 km south of the C & D canal confluence with the estuary and extends another 8 km south (Figure 1).The grid contains 2050 cells covering an area of 28.5 km 2 .Two sets of cell resolutions, identified as Mesh 1 and Mesh 2 in Figure 3, are tested.In Mesh 1, cell resolution varies from 20 m in the vicinity of the intake and outfall to 190 m away from them.In Mesh 2, cell resolution is set to 80 m near the intake or outfall and 190 m farther away from them.Sensitivity runs show that both meshes are able to generate almost identical results.Since Mesh 2 is less computationally demanding, it is used for all other model runs.Vertical resolution is provided by 10 sigma levels; spacing is 5% of the water column near the surface and bottom (upper and lower two levels) with a linear increase toward the interior.This vertical spacing allows reasonable model run time, while still resolving features of interest.
Bathymetry of the overall and nested domains is presented in Figures 2 and 3, respectively.Coastline and bathymetry (on the shelf and within the estuary) are extracted from the ADvanced CIRCulation (ADCIRC) Federal Emergency Management Agency (FEMA) Region III storm surge model [24]

Standard Setting
Delft3D uses the standard Alternating Direction Implicit (ADI) for solving the shallow water equations [23].Sensitivity tests are carried out in order to determine the largest time step for which the ADI-method still yields accurate results.Temporal differencing for the simulation is set at three minutes for the overall grid and one minute for the nested grid to avoid numerical instabilities.
The 3D turbulence closure is set to the second order k − formulation [25,26].Models are run with no background vertical viscosity or diffusivity.Since the salinity process is not considered in the analysis, a constant background salinity of 11 parts per thousand (ppt) is selected for the model grids.This average salinity is based on the analysis of salinity observations from National Oceanic and Atmospheric Administration (NOAA) stations in the vicinity of the nested grid.Delft3D uses the UNESCO equation for the state of seawater (1980) as the relation between salinity, temperature and density.The salinity of the discharged water is kept the same as that of the withdrawn water from the intake.A free slip condition is applied along land boundaries, setting zero tangential shear stress at the coast.This is a valid and safe approach for large-scale calculations similar to this study [23].Changing this boundary condition from free slip to no slip had little to no impact on the thermal plume spreading and magnitude, as both intake and discharge heads are located a safe distance from shore and from the bottom of the estuary.Both normal and momentum-type discharges are considered for establishing model scenarios.For the momentum type, the velocity of the outfall discharge is set to 4.14 m/s based on the size and number of ports.

Boundary Conditions
The wind process is considered for both the overall and nested models, while the temperature process is applied only to the nested model.Conservatively, heat exchange through the free surface is set to zero.For the overall model, the offshore boundaries, as well as the C & D boundary (Figure 2) are fixed to the tides by specifying tidal amplitudes and phases.The nested grid takes its hydrodynamic boundary conditions at the north and south (Figure 3) of its domain from the overall grid.The temperature boundary conditions are applied manually to the nested grid.Delaware River discharge is imposed as a vertically uniform input at Commodore Barry Bridge for the overall grid (Figure 2).

Model Forcing
The three main forcing agents included in the study are: tides, Delaware River discharge and winds.Effort is made to capture the main features and processes of the Delaware coastal currents in the vicinity of the power plant through proper application of these forces.Waves are not included in the hydrodynamic modeling.Even though waves are important in turbulent mixing processes, they do not contribute to the net fluid transport.Sensitivity analysis shows that a seven-day spin-up period is sufficient to achieve a quasi-steady state.Each of the forcings is described below.

Tides
As tidal currents are dominant in this region, they must be included.The dominant tidal constituent is the principal lunar semi-diurnal (M2) with a very large tidal volume to mean freshwater discharge.As a result, the stratification is weak, ranging from zero to a few parts per thousand (ppt) [16].In addition to M2, eight other constituents-S2, N2, K2, K1, O1, Q1, M4 and M6-are included in the tidal boundary condition.The Form Factor (FF) defined as (K1 + O1)/(M2 + S2) is close to 0.25, an indication of a semi-diurnal tidal condition.Furthermore, the estuary is classified as mesotidal with tides ranging between 2 m to 4 m.The combined effects of these tidal constituents generate variable tidal current along the estuary.Figure 4 shows the tide time series, at the middle of the nested domain, generated by the model after applying tides at north and south boundaries.The water level fluctuates between −1.1 and 1.2 m MSL.Delft3D uses the following equation for combining all of the constituents and constructing the variable water level.
where H(t) is predicted water level at time t, k the number of relevant constituents, i the index of a constituent, A i and φ i the local tidal amplitude and phase of a constituent, F i the nodal amplitude factor, T i the period and V i the astronomical argument.The values for A i and φ i for the selected constituents are input variables.Delft3D computes V i , T i and F i for each constituent (for the period of prediction).The subtidal fluctuations are primarily driven by a direct coupling to the continental shelf and an indirect coupling to the Chesapeake Bay through the C & D canal [15].Tidal heights are imposed along all open boundaries by specifying amplitude and phases.In general, four open tidal boundaries are considered, named as south, west, east and CnD , in Figure 2. The south boundary is further divided into three segments of west, middle and east.In Delft3D, the amplitude and phase at the beginning and end of each open boundary are required as input.These parameters are all extracted from the Western North Atlantic, Caribbean and Gulf of Mexico Tidal Databases [27] using ADCIRC.Table 1 presents the amplitude for each tidal constituent and boundary.

River Discharges
Approximately 70 percent of the freshwater inflow enters the estuary by the Delaware and Schuylkill Rivers [16].In addition to these two sources, Neshaminy and Ranconas Creeks are considered in this study.The U.S. Geological Survey (USGS) instantaneous river discharge records at these sources (see Table 2) for the two-month simulation period are summed and used as the discharge boundary condition.The total freshwater inflow is shown in Figure 5.The imposed freshwater is representative of an average flow condition and varies between 200 m 3 /s and 1200 m 3 /s.Table 2 shows the contributing percentages of each source.Statistical discharges shown in Figure 5 are estimated applying the method of moments on instantaneous observed flows from 1981-2007.Almost 88 percent of the flows in the two months of simulation (red line in Figure 5) are smaller than the flows with a 10-year return period.The average of the simulated flow is 420 m 3 /s, which is very close to the long-term (1984-2007) average flow of 450 m 3 /s.The work in [16] states that flood discharges impact mainly upstream parts of the estuary, which is farther from the study area that is the focus of this analysis.In addition, considering the location of the outlet with respect to the intake (Figure 3), higher river discharges produce a force against the tidal prism, which moves the thermal plume further away from the intake.Therefore, the application of an average river discharge can be an appropriate representation of the average forcing throughout the power plant life cycle.There are other tributaries flowing into the estuary between the mouth and upstream boundary, but not included, as their effects are considered localized [18] and not significant to the goals of this study.

Winds
For calibration runs, time-variable spatially-uniform winds are applied to the entire domain.Hourly surface meteorological observations for the NOAA Buoy 44009 (46 km southeast of Cape Henlopen) during July and August of 2003 are obtained from the National Data Buoy Center (NDBC).Figure 6 depicts the wind rose diagram and time series of the wind for the calibration simulation.The average wind has a speed of 5.39 m/s with a dominant south-northeast direction.For the actual performance runs (i.e., model scenarios), conducted to study the plume behavior, 10 years of historical wind data reported at Buoy 44009 are used.The theory-of-random waves is applied to the historical wind data to establish a typical and frequent rough weather condition.In this approach, the 1/3-highest wind speeds for summer months (April-November) and winter months (December-March) are averaged.The resulting wind values are 5.5 m/s for summer and 7.6 m/s for winter.The dominant direction is established using directional frequency distribution analysis.The along-shelf wind direction blowing from the power plant outlet the toward intake (southwest-northeast) is dominant, conservative and applied to all model scenarios.

Model Calibration and Sensitivity
The overall model is calibrated for water level and current, while the nested model is calibrated for temperature.A two-month period from 1 July 2003-1 September 2003 is selected as the calibration time frame.The first seven days of the analysis are considered as the spin-up period.

Water Level and Current
The hydrodynamic is calibrated for water level and current by comparing the model results against observed values.Water level observations and currents are extracted from NOAA stations (see Figure 1) and field studies [14,19].The bottom roughness can be parametrized by several formulas, including Manning, White-Colebrook and Chezy.The White-9Colebrook formulation, as implemented on the Delft3D-Flow using Equation ( 2) with a spatially-uniform and constant Nikuradse roughness length of 0.05, is applied to best match model results with observed tidal characteristics within the estuary [28].
where C 2D is the friction coefficient, H is the total depth of water column and k s is the Nikuradse roughness length.As stated, the final goal of this study is to evaluate the impact of power plant thermal outfall on the intake system, which is located within the nested grid.Therefore, the calibration is aimed at providing reasonable hydrodynamics for the nested grid by calibrating the overall grid for stations located at both sides of the nested grid.The Ship John Shoal and Reedy Point NOAA stations are located south and north of the nested grid (see Figure 1), respectively, and used for the hydrodynamic calibration.In addition to the visual comparison, to obtain an objective and quantitative measure of how well the calibrated model results compare to the observed data, the statistical parameter of Nash-Sutcliffe model Efficiency coefficient (NSE) is calculated.It is defined as: where Y obs i is the i-th observation for the constituent being evaluated, Y sim i is the i-th simulated value for the constituent being evaluated, Y mean is the mean of the observed data for the constituent being evaluated and n is the total number of observations.NSE ranges between negative infinity and 1.0 (one inclusive), with NSE = 1 being the optimal value.Model outputs are extracted every six minutes to synchronize them with the time stamps of the observations.The calibrated model resulted in an NSE value of 0.92 for the water surface levels, which is an indication of acceptable model calibration.Figure 7 presents a plot of the modeled and measured water level at the Ship John Shoal NOAA station with associated quality indices, such as bias, correlation and NSE.
Various current measurements are performed within the study area and presented in [14,19].Field observations show that the tidal current has a magnitude of approximately 1.0 m/s near the surface and 0.85 m/s near the bottom during both flood and ebb phases.River discharge during field measurements was 683 m 3 /s, which is within the range of the flows considered in this analysis.In the calibrated model, the velocities are synched with the tidal variations with a maximum value of 1.02 m/s at the surface and 0.81 m/s near the bottom.Modeled velocity results showed the same variability and magnitude as those presented in [19].A comparison between total cross-sectional flow between simulated and observed values was performed, and satisfactory results were achieved.Common to estuaries, the maximum velocity occurs at the beginning of the flood limb of the tide, which is supported by current analysis.

Temperature
The nested model takes its hydrodynamic boundary condition at its north and south boundaries (Figure 3) from the calibrated overall model and is further calibrated for the temperature.The horizontal eddy diffusivity is changed within an acceptable range until the model results are matched with the observed values.With lack of water temperature measurements in the vicinity of the nested model, a temperature gradient for each month of simulation is established using the field data extracted from NOAA stations at Ship John Shoal, New Jersey, and Brandywine, Delaware (Figure 1).
where ∆T is the temperature gradient in degree Celsius and D is the distance between two points along the axis of the estuary in kilometers.It should be noted that the temperature has a positive gradient toward the upper parts of the estuary; meaning that the water temperature is higher at Ship John Shoal compared with Brandywine.The synthetic temperature boundary conditions at the north and south side of the nested grid are established by adding their associated gradients to the temperate measurements at Ship John Shoal for each month of simulation.The south and north of the nested grid are 20.8 km and 28.5 km away from the Ship John Shoal. Figure 8 presents the surface temperature for the Ship John Shoal observations, synthetic south and north boundaries and model results at a point in the middle of the nested domain.As expected, the model results fall between upstream and downstream boundaries.This level of calibration is assumed satisfactory for the purpose of modeling since the main results are related to impact modeling and a relative comparison of two situations before and after the introduction of the cooling water system.Originally, a horizontal eddy diffusivity of 10 m 2 /s resulted in good agreement between modeled and observed values.However, in addition to calibrating the overall and nested models for the most important processes, sensitivity analyses are performed to help select the most conservative parameter.Horizontal and vertical eddy diffusivities, wind direction and turbulence closures were the sensitivity parameters.Analyses show that the excess temperature at the intake inlet has an inverse linear relation with the horizontal eddy diffusivity for values higher than one and no change for lower values.Therefore, a horizontal eddy diffusivity of one and lower can produce the most conservative results.Increasing the vertical eddy diffusivity from base value of zero resulted in decreasing the excess temperature, and therefore, it was kept as the default value of zero.Between the three tested turbulence closures of k − , k − l, and algebraic, the k − resulted in the highest excess temperature, and therefore, it is used in all of the analyses.As expected, winds blowing from southwest toward the northeast resulted in the highest excess temperature at the intake inlet, and this is used in all of the model scenarios, except the calibration scenarios.

Model Scenarios
Upon accepting the calibration results, the nested model is used to further evaluate the power plant withdrawal and discharge impacts.Intake withdrawal and outfall discharges are included as localized internal sinks and sources.Model scenarios are established, including the most conservative ambient condition and model parameters to capture the expected severe, yet reasonable, condition.Though the region is prone to hurricanes, these events are not considered in this analysis as their impacts on thermal plume are considered short-term with a low probability of occurrence, as hurricane events last for only a few days.The general and specific characteristics of the model scenarios are presented in Tables 3 and 4. Software default parameters are used for those not described here.A full power plant capacity of 2600 MW is considered in all models, with the total intake withdrawal and corresponding outfall discharge of 130 m 3 /s.As depicted in Figure 3, two locations for the heated discharge are tested; one located 1300 m and the other 800 m from the shoreline.Power plant system descriptions show that the corresponding temperature rise for the water passing through the plant condenser is about 10 • C. Wind speed varies from season to season with average summer winds of 5.5 m/s and winter winds of 7.6 m/s.Winds are applied blowing from southwest toward the northeast, which generates shear forces on the water surface from outfall toward intake.Model scenarios are run for 56 days from 7 July 2003-1 September 2003.The combination of season parameters with outfall location, type and direction results in establishing a minimum of 20 model scenarios, as outlined in Table 4.

Results and Discussions
For the purpose of better understanding the mixing processes and their dependability on various factors, 20 simulations were conducted.Scenarios present situations with various wind speeds, water temperatures, outfall locations, outfall types and port orientations.The surface temperatures during a tidal cycle on 21 August 2013 and for Scenario 3 are given in Figure 9. Analysis of the excess temperature time series at the intake location shows that one of the maximums occurs on this date.This can be associated with low river discharge (less than average flow; see Figure 5), which generates less resistance against the flow created by the tidal prism.In Scenario 3, the summer season is simulated with the outfall and the intake located almost 800 m and 500 m from the shoreline (see Table 4), respectively.Tidal currents are landward during the flood and seaward during the ebb.As expected, the thermal plume spreading direction changes from northward during flood tide to southward during ebb tide.Results from Scenario 3 show that the worst-case condition for the intake occurs at the beginning of flood tide (identified by Number 2 in Figure 9) and ebb tide (identified by Number 5 in Figure 9).This is reasonable as the highest current velocity and shear stress is known to occur during the first flood and last ebb in an estuarine setting [29,30].Surface temperature during the same period and season is shown in Figure 10 for a condition with outfall location almost 400 m farther from its location in Scenario 3. Comparing the results in Figures 9 and 10 shows that excess temperature at the intake is reduced for Scenario 1. Figure 11 graphs the maximum excess temperature for all scenarios.Combining values shown in the figure with parameters associated with each model scenario (see Table 4) reveals that the excess temperature criteria (i.e., less than 1 • C) might only be achieved by a combination of longer outfall pipe and proper port orientation.For example, the outfall in Scenarios 5 and 17 is located further away from the intake; however, the associated excess temperature at the intake is still higher than the defined project criteria.This is due to the orientation of the ports as they are either partially or completely pointed toward the intake, which results in higher excess temperature.For all scenarios with maximum excess temperature exceeding the acceptable threshold (i.e., less than 1 • C), their plumes are presented in Figures 12 and 13.Almost 60 percent of the model scenarios show higher excess temperature during summer compared with the winter season.The remaining 40 percent of model results show identical excess temperature between the two seasons.Reduction in excess temperature from summer to winter varies from 8%-57%.The exact reason for the difference in results between summer and winter seasons is not entirely known, but it appears that the higher temperature of the discharged water results in a lower density of the plume as expected from the equation of state.This causes the plume to rise faster, and therefore, it reaches the intake inlet with less dispersion.Similar behavior and results have been observed by other researchers and engineers [11,31,32].

Summary and Conclusions
This study evaluates the transport of power plant heated water discharged to the source water body (i.e., Delaware Estuary).The goal is to finalize the outfall location through limiting the rise (i.e., excess) in water temperature at the power plant intake location and port orientations to equal or less than 1 • C. Numerous parameters are used in the momentum, continuity and transport equations to simulate that the physical processes pose an essential level of uncertainty in estimating the shape and magnitude of thermal plume in complex environments, such as estuaries.The combination of hydrodynamic calibration, thermal calibration, sensitivity tests, conservative assumptions and statistical analysis is used to establish the most conservative model parameters and ambient condition.A three-dimensional hydrodynamic and thermal computational model is developed using Delft3D-FLOW for a portion of the Delaware Estuary, encompassing the power plant intake and discharge systems.The overall model (i.e., larger model) that encompasses the nested model is calibrated for hydrodynamics using existing NOAA tidal and current gauges.The nested grid takes its hydrodynamic boundary condition from the overall grid.A temperature calibration is performed for the nested grid to satisfy proper selection of transport parameters.The intake withdrawal and heated outfall discharges are considered only in the nested model during the model scenarios for the winter and summer seasons with average water temperatures of 5 • C and 24 • C, respectively.The average wind condition with different speeds during the summer and winter seasons in the conservative direction of southwest-northeast is applied to all model scenarios.
The two primary conclusions of this study are: the power plant outfall should be located at least 1300 m from the shore; and outfall ports should be oriented between 90 • and 180 • from north, to reduce the thermal impact of the heated water to the intake to 1 • C.

Figure 1 .
Figure 1.Geographical extent of the study area.

Figure 2 .
Figure 2. Overall domain grid cell resolution, bathymetry and boundaries.
developed for a coastal flooding study.The best available data are used to construct the bathymetry by the Coastal Process Branch of the Flood and Storm Protection Division, U.S. Army Engineering Research and Development Center (ERDC).All elevations are converted from the North American Vertical Datum of 1988 (NAVD88) to the MSL.The QUICKIN module of Delft3D coupled with ArcGIS is used to interpolate bathymetry data into model grids.

Figure 3 .
Figure 3. Nested domain grid cell resolution, bathymetry, boundaries and location of intake/discharges.

Figure 4 .
Figure 4. Time series of simulated water levels at the intake inlet locations.Left: entire 546 days; right: 3-day period.

Figure 5 .
Figure 5.Total freshwater discharge at the inland boundary of the overall domain at Commodore Barry Bridge.

Figure 6 .
Figure 6.Wind time series and rose plot of the study area.

Figure 7 .
Figure 7. Scatter plot of modeled and measured water level at the Ship John Shoal NOAA Station.

Figure 8 .
Figure 8.Estimated and modeled temperature at a point in the middle of the nested grid.

Figure 9 .
Figure 9. Thermal plume pattern during tide at the surface: Scenario 3.

Figure 10 .
Figure 10.Thermal pattern during tide at the surface: Scenario 1.

Figure 11 .
Figure 11.Maximum excess temperature at the intake location for all model scenarios.

Table 1 .
Tidal constituent amplitude (m) at each side of open boundaries (see Figure2for the location).

Table 2 .
Freshwater sources and associated USGS gauges.

Table 3 .
General settings of model scenarios.

Table 4 .
Special settings of model scenarios.N stands for Normal; M-A for Momentum Alternating north and south ports; M-S, M-W and M-NE for Momentum with ports pointing towards the South, West and Northeast, respectively.