Monitoring of Fine-Scale Warm Drain-Off Water from Nuclear Power Stations in the Daya Bay Based on Landsat 8 Data

Abstract

Monitoring the drain-off water from nuclear power stations by high-resolution remote sensing satellites is of great significance for ensuring the safe operation of nuclear power stations and monitoring environmental changes. In order to select the optimal algorithm for Landsat 8 Thermal Infrared Sensor (TIRS) data to monitor warm drain-off water from the Daya Bay Nuclear Power Station (DNPS) and the Ling Ao Nuclear Power Station (LNPS) located on the southern coast of China, this study applies the edge detection method to remove stripes and produces estimates of four Sea Surface Temperature (SST) inversion methods, the Radiation Transfer Equation Method (RTM), the Single Channel algorithm (SC), the Mono Window algorithm (MW) and the Split Window algorithm (SW), using the buoy and Minimum Orbit Intersection Distances (MOIDS) SST data. Among the four algorithms, the SST from the SW algorithm is the most consistent with the buoy, the MODIS SST, the ERA-Interim and the Optimum Interpolation Sea Surface Temperature (OISST). Based on the SST retrieved from the SW algorithm, the tidal currents calculated by the Finite-Volume Coastal Ocean Model (FVCOM) and winds from ERA-Interim, the distribution of the warm drain-off from the two nuclear power stations is analyzed. First, warm drain-off water is mainly distributed in a fan-shaped area from the two nuclear power stations to the center of the Daya Bay. The SST of the warm drain-off is about 1–4 °C higher than the surrounding water and exceeds 6 °C at the drain-off outfall. Second, the tide determines the shape and distribution characteristics of the warm drain-off area. The warm drain-off water flows to the northeast during the flood tide. During the ebb tide, the warm drain-off water flows toward the southwest direction as the tide flows toward the bay mouth, forming a fan-shaped area. Moreover, the temperature increase intensity in the combined discharge channel during the flood tide is lower than that during the ebb tide, and the low temperature rising area during the flood tide is smaller than that during the ebb tide.

1. Introduction

With the increase in energy demands and the development of nuclear energy technology, China has established a number of nuclear power stations in coastal areas, which have greatly promoted the economic development of those regions. However, the efficiency of nuclear power stations is lower than that of thermal power stations. Only 30% to 35% of nuclear energy is converted into electrical energy [1]. Most of the remaining energy is discharged with cooling water in the form of thermal energy, which causes the temperature of the surrounding water to rise rapidly, thus directly or indirectly affecting the growth and reproduction of aquatic organisms [2]. Algae, in particular, are sensitive to temperature variations. The rise of sea temperature may lead to an abnormal outbreak of algae. Driven by wind and tide, algae may invade the water-cooling system, block the discharge of the nuclear power stations, and eventually obstruct the safe operation of such stations [3]. Therefore, mapping and monitoring this warm drain-off water is of great significance to ensure the health of the coastal environment and ecosystem.

The warm drain-off water of coastal power stations is usually monitored and estimated by numerical simulation, in situ observations and thermal infrared remote sensing. The accuracy of numerical simulations mainly depends on numerical models and hydrothermal parameters. Field observations cost a great deal of manpower and resources. In addition, the limited sampling range and discrete data of such approaches cannot clearly display the spatial variation distribution of warm drain-off water. However, due to its advantages of rapid, economical and multi-scale analysis of heat distribution change regions, thermal infrared remote sensing technology can reflect the spatial and temporal differences and changes of water temperature. Therefore, this method is an effective means to monitor and estimate the warm plumes of nuclear power stations. Tang et al. [4] used the Advanced Very High-Resolution Radiometer (AVHRR) data for monitoring the warm plume of the Daya Bay Nuclear Power Station (DNPS) in different seasons, and the results showed that warm plumes could make the SST increase by 1~3 °C. Ahn et al. [5] indicated that the poor spatial resolution (1 km) of the sensors is insufficient to measure small-scale coastal warm plumes. Meanwhile, the authors used the Landsat 5 Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper Plus (ETM+) to seasonally observe the warm plume of the Younggwang nuclear power station. The winter thermal plume structure is significantly smaller than that during the summer, because of the input of cold water from the river in the winter. Zhou et al. [6] demonstrated the feasibility of HJ-1B Infrared Scanner (IRS) and FY-3A Medium Resolution Spectral Imager (MERSI) for estimating the sea surface temperature (SST) in coastal waters. They found that the area where SST had increased more than 3 °C in winter was much larger than that in the summer. The warm drain-off water spread along the coast in the summer and extended into a fan-shape area in the winter. By applying the single-channel algorithm to retrieve the SST from the Landsat 8 Thermal Infrared Sensor (TIRS), Cheng et al. [7] used tide currents and wind data to analyze the thermal diffusion characteristics of the Yangjiang nuclear power station in different seasons and pointed out that the tide was the main reason for the seasonal differences in the warm drain-off water. TIR data with a high spatial resolution have been used to investigate the warm drain-off water from coastal power stations by analyzing the spatiotemporal distribution of the SSTs retrieved from HJ-1B/IRS, FY-3A/MERSI, Landsat5/TM and Landsat7/ETM+. These sensors with a single-channel design are subject to an immature single-channel algorithm, which leads to uncertainty in their results. Landsat 8/TIRS consists of the thermal infrared Band 10 (10.60~11.19 µm) and Band 11 (11.50~12.51 µm). With more channels and a higher spatial resolution, Landsat 8/TIRS could be used as an ideal data source for monitoring the thermal characteristics of drain-off water.

The TIRS suffers from a serious stray light problem [8]. Non-uniform banding artifacts and the absolute radiometric calibration errors in Band 11 are two times higher than those of Band 10 [9]. Thus, previous studies have often used a single channel inversion algorithm for SST. In April 2017, the United States Geological Survey (USGS) achieved stray light correction for TIRS images, which greatly improved TIRS accuracy. The error in Band 10 was reduced from 2.1 K@300K to 0.3 K, while the error in Band 11 was reduced from 4.4 K to 0.19 K [10]. The split-window method has been applied to TIRS data and the feasibility of split-window atmospheric correction has been demonstrated [11,12,13]. However, these studies mostly focused on the land rather than the ocean.

Therefore, in order to explore the feasibility of split-window SST inversion and to select the optimal model to monitor warm drain-off water, this study estimates the Radiation Transfer Equation Method (RTM), the Single Channel algorithm (SC), the Mono Window algorithm (MW) and the Split Window algorithm (SW) based on Landsat 8/TIRS data and discusses the distribution characteristics of warm drain-off from nuclear power stations.

2. Materials and Methods

2.1. Study Area

The first large commercial nuclear power station in China, the Daya Bay Nuclear Power Station (DNPS), has been located at 22° 36′ 2.7″N, 114° 32′ 57.75″E (located in the southern part of Guangdong Province) since 1994 (Figure 1). The Ling Ao Nuclear Power Station (LNPS) is the eastern side of the DNPS. There are six million kW pressurized-water reactor nuclear power units at these stations. Two discharge channels are combined together towards the east, and the total cooling water discharge is more than 319 m³/s.

The DNPS and the LNPS are located on the western coast of the Daya Bay, which is a semi-enclosed area of about 600 km² in a subtropical monsoon zone with monsoon climate characteristics, four distinct seasons, and the characteristics of a maritime climate. The annual average air temperature is around 22 °C [14]. Affected by monsoons and terrain, the wind field in the Daya Bay has seasonal characteristics. Northeastern wind prevails in winter and southwestern wind prevails in summer. The wind speed in spring and early summer is relatively high, and the monthly-average wind speed is about 5.0~5.4 m/s. The wind speed in summer and winter is relatively small and the monthly-average wind speed is 4.6~4.8 m/s.

The tide is irregular semi-diurnal with an average height of 1.01 m and a maximum height of 2.5 m [15]. The tidal currents are the reciprocating movement, with velocities during the ebb tide greater than those during the rising tide. The current velocity in the Daya Bay decreases gradually from the mouth of the bay to the north, and the velocity at the top of the bay is smaller. The current in the western Daya Bay is dominated by tidal currents. Except for the strong velocity in the Eastern area, most of the sea area presents weak current characteristics. The flow velocity is about 30 cm/s in the Eastern area and 20 cm/s in the western area, while it is 5–10 cm/s near the station site [16]. The current is convergent and divergent near the nuclear power station, as shown in Figure 2. During the flood tide, the tidal current diverges near the site; one flows west to the Dapeng AO, and the other flows northeast along the coast. At the ebb tide, the western and northeastern tidal currents converge near the nuclear power station and flow south to the entrance of the Daya Bay. The water depth near the outlet of the nuclear power station is only 5 m.

2.2. Data

The data used in this study include Landsat 8 images, MODIS products, buoy data, the Optimum Interpolation Sea Surface Temperature (OISST), ERA-Interim data and tidal current data.

Landsat 8 was successfully launched on February 11, 2013. It offers a clearer view and better spatial resolution than most ocean sensing instruments, and is more sensitive to brightness and color than previous Landsat satellites. It carries two instruments: The Operational Land Imager (OLI) and the TIRS. Landsat 8/OLI collects data in visible, near infrared, and shortwave infrared wavelength bands, as well as in a panchromatic band. Panchromatic and multispectral data are taken at in 15 m and 30 m resolutions, respectively. Landsat 8/TIRS consists of the thermal infrared Band 10 (10.60~11.19 µm) and Band 11 (11.50~12.51 µm). The 100-m spatial resolution of the TIRS data is resampled to the OLI data to create radiometrically and geometrically calibrated terrain-corrected 30-m spatial resolution Level 1 data products [17]. The EarthExplorer (https://earthexplorer.usgs.gov/) offers Landsat 8 Level 1 Products and metadata. A total of 56 scenes of Landsat 8 satellite images in the South China Sea from January 2017 to January 2019 are used in the calculation of the split window model parameters. Figure 3 displays the locations of these satellite images. Eight Landsat 8 satellite images are used to analyze the distribution of the warm drain-off water under different tides and wind, as listed in

Table 1. Landsat 8 image list used in analyzing the warm drain-off water distribution of the DNPS and LNPS.

No.	Date	Time (UTC)	Landsat Scene ID
1	2017.02.18	2:45:51	LC81210442017049LGN00
2	2017.04.07	2:45:42	LC81210442017097LGN00
3	2017.08.29	2:46:02	LC81210442017241LGN00
4	2017.11.01	2:45:14	LC81210442017305LGN00
5	2018.03.09	2:45:34	LC81210442018068LGN00
6	2018.07.31	2:45:13	LC81210442018212LGN00
7	2018.10.03	2:45:41	LC81210442018276LGN00
8	2019.01.23	2:45:44	LC81210442019023LGN00

.

MODIS SST products with spatial resolutions of 1km are obtained from the Ocean Color (https://oceancolor.gsfc.nasa.gov/). The time difference between MODIS and Landsat 8 imaging is about 2 h. The operational algorithm of daily MODIS SST is a nonlinear SST algorithm (NLSST), which has 5 quality levels (0~4), corresponding to best, good, questionable, bad, and not processed. Only high-quality MODIS SST data (quality level marked 0 or 1) are selected in this study. The MODIS SST data are of good quality in the offshore area of China, with a mean error of 0.02 °C [18].

The SSTs from the buoy are provided by the Guangdong Ocean Monitoring Center for December 2017 to January 2019, with a temporal resolution of 30 min, which has lower temporal difference than MODIS SST. The locations of the buoy data are indicated in Figure 3 and

Table 2. Information of the buoy data.

No.	Station	Lat	Lon	Available Time
1	Leizhouxi	20.49	109.46	2018.03~2019.01
2	bohe	21.16	110.98	2017.12~2019.01
3	Zhapo	21.52	112.22	2017.12~2019.01
4	Honghaiwan	22.29	115.24	2017.12~2019.01
5	Zhelang	22.70	115.60	2017.12~2019.01

.

The OISST is an analysis constructed by combining observations from different platforms (satellites, ships, buoys) on a regular global grid, provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA. This study used the NOAA OI SST V2 High Resolution Dataset with a 0.25° spatial resolution using satellite SSTs from AVHRR (https://www.esrl.noaa.gov/psd/data/gridded/tables/sst.html).

ERA-Interim is the reanalysis data from ERA-40 provided by the European Center for Medium Range Weather Forecasts (ECMWF) [19]. It contains global reanalysis data since 1979 and was kept updated until August 2019. The data assimilation system used to generate ERA-Interim is based on cycle 31r2 of ECMWF’s Integrated Forecast System (IFS) published in 2006. This system includes a 4-dimensional variational analysis (4D-Var) with a 12-h analysis window. The spatial resolution of the data set is approximately 80 km (T255 spectral) at 60 vertical levels from the surface, up to 0.1 hPa. This study used SST, the near-surface temperature, and wind field data with a 0.125° spatial resolution at the surface level.

The Daya Bay tidal current data are based on the Finite-Volume Coastal Ocean Model (FVCOM). The FVCOM is a near-shore ocean numerical model developed by Chen et al. [20]. This model adopts an unstructured triangular mesh in the horizontal direction, σ coordinate transformation in the vertical direction, and the finite volume method for numerical purposes. The governing equation can be solved numerically by integrating the triangle control volume flux. Moreover, the model has the advantages of fitting a zigzag shoreline with the finite element method and efficient calculations with the finite difference method, which guarantees the conservation of regional mass and the momentum of the complex terrain and the shoreline, such as the areas near the shore and the estuary. The FVCOM model also includes 3D dry/wet grid processing modules, which can be used to deal with overland problems. Therefore, the FVCOM model is suitable for numerical simulation studies of the Daya Bay. The model simulations are consistent with the RMSE of the in situ data (about 5cm), and the relative error of current velocity is less than 15% [21]. Figure 2 shows the tidal current distribution of the Daya Bay simulated by FVCOM, which shows that the current is convergent and divergent with tides near the nuclear power stations.

2.3. TIRS Data Pre-Processing

The Landsat 8 Level 1 data available to users feature radiometrically and geometrically corrected images in the units of digital numbers. The SST inversion algorithm requires the clean spectral radiance and brightness temperature as input parameters, so we need to perform pre-processing steps, such as radiometric calibration, geometric correction, land-sea separation, and cloud detection. The methods and parameters involved in pre-processing can be obtained and calculated according to the Landsat 8 data user’s handbook. The removal of artifacts in Band 11 is also an essential step in studying small-scale environments using the SW algorithm.

In April 2017, the USGS achieved stray light correction for TIRS images. The amount of stray light in the scene is estimated using the scene obtained before and after the target scene and the edge pixels of the target scene. The stray light estimation is then subtracted from the target scene to obtain the corrected image. This method corrects the stray light and reduces the absolute radiation error. However, there are still the artifacts in the TIRS, especially in Band 11. The variability of the error in Band 11 is about two times higher than that in Band 10. The stripe in Band 11 is especially noticeable in the areas of a scene that are expected to be uniform (Figure 4). This leads to the stripe problem in the inversion results of the SW algorithm when using Band 11.

The artifacts in Band 11 are clearly visible compared to those in Band 10. In order to apply the SW algorithm, it is necessary to remove the obvious stripes while preserving the data information of Band 10 and Band 11 as much as possible. The conventional destriping method is not suitable for TIR data due to the oblique stripe. This study adopted the edge detection method to remove the stripes in Band 10 and Band 11. The main idea of this process is to determine the stripe boundary and the central position based on the image gradient so as to extract the stripe. The destriping of Band 10 and Band 11 includes 6 steps.

• Remove clouds and the land from images. Cloud and land detection are based on the Landsat Collection 1 Level-1 Quality Assessment (QA).
• Calculate the horizontal derivative using the Sobel operator. The Sobel operator first performs weighted smoothing and then engages in differentiation to enhance the edge mutation. It has the capabilities of noise suppression and accurate edge positioning. Since the stripes are vertical, and there are fine stripes and bright spots on the image, the horizontal gradient can be calculated by the vertical edge Sobel operator, as follows:

  $$\begin{gathered} G\left(i,j\right)=f\left(i+1,j-1\right)+2f\left(i+1,j\right)+f\left(i+1,j+1\right)-f\left(i-1,j-1\right)-2f\left(i- \\ 1,j\right)-f\left(i-1,j+1\right), \end{gathered}$$

  (1)

  where $f\left(i,j\right)$ is the gray value of image at position $i$ and $j$, and $G\left(i,j\right)$ is the horizontal gradient of $f\left(i,j\right).$
• Set threshold and extract the boundaries of the stripes. After repeated experiments, it is found that when the edge detection threshold is set to 27, the stripe edge in the image can be extracted completely.
• Determine the central position of the stripes. For each pixel, the gradients with respect to the two nearby pixels in a horizontal direction are calculated. Then, the average of the positive and negative horizontal gradients of all the pixels in the images are calculated, respectively. If the pixel located between the positive and negative horizontal gradients is larger than the corresponding average value, it is considered as the center of the stripe.
• Extracting stripes. If the center of the stripe exists between two boundaries within the certain range, the area between the two boundaries is the stripe.
• Fill in the stripes. When the cell is determined as the stripe, the pixel value is replaced with the average of the surrounding 5*5 window pixels.

2.4. Sea Surface Temperature Inversion Algorithm

In the 1960s, the United States began to study global SST quantitative inversion and its impacts on environmental changes. With the advent of AVHRR sensors in the 1970s, thermal infrared-based remote sensing technology has attracted increasing attention. The inversion algorithms have been investigated for the production of SST products using aerospace infrared radiometers. Many sophisticated infrared sensing devices and corresponding inversion algorithms have been developed. At present, common SST inversion algorithms include the RTM, the SC algorithm, the MW algorithm and the SW algorithm.

2.4.1. Radiation Transfer Equation Method

The thermal infrared radiation transfer theory assumes that the atmosphere is horizontally uniform and does not consider multiple scattering [22]. The blackbody radiation of the sea surface $B\left(\lambda ,T_S\right)$ is:

$$B\left(\lambda ,T_s\right)=\frac{L_{sensor}-L_{up}}{\tau \epsilon }-\frac{\left(1-\epsilon \right)L_{down}}{\epsilon }$$

(2)

where $T_s$ is SST, $\lambda$ is the wavelength of the band, $L_{sensor}$ is the spectral radiance received by the sensor (W*m⁻²*sr⁻¹*µm⁻¹), τ is the total transmittance of the atmospheric, ε is the surface emissivity, and $L_{up}$ and $L_{down}$ are the upwelling and downwelling atmospheric radiance respectively. The sea water ε in the Band 10 is 0.98 [23]. $\tau$, $L_{up}$ and $L_{down}$ can be obtained by MODTRAN simulation according to the atmospheric profile and surface meteorological data. For Band 10, these values were calculated using the NASA atmospheric correction parameter calculator (https://atmcorr.gsfc.nasa.gov/).

2.4.2. The Single Channel Algorithm

In 2003, Jiménez-Muñoz et al. [24] proposed a single channel method that only requires the atmospheric water vapor content. This algorithm performs a first-order Taylor expansion on a Planck function near a certain temperature value, and simulates the atmospheric effect by using a Gaussian triangular filter function as the thermal infrared band channel response function. The formula is

$$T_s=\gamma \left[{\epsilon }^{-1}\left({\psi }_1{L}_{sensor}+{\psi }_2\right)+{\psi }_3\right]+\delta$$

(3)

where $\gamma$ and $\delta$ are intermediate variables dependent on the Planck function, and ${\psi }_1$, ${\psi }_2,$ and ${\psi }_3$ are atmospheric functions. These parameters are given by:

$$\gamma ={\{\frac{c_2{L}_{sensor}}{T_{\mathrm{sensor}}^2}[\frac{{\lambda }_{\mathrm{eff}}^4}{c_1}{L}_{\mathrm{sensor}}+{\lambda }_{eff}^{-1}]\}}^{-1}$$

(4)

$$\delta =-\gamma L_{\mathrm{sensor}}+T_{\mathrm{sensor}},$$

(5)

$${\psi }_1=1/\tau ,{\psi }_2={-\mathrm{L}}_{\mathrm{down}}-L_{up}/\tau ,{\psi }_3=L_{\mathrm{down}}.$$

(6)

In Equation (4), $c_1$ = 1.19104*10⁸ W/m²/sr¹/µm, $c_2$ = 1.4388*10⁴ µm*K, and ${\lambda }_{eff}$ refers to the effective wavelength of the channel considered. For Band 10 of Landsat 8, the effective wavelength can be replaced by the central wavelength. The total atmospheric water vapor content ω comes from the MODIS near infrared water vapor product, MOD05. The function relationship between ${\psi }_1$, ${\psi }_2,$ and ${\psi }_3$ and the total atmospheric water vapor content $\omega$ is as follows:

$${\psi }_1=-0.0047{\omega }^3+0.1241{\omega }^2-0.3299\omega +1.3157$$

(7)

$${\psi }_2=0.0624{\omega }^3-1.4300{\omega }^2-2.7246\omega -2.5905,$$

(8)

$${\psi }_3=-0.0755{\omega }^3+0.8319{\omega }^2-1.1160\omega +0.9408.$$

(9)

2.4.3. The Mono Window Algorithm

In 2001, Qin et al. [25] proposed a surface temperature algorithm applicable to Landsat TM under the same conditions of upward and downward atmospheric radiance. This algorithm needs only three parameters: the surface emissivity ($\epsilon$), the atmospheric transmittance ($\tau$) and the effective mean atmospheric temperature (T_a). These parameters can be estimated based on in situ meteorological observation data, without using real-time atmospheric profile data. The algorithm is:

$$T_S=\left[a\left(1-C-D\right)+\left(b\left(1-C-D\right)+C+D\right)T_{sensor}-D{T}_a\right]/C$$

(10)

where $C, D$ are the intermediate variables, $C=\tau \epsilon ,D=\left(1-\tau \right)\left[1+\tau \left(1-\epsilon \right)\right]$, and $a$ and $b$ are the regression coefficients. Qin defined the temperature parameter L that is linear with temperature T; based on Band 10 of the Landsat-8 spectral response function and linearly fitting L and T to the temperature range of 273.15~313.15 K (as shown in equation (11)), a = –60.98, b = 0.4278, and correlation coefficient $R^2$ = 0.9998.

$$L=B\left(T\right)/\left[\partial B\left(T\right)/\partial T\right]$$

(11)

$$L=–60.98+0.4278T$$

(12)

$T_a$ is a parameter introduced to simplify the atmospheric radiation, which is related to the temperature of the atmospheric profiles. Qin et al. [25] proposed to estimate $T_a$ based on the near surface air temperature ($T_o$) of the four standard atmospheric profile models, as shown in

Table 3. The calculation model of the effective mean atmospheric temperature ($T_a$) in Equation (10).

Atmospheric Model	Ta Estimation Formula
Tropical	17.9769 + 0.91715To
Mid-latitude summer	16.0110 + 0.92621To
Mid-latitude winter	19.2704 + 0.91118To
Standard atmospheric	25.9396 + 0.88045To

.

2.4.4. The Split Window Algorithm

The SW algorithm was first proposed by Mcmilln [26] for NOAA/AVHRR sensors. The SW algorithm uses the atmospheric transmission difference of the AVHRR 4 and 5 channels to eliminate atmospheric influence. The split window method has been improved since then, producing a variety of algorithms. This study adopts the NLSST algorithm proposed by Walton [27], which has been widely used in MODIS, VIIRS, VIRR and other products. The NLSST solution is

$$T_s=a_1+a_2{T}_{11}+a_3{T}_{sfc}\left(T_{11}-T_{12}\right)+a_4\left(T_{11}-T_{12}\right)\left(\sec (\theta \right)-1)$$

(13)

where $T_{11}$ and $T_{12}$ are the brightness temperatures of the 11 µm and 12 µm channels, $a_i$ (i = 1~4) are the regression coefficients, $\theta$ is the zenith angle, and $T_{sfc}$ is a priori estimate of the SST. $T_{sfc}$ adjusts the difference between the $T_{11}$ and $T_{12}$. $T_{sfc}$ can refer to the climate SST or the satellite SST estimated by MCSST. The Landsat-8/TIRS is a narrow field sensor, and the maximal zenith angle is only 7.5°. Thus, we ignored the impact of the zenith angle. $\theta$ is regarded as 0. Thus, the SST can be calculated as follows:

$$T_s=a_1+a_2{T}_{11}+a_3{T}_{sfc}\left(T_{11}-T_{12}\right)$$

(14)

Since the collocations between Landsat and the buoys are too few to develop SST retrieval algorithms, this study developed an SW algorithm based on MODIS SST and validated it with the buoy SST. A total of 56 cloudless Landsat 8 images from January 2017 to January 2019 (after destriping) and the MODIS SST were selected to derive the coefficients $a_1$, $a_2$, and $a_3$ in Equation (12) (

Table 4. The regression analysis coefficients $a_i$ and the correlation coefficient.

	${\mathit{a}}_1$	${\mathit{a}}_2$	${\mathit{a}}_3$	R2
Spring	−18.4206	1.0619	0.0080	0.96
Summer	81.6599	0.7157	0.0080	0.66
Autumn	−0.6963	1.0013	0.0083	0.93
Winter	−33.3589	1.1156	0.0073	0.98

). The Landsat 8 brightness temperatures are averaged over 33*33 cells (about 1 km², close to the resolution of MODIS) in order to be comparable with the MODIS SST. More than 30,000 measurements are collocated in each season. The coefficients $a_1$, $a_2$, and $a_3$ are obtained via the linear regression.

3. Results

3.1. TIRS Image Destriping Results

Figure 5 and Figure 6 show the results of destriping in the open ocean and coastal areas. The results show that the Sobel edge detection method can completely extract the obvious stripes in the Band 11 without contaminating other valid information and can be adapted to open ocean and coastal areas.

As shown in Figure 7, the SST distribution changes after destriping are more uniform than those before destriping, which is more conducive to warm drain-off water studies. Because of the lack of buoy data for most of the striping area, the MODIS SST product is used to estimate the influence of the destriping processing. A total of 29,567 data samples were matched from January 2017 to January 2019 for the comparisons between MODIS SST and Landsat SST.

Table 5. Comparison of the SST inversion results before and after destriping at the stripe.

	Number	Bias (°C)	MAE (°C)	RMSE (°C)	STD (°C)
Destriping Before	29567	0.24	0.67	0.86	0.82
Destriping After	29567	0.16	0.56	0.72	0.70

shows the difference between the MODIS SST and the destriping Landsat SST. The bias reduces from 0.24 to 0.16 °C, and the standard deviation (STD) reduces from 0.82 to 0.70 °C. This indicates that the proposed destriping method can very effectively improve the stripe problem.

3.2. Validation of Sea Surface Temperature Algorithms

In this section, the buoy data and MODIS SST products are used to validate the SST inversion algorithms.

Table 6. Comparison of the four inversion methods and the buoy data (standard data).

	Number	Min (°C)	Max (°C)	Bias (°C)	MAE (°C)	RMSE (°C)	STD (°C)
RTM	31	−0.96	2.36	0.55	0.55	0.92	0.74
SC	31	−0.87	3.87	1.04	0.83	1.44	0.99
MW	31	−3.51	2.44	−0.82	1.05	1.81	1.61
SW	31	−1.70	1.37	−0.01	0.12	0.71	0.71

and Figure 8 show the validation results of the algorithms compared with the buoy data. Thirty-one buoy data were collocated with Landsat SST. The results demonstrate that the SST from the SW algorithm is significantly better than the SST from the RTM algorithm, the SC algorithm, or the MW algorithm. The bias from the SW algorithm is −0.01 °C, which is notably smaller than the values of 0.55, 1.04 and −0.82 °C from RTM, SC, and MW, respectively. Moreover, the $R^2$ is 0.98, and the RMSE is 0.71 °C in the SW inversion results.

In order to ensure the accuracy of the inversion results in the study area without buoy data, this study used the SST products (MODIS SST, ERA SST and OISST) from the Daya Bay region to validate the Landsat inversion results. The Landsat SST are averaged to follow the resolution of the ERA SST and OISST the same way as the MODIS SST.

Table 7. Comparison of the four inversion methods around the Daya Bay with MODIS SST.

	Number	Bias (°C)	MAE (°C)	RMSE (°C)	STD (°C)
RTM	13950	1.65	1.67	1.86	0.86
SC	13950	1.66	1.89	2.80	2.25
MW	13950	0.73	1.47	1.60	1.43
SW	13950	0.18	0.32	0.47	0.43

,

Table 8. Comparison of the four inversion methods around the Daya Bay with ERA SST.

	Number	Bias (°C)	MAE (°C)	RMSE (°C)	STD (°C)
RTM	594	0.58	1.44	1.78	1.68
SC	594	3.37	3.44	3.74	1.62
MW	594	−0.43	2.36	2.77	2.74
SW	594	−0.30	1.19	1.60	1.57

and

Table 9. Comparison of the four inversion methods around the Daya Bay with OISST.

	Number	Bias (°C)	MAE (°C)	RMSE (°C)	STD (°C)
RTM	158	1.50	1.76	2.11	1.78
SC	158	3.53	3.55	3.95	3.95
MW	158	1.38	1.88	2.35	1.91
SW	158	−0.50	0.87	1.13	1.02

outline a comparison between MODIS SST/ERA SST/OISST and the inversion results of each algorithm around the Daya Bay. The results show that the SST from the SW algorithm is more consistent with the MODIS SST, ERA SST and OISST than the other algorithms (Table 7, Table 8 and Table 9). Landsat 8 SST in the Daya Bay area provide more fine-scale details of the distribution of SST compared with other low-resolution SST products like MODIS, ERA and OISST.

In summary, the SST from the SW algorithm is the most consistent with the buoy, the MODIS, the ERA and the OISST SST among the four algorithms. Therefore, we use the SST from the SW algorithm for the subsequent research of the distribution of the warm drain-off from the two nuclear power stations.

3.3. Distribution of Warm Drain-Off Water

3.3.1. Characteristics of SST

A set of eight Landsat 8 SST images (two from each season from February 2017 to January 2019) illustrates the water temperature distribution in the Daya Bay (Figure 9a–f). The red pentagram in Figure 9 indicates the location of the DNPS and the LNPS. For a better display, we used color bars with different scales to show the SST.

Figure 9 shows the clear warm drain-off water, which was mainly distributed in the fan-shaped area from the nuclear power station to the center of the bay. The warm plume is about 1–4 °C warmer than the surrounding water, and the outlet plume can reach 5~6 °C. The influence radius of the warm plume can reach 3~6 km.

On August 29, 2017 and June 31, 2018, the plume was limited to a small area on the western side of the power station, where it extended six kilometers to the western part of the central archipelago on January 23, 2019 and February 18, 2017. This result may be related to the environmental characteristics of the Daya Bay. The Daya Bay is a weak tidal field; most of its water originates from the South China Sea, and its water exchange conditions are poor. Therefore, the huge energy brought into the seawater by the warm drain-off water is mainly produced through the heat exchange between the sea surface and the atmosphere. The heat transfer coefficient in the summer is larger than that in the winter, so the temperature increase area in the summer is smaller than that in the winter [28]. A study on water cooling engineering during the preliminary design stage of Guangdong LNPS calculated the heat transfer coefficient of the Daya Bay under average hydrometeorological conditions. Without considering the influence of the adverse wind direction, the heat transfer coefficient of the water surface in the summer is 47 W/(m²·°C); the heat transfer coefficient in the winter is 28.5 W/(m²·°C). Considering the influence of the adverse wind direction, the heat transfer coefficient of the water surface in the summer is 57 W/(m²·°C); the heat transfer coefficient in the winter is 37.2 W/(m²·°C).

3.3.2. Distribution of SST Rising Intensity

To reveal the specific intensity and range of the SST increase caused by warm drain-off water, the SST that was not affected by warm drain-off water from DNPS and LNPS was computed as the background temperature. As a semi-enclosed bay, the background temperature of the Daya Bay can refer to the average temperature of the entire Bay area. However, the two sewage outfalls of the Petrochemical New Town and the Pinghai Power Station in this bay also increased the surface temperature of the inner sea (as shown in Figure 1). Figure 9 shows that there is a temperature difference between the Eastern and the western Daya Bay due to a stray light. Therefore, the fan-shaped area with a 15 km radius from the nuclear power stations is taken as the study area. We calculated the average of the SST in the study area and then excluded the area where the SST is 1 °C higher than the average temperature. Finally, we counted the average SST of the remaining area as the background temperature (

Table 10. The Background temperatures, tides and wind information.

Date	Background Temperature (°C)	Tide	Tidal Current Speed (cm/s)	Wind Direction	Wind Speed (m/s)
2017.02.18	19.27	Flood tide	5.19	Northeast	1.81
2017.04.07	25.30	High slack	5.65	Southeast	0.41
2017.08.29	28.64	Flood tide	5.24	Northeast	1.92
2017.11.01	26.98	Ebb tide	7.18	Northeast	4.67
2018.03.09	19.25	Flood tide	6.28	Northeast	6.71
2018.07.31	32.19	High slack	6.54	Southwest	1.65
2018.10.03	29.79	Low slack	4.21	Northeast	3.98
2019.01.23	19.16	Flood tide	5.26	Northeast	3.09

). The warm plume effected area was divided into seven levels, as shown in

Table 11. Grade scale of various temperature increase intensities.

Temperature Rising Range	Levels
< 1 °C	< 1 °C
(1 °C, 2 °C)	+ 1 °C
(2 °C, 3 °C)	+ 2 °C
(3 °C, 4 °C)	+ 3 °C
(4 °C, 5 °C)	+ 4 °C
(5 °C, 6 °C)	+ 5 °C
>6 °C	+ 6 °C

, and annotated by different color.

Figure 10 shows the spatial pattern of the SST increase intensity around the nuclear power stations. In the process of warm drain-off water diffusion, the high temperature increase area (+4, +5, and +6 °C) is relatively small, extending outward in the circle. The low temperature increase area (+1, +2, and +3 °C) is relatively large, spreading around the drain-off outlet with different forms. This appearance may be due to tides and wind. The area affected by warm drain-off water shows a fan shape, and the temperature increase intensity decreases from the outfall to the northeast. The tidal current and wind field shown in Figure 10 promote or suppress the diffusion of warm drain-off water and even control the diffusion direction.

On February 18, 2017, August 29, 2017, March 9, 2018 and January 23, 2019, when the satellite passed over, the tide was a flood tide. Figure 10a,c,e, and h shows that the tidal current diverged near the nuclear power station and the warm plume drained into the eastern branch of the tidal current at Changshanwei. The warm drain-off water extended northeastward, and the warm drain-off water distributed across a long and narrow area due to the flood tide. On April 7, 2017 and July 31, 2018, the tide was a high slack tide, during which the flood tide became an ebb tide. Figure 10b,f shows that the tidal current changes from divergence to convergence. The warm drain-off water first extended northeastward, and then it turned southward when the tide flowed offshore. At the same time, the decrease in the current velocity in the southeastern direction led to the expansion of the small branch separated by warm drain-off water to the southeast. On October 3, 2018, the tide was a low slack tide at the time of satellite pass. During that day, the ebb tide became a flood tide. Figure 10g shows that the tidal current changes from convergence to divergence. The warm drain-off water expanded southward, but turned northward due to the slack tide. On November 1, 2017, the tide was an ebb tide at the time of satellite pass. The current converges near the nuclear power station. Since the tide inhibits the spread of the current to the northeast, the warm drain-off water flowed in a southwestern direction as the tide flowed toward the bay mouth, forming a fan-shaped area (Figure 10d).

Figure 11 shows that the maximum area affected by the warm drain-off water is about 16.89 km². The maximum warm area exceeding the background temperature of 4 °C is about 1.08 km², which is mainly concentrated in the drain-off channel. Therefore, the warm drain-off water of the nuclear power station will not affect the surrounding environment.

Based on Figure 10 and Figure 11, we find that the characteristics of the warm drain-off water diffuse with tidal states. During the flood tide, the warm drain-off water is moved to the northeast by the tide, the high temperature increase zone moves from the drain to the northeast, and the combined drain-off channel temperature increase intensity is +3~4 °C, which can be seen in the channel sections (Figure 10a,c,e,h; green or orange color). The +4 °C area increases, and the +6 °C area decreases or even disappears (Figure 11). At the same time, the water body is subjected to the action of the tide in the direction of the shore, resulting in a small area of a low temperature rising zone. On the contrary, during the ebb tide, the warm drain-off water is affected by the convergence of the current. The high temperature increase zone is very small and is mainly concentrated in the drain (Figure 10d, red color). The water body is dragged by the tide to the bay mouth, resulting in a large area for the low temperature increase zone. In summary, the warm drain-off water diffusion shows a distinct change with the tide state, which is different under different tidal conditions and relatively similar under the same tidal conditions.

In addition, wind may also affect the diffusion of warm drain-off water from the nuclear power station. Table 10 shows the wind field information. In the study area, the northeast monsoon is strong in winter, while the weaker southwest monsoon is dominant in summer. When the southwest wind prevails, the wind blows warm drain-off water away from the shore, which promotes the diffusion of warm drain-off water. When the northeast wind prevails, the wind forces warm water to flow toward the shore which is not conducive to the diffusion of warm drain-off water. Meanwhile, waters are vertically mixed in the Daya Bay under the influence of the strong northeast monsoon. Therefore, wind has little influence on the sea surface distribution of warm drain-off water. Therefore, the shape and distribution characteristics of the temperature increase area are mainly determined by tide, not wind.

4. Discussion

The SST retrieved from TIR data is mainly affected by factors such as the inversion algorithm, spectral radiance received by the sensor and the surrounding environment.

Common thermal SST inversion algorithms include the RTM, the SC algorithm, the MW algorithm and the SW algorithm. The accuracy of the inversion algorithm is relevant to the accuracy of the water vapor attenuated by the algorithm and the accuracy of the input parameters. The RTM is a priori the best option for SST retrieval from one thermal band, since it does not involve additional approximations [29]. However, this method requires accurate atmospheric parameters such as atmospheric transmittance and atmospheric upwelling and downwelling radiances to eliminate atmospheric effects, which is not always possible. The SC algorithm relies on atmospheric functions, which were assumed to be dependent only on atmospheric water vapor content [29]. Meanwhile, previous studies [6,21] have indicated that the variation of the SST retrieved from the single channel algorithm is inversely proportional to that of the atmospheric water vapor content, and the inversion error is proportional to the atmospheric water vapor content. With the assumption that the atmospheric upwelling and downwelling radiation values are the same, the MW algorithm introduces effective atmospheric temperature to remove atmospheric effects. It needs three parameters as input: atmospheric transmittance, surface emissivity and effective atmospheric temperature. These three algorithms need accurate atmospheric parameters as input, which are difficult to obtain, so the resulting error is large. The SW algorithm is a more mature and accurate temperature inversion algorithm than the single channel algorithm [30]. Li [31] indicated that the accuracy of the SW algorithm was invariant in high atmospheric water vapor. The SW algorithm uses an atmospheric transmission difference of 10~12 µm channels to eliminate atmospheric influence. Empirical regression only needs a large amount of fitting data to ensure high accuracy without any meteorological data, which has been widely used in MODIS, VIIRS, VIRR and other products.

The noise of the thermal IRS leads to errors in the spectral radiance received by the sensor, which influences the accuracy of the SST retrieved from the remote sensing data. The absolute radiometric calibration errors in the Landsat 8/TIRS Band 11 are two-times higher than those of the Landsat 8/TIRS Band 10. Yu et al. [12] found that there is a general overestimation trend for the LST estimated from Band 11 compared to the LST estimated from Band 10. Based on the Landsat 8/TIRS without a stray light correction, Cahyono et al. [32] used the SW algorithm to estimate the SST in coastal area and validated their results with in situ data. The result has a low accuracy with a bias larger than 1 K. After a stray light correction, Aaron [13] indicated that the temperature errors of over 6 K in the water before stray light corrections were reduced to 0–1 K for most of the scene.

The SSTs inversed from the infrared instruments have an effective measurement depth of less than 10 µm and represent the so-called surface skin temperature. The buoy SST is the bulk temperature at a depth of about 0.5 to 1.5 m. The skin water is usually cooler (0.1–0.3 °C on average) than the water one millimeter or more below the surface [33]. Meanwhile, this water is sensitive to the surrounding environment, such as sunlight and wind fields. Under low wind (<6 m/s) conditions, the sunlight warms the surface sea water, forming a warm layer within a few meters of the surface. Under high wind (>6 m/s) conditions, strong vertical mixing keeps the upper seawater from becoming the same temperature, but the skin temperature is still lower than the bulk temperature due to the heat loss caused by sensible and latent heat fluxes [34]. Moreover, sea water in the coastal zone exchanges violently on the beach through tidal movements and is gradually heated by solar radiation [35]. This study subtracts the background temperature, which reduces the above-mentioned SST differences caused by water depth and wind to discuss the spatial distribution of SST increases in the Daya Bay. More data are needed to study the effects of water depth and wind on the SST in the bay.

5. Conclusions

The development of Landsat 8 TIR provides an effective and high spatial resolution data source for the dynamic monitoring of warm drain-off water from nuclear power stations. This study is the first time attempt to apply the fine-scale SST inversion algorithm together with the edge detection destriping method to Landsat 8/TIR data to analyzing warm drain-off water from the two nuclear power stations in the Daya Bay. We use the edge detection method to remove stripes, which reduces the bias of the SST from 0.24 to 0.16 °C and the std of the SST differences from 0.82 to 0.7 °C compared with buoys. However, the stripe removal algorithm can only deal with stripes that have a large variation and a narrow width and cannot distinguish artifacts with large widths, which need further study.

After destriping, we compare and estimate the accuracy of the RTM, the SC algorithm, the MW algorithm, and the SW algorithm in the South China Sea. The results from the SW algorithm are very consistent with the results of the buoy, the MODIS SST, ERA SST and OISST among the four algorithms. The bias between the inversed SST and the buoy SST is −0.01 °C and the RMSE is 0.71 °C. Based on the SST results retrieved by the SW algorithm, the distribution characteristics of warm drain-off water from the DNPS and LNPS are analyzed. The warm drain-off water is about 1 to 4 °C higher than the surrounding water, and the outlet plume can reach over 5 to 6 °C. The influence radius of the warm plume can reach 3 to 6 km. Based on the SST invention results, the FVCOM tidal currents and ERA-Interim winds, the tide is shown to determine the shape and distribution characteristics of the temperature increase area and to control the diffusion direction. Notably, when the tide rises, the warm drain-off water flows to the northeast and extend over a long and narrow area. When the tide goes out, the warm plume flows to the bay mouth under the pull of the tide, forming a fan-shaped area. Moreover, the temperature increase intensity in the combined discharge channel during the flood tide is lower than that during the ebb tide and the +1 °C area during the flood tide is smaller than that during the ebb tide.

In general, our results are consistent with those of previous studies. The new SW algorithm together with the edge detection destriping method improves the accuracy of Landsat 8 SST in the Daya Bay area and provides greater fine-scale details on the distribution of SST compared with other low-resolution SST products like OISST and ERA SST.