Characterising industrial thermal plumes in coastal regions using 3-D numerical simulations

Coastal power stations use sea water as a coolant, releasing it back into coastal environments at a higher-than-ambient temperature. Due to the possible ecological impacts on sensitive coastal zones, thermal plume formed by warmer coolant waters needs to be monitored, which is typically done through field campaigns. This paper assesses the use of simulations and remotely sensed observations as complimentary methods to characterise plume behaviour for a chosen coastal power station located within an inter-tidal embayment. Simulations of the thermal plume for two main tide phases and associated sea current conditions are validated against the high-resolution satellite observations. Simulated plume temperatures are higher than the observed values, with the biggest difference of 2 °C. The direction of the simulated plume dispersion is in agreement with observations and depends on the strength and direction of sea currents associated with the phase of the tide. The plume stretches most at the surface with limited impact on the benthic temperatures.


Introduction
Nuclear power stations are usually located next to large rivers or on the coasts and use river or sea water as a coolant for their installations. The effluent water is released back into the environment at a higher temperature than the ambient intake water. This coherent outfall water-mass at higher-than-ambient temperature will be referred to in this paper as a thermal plume. One of the impacts thermal plumes can have on the environment is the reduction of the soluble oxygen (20% for every 10°C increase in temperature) [1], which in turn changes the rate of photosynthesis, eutrophication and metabolic rate of fish [2]. Some aquatic organisms have narrow temperature tolerances and increase in water temperature causes deformation in fish species [3,4] or forces them to relocate to colder regions.
In order to fulfill environmental requirements, industrial power plants carry out studies prior to operation and continue to monitor thermal discharges during the life of the power plant. Monitoring of the industrial thermal plume has typically been done through field campaigns, which involve a survey vessel taking point measurements in the vicinity of the power plant outfall pipes. Such campaigns are expensive and can only be carried out during certain months with good weather conditions. Any issues with the measuring probes or sudden weather changes cause the campaign to be cancelled or postponed, which generates new costs.
With recent advances in technology, alternate ways of plume monitoring are emerging. High-resolution remote sensing observations provide good insight into spatial distributions of the sea surface temperatures (SSTs) in the vicinity of the industrial power stations [5][6][7][8]. Satellite observations provide a snapshot in time of the area of interest during the satellite overpass, during specific atmospheric and tidal conditions and outfall values from the power plant. Although good for understanding surface processes, satellite observations are not representative of the subsurface development of the thermal plume and hence its possible impact on the benthic environment. In order to gain broader understanding of possible impacts of industrial thermal plumes, high-resolution observations can be integrated with thermal plume simulations. Modelling of thermal discharges into aquatic environments allows us to gain understanding of the thermal plume behaviour sub-surface and its 3-D nature.
Calibrated numerical simulations have opened the opportunity to reproduce a specific instance of a thermal discharge [9,10], analyse multiple scenarios of thermal plume dispersion [11][12][13], or predict possible future thermal pollution prior to building a power station [14]. Shallow-water models, due to lower computational costs, are often used to investigate temperature distribution and thermal plume dispersion into river [15] or coastal environments [16]. Three dimensional (3-D) modelling is used for deep reservoirs [17] or to obtain a more detailed analysis of thermal discharge into coastal environment [12].
The area and the direction of the simulated thermal plume depend on multiple prescribed factors, such as: power station capacity [12], discharge velocity [15], sea currents [13], wind speed [9,12,13], tidal conditions [12]. In order to understand the plume behaviour and its possible effects, it is worthwhile to conduct multiple simulations with changing conditions and relate outfall temperature to the intake one. Due to long run times and simulation costs, modelling studies often need to focus on a few controlling variables that influence thermal plume dispersion.
The thermal plumes released into stagnant environments tend to rise to the surface and spread at the surface depending on the volume outfall rate, minimising the impact on the bottom temperatures [18]. For such sites it is worthwhile to run a set of simulations with varying discharge rate and power station capacity. Thermal plumes ejected into larger water bodies tend to be more reactive to winds [19,20] and tides [12,20] rather than operational capacity [19]. For those locations performing simulations with different current strengths and directions or varying wind speeds is of advantage.
Once the impact of the thermal plume is well understood, the simulations can become a tool for the mitigation strategy planning for existing power stations. In case the plume is found to critically impact the environment, future developments to decrease the impact of the plume, such as: barriers [21], relocation of the discharge pipe [22] or building an extra discharge pipe [15], can be simulated before undertaking any engineering works. In order to assess the effectiveness of proposed mitigation strategies and find the most costeffective solution, thermal plume simulations before and after the mitigation development can be compared.
In order to assume that the simulations are producing accurate results, of either current power station set up or future scenarios, it is vital to compare the output to observations or laboratory experiments. Most simulations are compared against field measurements [23][24][25], thermal sensor observations [18] or experimental data [19].
The purpose of this paper is to report additional means of characterising industrial thermal plumes in the coastal regions other than traditional field measurements. We use high resolution remote sensing imagery combined with a computational fluid dynamics (CDF) model. We detect thermal plumes in the high-resolution satellite imagery [5] and create 3-D simulations for selected days, simulating the tidal conditions during the satellite overpass. We use skill scores to compare the plume at the surface obtained through numerical simulations with the plume observed by the satellite sensor. We then assess any possible impact on surrounding sensitive habitats. This paper proceeds as follows: section 2 describes the methodology, section 3 presents the results, section 4 provides discussion and limitations, section 5 summarizes our findings.

Comparison of different plume characterisation methods
Traditionally, thermal plumes have been characterised through field campaigns. Boat surveys combining horizontal and vertical profiling gives an impression of how the plume is behaving with depth and the three dimensionality of the plume can be inferred. Field campaigns usually occur in the summer season. Moreover, field campaigns are not carried out regularly and are typically done once for each power station. The measurements are usually undertaken for one neap tide and one spring tide. The plume behaviour captured during those tide phases is assumed to be consistent throughout the whole year.
Taking measurements over the whole day, does not allow the boat survey to capture a complete instantaneous snapshot of the sea surface temperature (SST) distribution during e.g. mid-ebb tide. An uneven distribution of the profiles coupled with shallow bathymetry can compromise the extent of the plume sampling as measurements are only collected in the places deep enough for the survey vessel to move freely.
Freely available high-resolution satellite data provides complimentary information. Due to the large spatial domain, a single high-resolution (30-100 m) snapshot of clear-sky remotely sensed data provides information the surface temperature distribution over the whole region of interest during a single moment in time. Satellite records over several years can be analysed together to give information on plume behaviour over changing tidal phases, varying sea currents and seasons. Since satellite observations do not observe temperature distribution subsurface, using complimentary methods like 3-D numerical simulations is useful to understand the temperature distribution throughout the water column and at the seabed.

Satellite observations of thermal plumes
We use high-resolution infrared and visible observations from the Operational Land Imager (OLI) and the Thermal Infrared Sensor (TIRS) on-board Landsat 8 focused on our Region of Interest (ROI), the environs of Heysham. Heysham power station (54.0302°N, 2.9153°W) is located in Heysham, on the west coast of the UK, in the Morecambe Bay. This is one of the largest macro-tidal embayments in the UK with regularly exposed sandbanks and mudflats. We have chosen this ROI due to the bay's tidal regime and its associated dynamic landwater boundary.
Landsat 8 provides high-resolution observations from February 2013 until the time of writing, with a return time of 16 days. The number of observations is further limited by cloud cover, since only observations with clear sky or low cloud cover are useful in this application. We detected the extent of water using the Modified Normalised Water Difference Index (MNDWI) [26] due to the dynamic land-water boundary and presence of intertidal flats.
The TIRS instrument onboard Landsat 8 had large radiometric errors caused by the stray light entering the field of view (FOV). In 2017, a new calibration process with stray light correction has reduced the uncertainty from fluctuations in calibration to 0.51 K at 300 K for Band 10 (10.6-11.2 μm) and 0.84 K at 300 K for Band 11 (11.5-12.5 μm) [27]. Since the observations carry some uncertainty, there have been multiple validation practices developed and metrics used to quantify the uncertainty related to the satellite observations [28,29]. In order to investigate thermal contrasts in the imagery, we retrieved sea surface temperature (SST) over the water regions in every image using Optimal Estimation [30]. We applied a plume detection algorithm on the SST scenes using a combination of cloud and land mask expansion, surface temperature thresholds, warm water areal size and distance from power station outlet calculation. Full details of these steps can be found in Faulkner et al (2019) [5]. We obtained 46 clear-sky images over our region of interest for the time period 2013-2019. From the 46 Landsat 8 satellite scenes, we chose two cases that fitted different scenarios: mid-ebb tide and midflood tide.

Setting up thermal plume simulations
Our thermal plume simulations address clear-sky images obtained on 08/12/2017 and 17/03/2016. We performed the simulations for two different tide scenarios using FLOW-3D software [31,32]. Setting up simulations under the same conditions as the observations required the information on phase of the tide, which was obtained from National Tide Gauge Network (NTNG). For each observation and its tidal phase, we required an estimation of ambient sea currents associated with the tide. For this, we accessed sea current data from the Norwegian Meteorological Institute (MET Norway) [33].
In order to reproduce the local topography, we obtained the bathymetry from European Marine Observation and Data Network (EMODNet) and land elevation from ASTER Global Digital Elevation Model (ASTER GDEM). Land elevation was merged with the bathymetry and cropped to the area corresponding to the ROI using QGis software [34] and converted into a 3-D object using Rhinoceros software [35]. The 3-D object was then set up as the area's geometry for the simulations. Next, details such as outfall pipes were added to the geometry elements. The cooling water from Heysham Power Station is ejected through two parallel channels into the coastal waters. The information on the location and dimensions of the outfall channels is available through Google Maps open source software [36]. Because of the large tidal variations in the height of the surface water level and the water retreating during the ebb tide from a couple hundred meters away from the shoreline, the outfall channels are dug deep and long enough to make sure that the ejected plume reaches the ambient coastal waters under any conditions. In order to be computationally cost-effective the simulation domain was limited for each plume simulation, depending on where the observed plume was located in the Landsat 8 imagery (figure 1). For the ebb tide the thermal plume always disperses south-westward off the discharge pipes, which allowed us to position the active simulation domain in that region. For the flood tide, the simulation domain was positioned northwards of the discharge pipes, based on the plume location observed in the satellite imagery.
The computational mesh specifies the 3-D region for which the simulation will be performed, the domain outside the mesh is not included in the simulations. The mesh was comprised of 5 m × 5 m cells in the horizontal (x and y dimension) and the smallest cell size in the z direction that is allowed by the software solver, i.e. 2 m cell size for ebb simulations and 1.4 m for flood simulations. The size of the two domains in horizontal is similar, however, the ebb domain is much larger in the vertical compared to the flood domain. Hence reducing the cell size for the ebb domain would result in number of cells within the domain exceeding the permitted limit. The smaller cell size in the z direction for flood domain was essential to get the best possible depth resolution in the very shallow regions present north of the power station (identified by the blue box in figure 1). Once the geometry in the simulation domain was representing real life conditions, we added the ambient fluid, which as default for hydrodynamic simulations in FLOW-3D, is water at 20°C. In this study we were practically constrained to use the default ambient fluid. We specified the water surface elevation in the simulation based on the tidal observations. First, simulations were run to establish a flow of the ambient water through the domain. Values of the velocities associated with the tidal currents, obtained from the Norwegian Meteorological Institute (MET), were prescribed at the boundaries of the simulation mesh, where the flow was entering the domain. This means that the velocities were prescribed at the upper and the right boundary of the mesh for the ebb flow and at the lower and the left boundary of the mesh for the flood tide. The other two mesh boundaries in the simulations were set to enable the flow to enter and leave the domain, which in turn removes any possible effects of the shear velocity associated with the solid wall boundary. The flow simulations were run until the solution was close to reaching steady state with respect to those boundary conditions, and thereby approximating the dynamic flow for the corresponding phase of tide. The simulations were then continued with thermal plume injection prescribed at the location of the outfall pipes. The thermal plume injection represented the outfall volume rates and the temperatures of the outfall plume, information provided by the power station operators.
To address the constraint that the ambient water temperature was 20°C, in the simulations the outfall temperatures were scaled to obtain the same density contrast as that between measured outfall water and the observed ambient water temperature (figure 2).
First, we calculated corresponding densities for the observed ambient temperatures and the temperatures measured at the outfall pipes using an equation after McCutcheon et al (1993): where ρ is density and T is temperature of water. Then, we calculated the density difference for the ejected plume with respect to ambient water density, as shown by: Next, the density of simulated ambient water at 20°C was calculated using equation (1). The density difference δρ was then added onto the density of simulated ambient water density: Newly calculated density at outfall pipes (new ρ o ) in the simulation was then recalculated into corresponding temperature. By doing this we ensured that the density difference between ambient water and the outfall thermal effluent was consistent between observations and the simulation.
After the simulations have finished running, we used the density difference information to turn the simulations into a temperature field that can be compared with observations. For all simulation cells, where δρ = 0 we assigned the observed ambient temperature values. For all other simulation cells, temperature values were recalculated based on density difference from the observed ambient values. The caveat of this approach is that the relationship between temperature and density for water is not linear, hence the same density difference will correspond to different temperature differences. Figure 2 shows an example scenario, where the observed ambient temperature and measured temperature of the outfall water are both colder than 20°C. The difference in the temperature between the ambient waters, 5.6°C, and the measured outfall waters, 17.0°C , is 11.4°C . The corresponding density difference is 1.17 kg m −3 . To obtain the same density difference in the simulations, with the constraint of ambient waters being 20°C, the outfall water temperature is set to 25.05°C. Even though the density difference remains the same between the observed and the simulated temperatures, the gradient representing the change in density will be different (see figure 2).

Results
In sections below, we discuss our results for the fully developed simulated plumes during the ebb tide on 17/03/ 2016 ( Figure 3) and the flood tide on 08/12/2017 (figure 4). Land in the figures is shaded by the dark grey colour with ambient water regions in blue colour. The plume is marked by the black solid line. The contour lines within the plume represent the temperature gradation. In order to compare the simulations of the surface plume to the observations, the simulation results were regridded onto observation grid (50 m × 50 m) by taking an average over all simulation points that fall into one observation grid. The simulated plume was defined as ambient temperature + 1.5°C, to keep the criterion consistent between detected plumes in Landsat 8 observations and simulations. This criterion was first used for plume detection in the satellite imagery and was chosen in view of the uncertainty associated with the satellite instruments and the sea surface temperature retrieval method [5]. figure 3 presents horizontal cross sections of the simulated plume during ebb tide for: the water surface, −3.0 m, −5.0 m and −7.0 m depth contours. The plume is dispersing in the south-westward direction, which is consistent with the sea currents during ebb tide. At the time of the observation, the tide was in the mid-ebb phase, which means that the tidal currents were the strongest and consistently moved in the same direction. The tidal current velocity estimated during the observation time was 0.65 m s −1 .

Simulations of thermal plume dispersion for ebb tide
The simulated plume at the surface (figure 3 top left) disperses as a narrow band of 2.2 km in length and approximately 0.3 km in width. The plume does not disperse much laterally in the ambient waters, which is indicated by the steep temperature gradients along the length of the plume. The area, over which the temperature is elevated over 5°C higher than ambient, is spread south-westward up to 1 km away from the outfall pipes. The highest temperature noted in the simulation is 17.5°C, which is 7°C higher than ambient, and is located within 0.3 km of the outfall pipes. There is a small circular area in the south-east part of the plume, where a slight temperature rise of 3°C warmer than ambient is noted ( figure 3). Comparing the simulation results with the bathymetry map in figure 1, the same circular feature can be observed. This circulation feature is not shown in the -3.0 m cross section. The plume thus follows the bathymetry and does not disperse into shallow-water regions.    The plume observed at the surface (figure 4 top left) disperses as a narrow band of 2 km in length and approximately 0.4 km in width. The plume is mixing well with the ambient waters with distance away from the source due to the strong currents, which can be observed by the temperature gradient within the plume with highest temperatures located close to the outfall pipes and plume temperature decreasing away from the outfall pipes. The area of highest temperature rise (over 5.5°C higher than ambient) is limited to 0.4 km away from the outfall pipes. The highest temperature noted in the simulation results is 15.2°C, which is 9.6°C higher than ambient. The temperature difference between the ambient water and the plume in the simulation is higher than in the ebb case, which is related to higher outfall temperatures relative to the ambient waters recorded during 08/12/2017. The area with the highest temperature rise is limited to 0.1 km away from the outfall pipes.
The simulated plume is hugging the coastline and the extent of the plume is limited by the bathymetry in the shallow regions, which can be seen on the cross sections at −0.9 m (figure 4 upper right), −3.7 m (figure 4 lower left) and −6.5 m ( figure 4 lower right). The region of highest temperature rise is located close to the outfall pipes and its area decreases with depth, suggesting that the plume rises to the surface and spreads more at the surface. For −3.7 m and −6.5 m depth, the plume's area is limited to 1 km and 0.4 km in length respectively, with warmest temperatures around 5.4°C above ambient water temperature.
The plume is well mixed and the plume temperatures further from the outfall pipes are lower. The mixing is caused by the strong north-eastward tidal currents. There is a small recirculation area, where the south-east part of the plume enters the Heysham harbour. The coolant water intake point is located in the harbour, so this could cause potential challenges with intake of warmer plume waters for cooling purposes. However, the part of the plume entering the harbour is very narrow and the temperature rise in the harbour is only 1.5°C warmer than ambient, which is on the threshold of plume detection. At deeper levels, no plume recirculation into the Heysham harbour is present.

Comparing surface SST distribution from simulations with satellite observations
In order to assess the simulations, we compared the simulated surface thermal plume with the corresponding Landsat 8 observations for ebb ( Figure 5) and flood (figure 6) tides. The exposed sand banks are marked in light brown and the land in dark green with the power station as a red square. For both ebb and flood tide, there is good agreement between the observations and simulations in the direction of the plume dispersion (see table 1), suggesting that the plume is embedded in the tidal currents. The temperature distribution within the plume is similar between the observations and the simulations. However, the absolute temperature values are higher in the simulations than the observations, suggesting that the simulations overestimate the temperature (see table 3). The difference between the temperatures can be related to the density approximation.
For the ebb tide ( Figure 5), the difference in the plume area between the simulations and observations is greater than for the flood tide ( figure 6). The observed plume is more laterally dispersed in the observations, having approximate width of 0.9 km. The simulated plume is narrower with warmer temperatures highly concentrated over an area approximately 0.3 m wide. The area of highest temperature rise ( 15°C) is similar in both observations and simulations, however, the maximum temperature values are higher in the simulations by 2.3°C compared to the observed maximum of 15.2°C. The simulated plume is constrained by the bathymetry and the location of the coastline, however, it turns westward further away from the outfall pipes, which is not observed in the satellite image. The plume turning away from the bathymetry further south is not caused by boundary effects, since the currents were prescribed only at two boundaries. However, the flow development in the simulation is limited to the area filled with water, which is represented by a channel widening to the west (see figure 3 top left), and the velocity vectors redirect towards the west from their prescribed south-westward flow.
Looking at the thermal plume direction for the flood tide (figure 6), the simulations and observations are in good agreement. The plume is dispersing as a narrow band north-eastward, reaching over 2 km in length. The simulated plume is narrower than the observed plume and the difference is greatest further away from the outfall pipes. A possible explanation is that at the time of the satellite overpass the wind was blowing north-westward at 8.5 m s −1 , and the influence of wind stress was not included in the simulations. The temperature distribution within the simulated surface plume resembles the one of the observed plume, with highest temperatures located close to the outfall pipes (0.3 km away from the outfall pipes in the observations and 0.2 km in the simulations). The highest temperature in the observed plume is 13.0°C, while in the simulations the maximum temperature within the plume is 15.2°C. The highest temperature region is limited in area and further away from the source the surface plume temperatures are cooler, which is caused by stronger tidal currents.
In order to compare the surface plume between the simulations and the observations we have used skill scores. Skill scores are a common measure used to measure the degree of accordance between the observations and the prediction, usually achieved through a model or simulation. Skill scores used in this paper are: Hit Rate, which indicates the percentage of simulated plume agreeing with the observed plume, False Alarm Rate, which indicates the percentage of simulated plume falling outside of the observed plume, and Misses, which show the percentage of observed plume, which is absent in the simulations. All skill scores are summarised in table 2.
In case of the ebb tide ( figure 7 a), the simulated plume area is located almost completely within the observed plume area, with 13% of the simulated plume not coincident with the observed plume. However, the simulated plume area is almost two times narrower, resulting in 55% of the observed plume lying outside the simulated. For the flood tide (figure 7 b), the simulated plume is dispersing slightly closer to the coast than the observed plume, with 19% of the simulated plume not consistent with the observed plume. As in the previous case, the simulated plume has a smaller area, resulting in 54% of the simulated and observed plume being in agreement and 46% of the observed plume being missed in the simulations.
The simulations capture the direction of the plume dispersion, but underestimate the areal extent at the surface ( figure 7). Because the plume is injected into an already steady state flow, it will continue dispersing in the same direction. In reality, the plume is dispersing constantly into a changing tidal environment and during transition phases it is forced away from the coast and disperses horizontally. It is possible that due to the history Figure 6. Observed (left) and simulated (right) surface plume during mid-flood tide for Heysham. The exposed sand banks are represented by the light brown colour and land by the dark green with the power station marked as a red square. Bathymetry is marked in thin black lines. The plume's temperature contours are thick lines marked also on the colour bar. of the plume over time, the plume will be more spread out at the surface in observations than in the simulations, which assume only equilibrium state.

Discussion
The simulated results are in part shaped by the assumptions and approximations made, and these considerations inform the interpretation of the results in combination with the satellite observations. The principle assumptions and constraints are: • Bathymetry, which impacts the development of the flow, is represented sufficiently by the highest possible simulation resolution.
• Simulated dynamics are adequately represented by matching the density contrast between discharge and ambient waters to that of the observed temperatures.
• Heat transfer in the simulations is limited to the exchange between the ambient water body and the discharged plume, air-sea heat exchange is neglected (radiative heat fluxes, latent heat fluxes and sensible heat fluxes).
• The simulations capture simplified tidal conditions with tidal currents prescribed in one direction for each simulation.
These are now discussed in turn.

Simulation analysis
Our choice of simulation mesh size permits most bathymetry features to be resolved, especially in the deeper regions. Discrepancies between the actual seabed and the simulated seabed will only occur in the very shallow regions, where the depth difference across an area is smaller than 1.4 m. The simulated geometry of the domain will influence the flow development, for which velocity vector directions and strengths are prescribed at  boundary as prior. This can be seen during the simulated ebb tide, when the water is diverging slightly westward in further parts of the plume, however, the main axis across the most part of the plume in both cases is following fairly close. This shows that for the complicated bathymetry, like Heysham, the simulations are able to capture most seabed features and simulate the flow plausibly. In order to capture the shallow coastal inter-tidal environment of Heysham, we have set up the highest possible horizontal and vertical resolution, which does not simplify the geometry to flat planes and allows for replicating the bathymetry of the region. It is important to remember that idealised simulations do not fully reproduce the actual conditions on the observed day and the reasons for subtle differences between the simulated and observed state should be analysed. The plume and ambient temperature are based on density contrasts rather than actual observed values. Conversion between the temperature and density was possible due to the governing equations employed in the FLOW-3D, where the density is evaluated as a function of temperature and fluid fraction. However, as shown in the Methods section, the relationship between density and temperature is not linear. This means that the observed density difference will correspond to a different range of simulated temperatures during the simulation time. This assumption does not fully capture the seasonality, but it enables account to be taken of the impact of seasonality. The ambient temperature is always the same value and the discharged plume temperatures are recalculated by preserving the density contrasts with respect to the fixed ambient temperature off the simulation software. For seasons with warmer observed ambient temperatures, the gradient on the temperature-density curve will be more similar to the gradient in the simulations, hence capturing the density contrasts and temperature contrasts closer to the observed values. For winter months, where the observed ambient water temperature is lower than the fixed simulated ambient temperature, the temperature-density gradient will be flatter, corresponding to a larger temperature difference. Table 3 below summarises the conditions for both days, for which the simulations were carried out. For the ebb tide simulation (17/03/2016) the ambient water temperature was 10.5°C, the average temperature rise within the plume was 15.5°C, the air temperature was 7.0°C and the wind speed measured closest to the satellite overpass was 1.0 m s −1 . For the flood tide observations the ambient water temperature was 5.6°C, the area of significant temperature rise within the plume was 11.1°C, the air temperature was 3.0°C and the wind speed measured closest to the satellite overpass was 8.5 m s −1 .
In order to investigate the differences in the plume's temperature between the simulations and the observations, we looked into heat transport within the water body and the heat exchange between the water and the air. A limitation of the simulations is that heat transfer is occurring between the ambient waters and the discharged plume through turbulent diffusion and advection of the plume. The simulations do not include airsea heat exchange, since they are performed for a single fluid (i.e. water) and hence surface heat loss is not captured. Surface heat loss between the water and the atmosphere occurs through radiative heat fluxes, latent heat fluxes and sensible heat fluxes, none of which are included in the simulations. We looked into the rate of heat loss associated with each of those fluxes and calculated the excess rate of heat loss associated only with the thermal plume. We calculated the excess fluxes to assess the effect the injection of the thermal plume has on the surface heat loss compared to a scenario including only ambient water and no thermal effluent.
The radiative heat flux is related to the water temperature and can be parametrised as: where Q R is the emissivity of the water surface (we used the value of 0.97), T s is the sea surface temperature, ε a is the broadband effective emissivity of the atmosphere assumed to be 0.75 after [37] , T a is the pressure-weighted average atmosphere temperature and σ is the Stefan-Boltzman constant.
In order to calculate the excess radiative heat flux, we modified the above equation to get the difference between the radiative flux over the plume and the radiative flux over ambient water: where ¢ Q R is the excess radiative heat flux caused by the presence of the thermal plume, T p is the average temperature across the plume and T w is the ambient water temperature. The excess radiative heat flux caused by the excess heat from the thermal plume for the ebb tide on 17/03/ 2016 was 25.51 W m −2 . The excess radiative heat flux value calculated for the flood tide on 08/12/2017 was 26.71 W m −2 .
We also calculated sensible and latent heat fluxes, which describe the energy being released from the Earth's surface and absorbed in the atmosphere. Sensible heat flux is related to changes in temperature by heating an object directly, while latent heat flux is related to changes in phase associated in this case with evaporation. Sensible heat flux, Q S , is dependent on wind speed and temperature [38] as: where c s is the stability dependent bulk transfer coefficient, U 10 is the wind at 10 m height, T s is the SST and T a is the temperature of the air at 10 m height. For c p we used 1004.83 J kg −1 K −1 . For c s we used a typical value of 0.001 75 after [39].
In order to calculate the excess sensible heat flux, we modified the above equation to: which yields: where ¢ Q S is the excess sensible heat flux caused by the presence of the thermal plume, T p is the average temperature across the plume and T w is the ambient water temperature.
The calculated excess sensible heat flux for the ebb tide was 10.86 W m −2 . For the flood tide the excess sensible heat flux was higher (112.72 W m −2 ), which is related to much stronger winds recorded during that time (see table 3).
Latent heat flux, Q E , is related to wind fluxes and differences in specific humidity, using bulk aerodynamic parametrisation, as: where C E is the bulk transfer coefficient for water vapour, U 10 is the wind speed at 10 m, q a is the near surface airspecific humidity, q s is the air-specific humidity just at the air-ocean interference.
In order to calculate the excess latent heat flux, we modified the above equation to: 10 10 which yields: where ¢ Q S is the excess sensible flux caused by the presence of the thermal plume, q p is the air-specific humidity across the plume and q w is the air-specific humidity over ambient water.
For the ebb tide the excess latent heat flux was 33.67 W m −2 . For the flood tide the excess latent heat flux was 118.18 W m −2 . Higher values of the excess latent heat flux for 08/12/2017 compared to 17/03/2016are, like in the case of sensible heat flux, related to much stronger winds observed on that day (see table 3).
Next, we looked into the excess rate of cooling of the water column related to the presence of the thermal plumes. We related the net rate of heating out of the water column to the rate of change of heat content by: where F is the sum of the fluxes out of the water column, A is the area, ρ is the density of water, h is the height of the water column, c p is the specific heat capacity of water, ¶ ¶ T t is the temperature change with time. After rearranging, the above equation yields:  (see table 4), hence not including the heat fluxes in the simulation will not impact the results substantially. Our simulation set up works well for domains, where the driving factor for heat dispersion are the sea currents. The effect of air-sea heat exchange for a small and shallow domain plays a much smaller role compared to diffusivity of the plume within the ambient water and therefore can be neglected. The simulations capture the plume evolution and dispersion in the ambient waters well and provide information on the vertical structure of the temperature and plumes impact on the seabed temperatures. Such analysis is useful for environmental monitoring the plume.
In the future, setting up a simulation with varying tidal and sea currents conditions would be an improvement. In reality the plume is ejected constantly into a changing environment and capturing the history of the plume over a changing tidal cycle would be more complex. Simulating plume discharge into a changing tidal cycle would result in a more spatially dispersed plume compared to the simulations of one steady flow tidal phase, due to tidal currents changing direction. However, for the operational characterising purposes, the simplified simulation is a good approximation.

Exploring the impact of seasonality in the region of interest (ROI)
Given the seasonal cycle of weather, there is potential for seasonal effects in plume dispersal. The seasonality of plume dispersal can be explored using the satellite observations. Figure 8 presents the differences in the water temperature between the ambient coastal waters and the released thermal plume for observations collected using Landsat 8 [5] over four seasons. The ambient water temperature for each observation was calculated as the average across an area of 5km × 5 km of shallow coastal waters not affected by the discharged plume. The maximum temperature within the plume was chosen by identifying the correct area of raised temperatures within the outfall pipes and finding the highest value. Blue points represent the temperature differences observed during winter months (December, January, February), green-spring months (March, April, May), yellow-summer months (June, July, August), red-autumn months (September, October, November). The ambient water temperatures are highest (above 11°C) in the summer time, late spring and early autumn, while during the winter the ambient coastal waters have temperature oscillating between 3°C and 6°C. In 85% of the observations the maximum temperature difference between the plume and the ambient waters was 4-8°C irrespective of the season. The highest temperature differences occurred when the plume spread was limited.  Since the study focuses on simulations for two tidal phases: mid-ebb and mid-flood, we processed observational data to identify how plausible the presented plume dispersion is considering other seasons. Out of the 46 Landsat 8 scenes, we have selected 22 images that captured the plume dispersion during the mid-ebb and mid-flood phases captured from 2013 until 2019 during different months and seasons [5]. There was a clear bimodal behaviour of the plume dispersion. The plume was observed to always disperse south-westward during the ebb tide, following the deep water channel in that part of the bathymetry. During the mid-flood tide the plume dispersed north-eastward into the shallow water regions. The extent of the observed plume was predominantly related to the strength of the currents and no seasonal modulation of that dominance is apparent in the data.

Plume's impact on neighbouring infaunal species
The simulations prior to the steady state capture the evolution of the plume from the injection of the plume into the ambient waters. Typically, the plume is injected continuously into the ambient waters, but this is useful for occasions, when the power station is turned off and then returns to operation. Such simulations also prove useful in assessing the possible impact of the thermal plume prior to the power station operation, which is often essential for obtaining environmental permits.
Morecambe Bay is the meeting point of major estuaries: Kent, Leven, Lune, Wyre, and multiple smaller estuaries. Together they form the largest continuous intertidal area in the UK [40]. The tidal currents are the main driving force within the bay. The 3-D thermal plume simulations for the chosen mid-ebb and mid-flood scenarios capture the direction, extent and subsurface plume dispersion and are plausibly representative of other ebb and flood tides. The main variable factor is the strength of the tidal currents, as during mid-tide phases the currents are the strongest. This suggests that during the beginning and towards the end of each tide phase, when sea currents are weaker, the plume does not disperse as far away from the discharge pipes.
The whole of Morecambe Bay is under protection as the Special Area of Conservation (SAC) for habitats. The intertidal sandflats located in the shallow parts of the bay are part of the Special Protection Area (SPA) for bird species. This means that the bay area is valuable in terms of British wildlife and diverse habitats. Due to the high ecological importance and its location with respect to SAC and SPA, any future development proposals, outside the already existing energy generation, are likely to have constraints on type and scale of the project [41]. Subtidal boulder and cobble skear communities as part of the SAC are located within 5km north of the power station. Figure 9 presents habitats of common mussel, common cockle, as well as boulder and cobble skear. The circular rings in the image present area with a radius ranging from 1 km to 5 km away from the outfall pipes. Mussel and cockle beds are widely distributed around the UK coast, preferring sandy bays and estuaries as their habitat.
Boulder and cobble skear in the Bay are habitats for such communities as serrated wrack, sponges, sea squirts and red seaweeds, as well as a good settlement for honeycomb worm (Sabellaria alveolata) reefs [42]. Honeycomb worm is is an eurythermic organism and can survive in temperatures ranging from 5°C up to 25°C. Exposing species to a stepwise temperature increase from 15°C to 25°C indicated that the worm is able to undergo biological changes in phospholipids to acclimatise to higher temperature locally [43]. For both simulated scenarios the temperature within the plume did not exceed 17.5°C and therefore it is unlikely that it had any impact on the honeycomb habitats.
The intertidal mudflats and sandflats support a range of infaunal species, such as Mytilus Edulis mussels and Cerastoderma Edule cockles, which prefer shallow depths up to 5 m. Common mussels (Mytilus edulis) have a large temperature tolerance, with temperature ranging between 10°C and 20°C not impacting the species [44], with other studies suggesting even larger tolerance of 5°C-0°C with upper sustained limit of 29°C [45]. Common cockles (Cerastoderma edule) are widely distributed around the British coast, and similarly to mussels they are eurythermic, which means they have a wide range of temperatures they tolerate with an upper lethal limit of 35°C after being exposed for 24 hours [46].
Mussel beds are located from 3 km to 10 km away from the power station. Cockle beds are located closer, present vastly in the shallow water regions, however, with none present in the deep water channel. Both mussel and cockle beds are located further away from the power station than the extent of the largest part of the thermal plume. In the south, the species are located in the shallow water regions, whereas the plume follows the southwestward deep water channel, so there is no intersection of the thermal plume and the habitats. The only possible risk for the infaunal species is in a scenario, where the plume dispersing northwards would extend further away from the outfall pipes. Since the plume dispersing northwards is well mixed and temperatures within the plume decrease rapidly away from the outfall pipes, the plume temperature would most likely not exceed 1.5°C warmer than ambient in those areas. Such temperature rise would not result in exceeding the thermal limits and therefore no lethal harm to the mussel and cockle beds.
Apart from the stationary habitats located in the bay, due to the presence of multiple estuaries, the bay is abundant in fish species, such as: plaice, flounder and dab. Moreover, it is also a vital fish nursery area and habitat for migratory fish species: salmon (Salmo salar), sea trout (Salmo trutta) and eels (Anguilla anguilla) [47]. Optimal temperature ranges between 0°C-28°C and 0°C-26°C for salmon and trout respectively, with optimal growth temperature of 16°C-17°C [48]. For eels [49] proposed that the European species can adapt to temperatures ranging from 0°C to 30°C. Since all river estuaries are located further than 5 km away from the thermal plume discharge point, the presence of the plume should not impact the fish species.

Conclusions
We have simulated the dispersion of the thermal plume at Heysham for two predominant, reoccurring scenarios: flood and ebb tide. In order to validate the simulation outputs, we have used Landsat 8 thermal imagery and a sea currents reanalysis dataset from Norwegian Meteorological Institute. The simulations were run separately for each tidal scenario. The main driver of the heat dispersion were the sea currents, with the airsea heat exchange playing a minor role. Simulations of thermal plume ejection into the coastal environment, validated by high-resolution observations, can be used as a good approximation for plume characterisation. Figure 9. Location of sensitive habitats in the vicinity of the Heysham power station. Boulder and cobble skear is marked in brown, mussel beds in grey, cockle beds in yellow. Left panel shows location of all communities, middle panel-boulder and cobble skear and right panel mussel and cockle beds. Distance away from the outfall pipes is marked by different sized circles with the radius written for each circle, ranging from 1 km to 5 km.
The simulations provide knowledge on plume development from the start of the operation until the simulation reaches steady state for each tide type. This includes specifying the main direction of the plume dispersion under different tidal conditions and possible extent of the plume. Performing simulations in 3-D provides information on subsurface dispersion and the 3-D nature of the plume. Such information can be used in understanding possible ecological impacts on the close-to-coast and benthic regions, and ensuring environmental requirements are met.
For small domains, like Heysham, idealised simulations driven by the tidal currents present a sensible solution compared to single-day field surveys. An advantage of simulations of thermal plumes over field campaigns is that the simulations enable investigating multiple scenarios for plume dispersion. Simplified simulations are a cost-effective way of understanding plume behaviour during different outfall rates and varying tidal conditions. A one day boat survey is only able to capture the plume only during specific atmospheric and hydrologic conditions and thermal plume outfall rates.
Lastly, simulations open a possibility of exploring scenarios with highest permitted outfall temperatures and volume outfall rates without the need for those worst case scenarios to occur in real life. Setting up simulations for maximum flow rates during different tidal conditions is useful in assessing areas at risk of being influenced by the thermal plume and explore potential ecological impacts on neighbouring sensitive habitats. Understanding possible impacts of the thermal stress on the temperature sensitive aquatic ecosystems in the vicinity of the power station is important from an environmental and sustainability point of view, to ensure the operation of the power station does not have a negative impact on the biodiversity.