Two field stations, S1 and S2, were selected to conduct field work ( Figure 1D ). These two field stations were arranged along the coast of Hengshan Island ( Figure 1D ), which is opposite to the Qiangjiao Power Plant ( Figure 1C ). Thermal discharge from the Ninghai Power Plant can affect seawater near Hengshan Island. The seabed was sandy with gravel at station S1 and muddy at station S2.

The field work was conducted between August 31 and September 15, 2019, covering an entire spring-neap tidal cycle. This period was in late autumn, which facilitated the capture of the thermal discharge. The field work was conducted in the wet season, and the signal of super Typhoon Lingling (no. 1913) was captured ( Figure 1F ), which impacted the bay (23−28°C, moderate rain) around September 6, 2019. On August 30, 2019, Lingling (tropical disturbance) occurred near the island of New Guinea, and was upgraded to a tropical depression on September 2, after which the typhoon began to intensify rapidly. The storm was upgraded to a severe tropical storm on September 3. Subsequently, it was upgraded to a super-typhoon. The maximum wind speed exceeded 50 m/s and the central pressure was 900 hPa when Typhoon Lingling passed by the east coast of China (over 2000 km away, in the open East China Sea). From September 5 to 6, 2019, moderate to heavy rains were observed in the central and northern coastal areas of Zhejiang. During the typhoon transit, the rains in the coastal areas of Zhejiang were 30–70 mm, with some over 100 mm.

The instruments were arranged on two frames, and both of them were 0.5 × 0.5 × 1.5 m. Every fame had a tidal gauge, an optical backscatter sensor (OBS), and an acoustic doppler velocimeter (ADV). The probes of the ADV were set 15 cm above the sea bed. The tidal gauge and the OBS were set at the same level above the seabed. A tidal gauge was set to measure the data every minute at a frequency of 4 Hz. The OBS and ADV were set to operate every 10 min. The tidal flat temperature was measured using temperature sensors near S2 station ( Figure 1D ) in the tidal flat. Temperature sensors were used to measure the temperature of the tidal flat every 10 min. Three sensors were arranged vertically, that is, near the surface of the tidal flat, 20 cm below the surface, and 82 cm below the surface ( Figure 1E ).

The tidal gauge data were directly extracted and temporally averaged every 10 min to obtain tidal level data. Wave height data were obtained from the tidal gauge fluctuation data. The wave parameters were then calculated.

The OBS data were calibrated with laboratory experiments before the field work. The correlation of the turbidity with the suspended sediment concentration (SSC) was high, indicated by the R² value of 0.92 and 0.96 at station S1 and S2, respectively ( Figure 2 ). The turbidity data were then converted to suspended sediment concentration values.

The ADV data were first filtered by the ratio of the signal to noise (SNR) and the correlation of the acoustic signals. Data with SNR values lower than 5% or correlation values lower than 70% were excluded from the series. The phase–space thresholding method (Goring and Nikora, 2002) was then used to remove singular points from the data. A method similar to the Improved Phase-Space Thresholding Method (Wahl, 2003; Parsheh et al., 2010) was also used to calibrate the accuracy of the results.

The turbulence data were obtained using the differences between the instantaneous current velocity and the averaged velocity. Subsequently, the turbulent current and turbulent kinetic energy were calculated.

To determine whether the thermal discharge from the Ninghai Power Plant affected the detectors, we used Landsat 8 data for temperature inversion (http://ids.ceode.ac.cn/query.html, Chen et al., 2018). Landsat 8 includes an Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS). The OLI land imager includes nine bands with a spatial resolution of 30 m, including a 15-meter panchromatic band with an imaging width of 185 × 185 km. To avoid atmospheric absorption features, OLI readjusts the bands as follows:

where T_s is sea surface temperature, a and b are empirical coefficients.

Landsat-8/TIRS has two thermal infrared channels, and sea surface temperature (SST) remote sensing inversion can use both single- and dual-channel imaging data. The single-window algorithm and the universal single-channel algorithm proposed by Jiménez-Muñoz and Sobrino (2003) were used in this study.

Field stations S1 and S2 were located in the tidal flat of Hengshan Island, which is opposite to the Ninghai Power Plant. Therefore, the sea surface temperature (SST) near the Ninghai Power Plant was retrieved using remote-sensing data ( Figure 3 ). SST data for summer (on July 29, 2019) and winter (on December 2, 2019) were retrieved at the head of the bay. The thermal discharge of the power plant diffused to stations S1 and S2 in the tidal flat of the Hengshan Island. The influence of thermal discharge on the water temperature in winter was more obvious ( Figure 3B ) than that in summer ( Figure 3A ). The temperature near the power plant was approximately 8°C higher than that in the vicinity.

Macro-tidal and semi-diurnal signals during the spring-neap tidal cycle were well captured during the field work ( Figure 4A ). The tidal ranges at the two stations were nearly 6 m during spring tides and approximately 1 m during neap tides. The storm surge signal was weak and almost undetectable from the water level data because the typhoon was more than 2000 km away from the bay, and the bay was sheltered by its narrow and semi-enclosed morphology. Double peaks during low-slack waters were also captured. According to the results of the harmonic analysis, the interaction of the M2 and M4 tides contributed to the double peaks of sea surface elevation.

The currents were weak at the two stations ( Figure 4B ), with a maximum current speed of approximately 0.5 m/s. The current magnitude and direction at station S1 fluctuated and were disordered. Around September 6, the current speeds at station S1 were high. This is most likely due to the impacts of the Typhoon Lingling. The maximum wind speed near the station during this period exceeded 15m/s, with a northerly wind most of the time. The currents at station S2 were stronger than those at station S1. Current magnitudes at station S2 showed spring-neap and flood-ebb tidal signals, with strong currents during spring tides and weak currents during neap tides. The current directions at station S2 were mainly along the coast in the northwest-southeast direction.

Water temperature at both stations showed an increasing trend during the field work and a periodic pattern in a similar period with the tides ( Figure 4C ). The sea water temperature fluctuated during the field work, with occasional low temperatures during the first several days and occasional high temperatures during the last several days. This fluctuation was due to the exposure of the temperature sensor to air at this moment, as seen from the zero depth in Figure 4A . The maximum temperature of the seawater at stations S1 and S2 was near 34.5°C, and the minimum temperature was approximately 21°C. The maximum and minimum values of seawater temperature occurred at peak ebb during spring tides when the tidal flats were exposed. Typhoons occurred around the neap tides (red box in Figure 4 ). When sensors are in seawater, their temperature is affected by the seawater temperature. The water temperature remained near 30°C during most of the field work. However, the temperature fluctuated with large amplitudes from approximately 22°C to 34.5°C during the following spring tides. During the neap tides, the water temperature was stable and remained at approximately 29°C. At station S1 ( Figure 4C ), the water was occasionally cool (22°C) during high waters and warm during low waters before the typhoon. After the typhoon, the water was warm (>28°C) throughout the observation period. The peak values of water temperature occurred periodically during the slack waters and showed an increasing trend. A similar trend was observed for water temperature at station S2 ( Figure 4C ). Different patterns of water temperature occurred during the first and second spring tides at both the stations.

The suspended sediment concentration (SSC) at station S2 (muddy seabed) was highly correlated with the tidal dynamics, with large SSC values (~0.9 kg/m³) during spring tides and small values (~0.1 kg/m³) during neap tides ( Figure 4D ). The SSC values peaked after the high slack waters and were directly affected by the current speed. During spring tides, the SSC values ranged from approximately 0.2 kg/m³ to about 0.9 kg/m³, whereas during neap tides, the SSC values were around ~0.2 kg/m³. SSC values at station S1 (gravel seabed) were much smaller (mostly<0.2 kg/m³) than those at station S2, and showed a fluctuating pattern ( Figure 4D ). The SSC values were generally smaller than 0.2 kg/m³ through the field work.

Tidal current statistics ( Figure 5 ) during the spring-neap tidal cycle showed that currents at station S1 ( Figure 5A ) were largely in the southwest direction. At station S2, the currents were in the northwest-southeast direction most of the time ( Figure 5C ). During the spring tides and neap tides, the current directions had similar distributions to those during the spring neap tidal cycles.

The residual currents at S1 were 0.01 m/s (181 degree), 0.02 m/s (209 degree), and 0.02 m/s (196 degrees) during the spring tide, neap tide, and spring-neap tidal cycle, respectively. Station S2 had smaller residual currents of less than 0.01 m/s during the spring tide, neap tide, and spring-neap tidal cycle, compared to those at station S1.

The wave parameters calculated using the ADV data are listed in Table 1 . The maximum wave height was 0.13 m and 0.14 m and the maximum wave period was about 9.1 s and 11.8 s at station S1 and S2, respectively. The mean wave height was about 0.02 m and mean wave period was about 2.5 s at both of the two stations ( Table 1 ).

The time series of Hs in Figure 6A shows that the significant wave heights at station S1 (blue lines) were slightly smaller than those at station S2 (red lines). The Hs values during the typhoon (around September 5–6, red box) were approximately double (~0.06 m) those at the other time (~0.03 m) at both stations. Large Hs values occasionally occurred at S2.

The average wave period was mainly distributed between 1 s and 3 s ( Figure 6B ). The two sites were located at the head of the bay, where the impacts of winds and swells are small owing to the topography. The significant wave height was small, approximately 0.03 m during calm weather, with a wave period of 1.5~2 s. During the typhoon, the maximum significant wave height was about 0.1 m ( Figure 6A ).

The distributions of the wave periods and normalized wave heights (H/Hrms) at S1 are shown in Figure 6C . Hrms is the root mean square of wave height H, which indicates the mean wave energy. The histogram of the wave periods ( Figure 6B ) indicates a skewed wave period distribution, with a peak at 2 s. The normalized wave heights ( Figure 6C ) show that the wave energy was asymmetrically distributed, with a peak near H/Hrms of 0.6.

The turbulence intensity (I) and turbulent kinetic energy (TKE) were calculated as follows:

where u ¯ is the averaged current velocity, u’_rms is the room mean square of the turbulent current velocity, and u’, v’, and w’ are the turbulent current velocities in the x, y, and z directions, respectively.

As shown in Figure 7 , the turbulent intensity at stations S1 and S2 was high, as indicated by the high values of I. The turbulent TKE at station S1 ( Figure 7A ) was higher than that at station S2 ( Figure 7B ), and the TKE values at station S1 were largely double those at station S2. The TKE magnitudes showed obvious spring-neap signals and flood-ebb variations at station S2, with higher TKE during spring tides and lower TKE during neap tides. There were small spring-neap variations in TKE at station S1.

Utilizing the power spectrum analysis, the properties of the velocity fluctuation (u'/v'/w') in the frequency domain wer eexamined (Li et al., 2022; Chang et al., 2019). As shown in Figures 7E , F , the typhoon increased both the horizontal and vertical energy. In burst 98 (Hs = 0.05 m, T = 2.75 s, h = 2.81 m), the horizontal energy in the wave frequency range (0.04-0.5 Hz) was greater ( Figure 7E , green and red lines), demonstrating that the horizontal wave orbital velocities predominated the seawater turbulence. During the typhoon (e.g., burst 694, Hs = 0.10 m, T = 5.42 s, h = 2.30 m), the energy of the horizontal and vertical velocity fluctuations and their corresponding peak periods increased ( Figure 7C, D ).

Waves in shallow-water areas flowed elliptically. The horizontal and vertical velocity components of the waves are given by:

where H, σ, k, and h represent the wave height, angular frequency, wavenumber, and water depth, respectively. Coshkh and sinhkh were simplified to 1 and kh, respectively, for shallow-water waves. The wave ellipses grow flatter toward the bottom, indicating that the vertical orbital velocities degenerated faster than the horizontal ones (Dean and Dalrymple, 1991; Bian et al., 2020; Li et al., 2022). As a result, the horizontal energy (red and green lines in Figures 7E , F ) was greater than the vertical energy in the wave frequency range (blue lines in Figures 7E , F ).

The temperature variation in the tidal flat was largest near the seabed surface ( Figure 8 , red lines), followed by the second layer (20 cm below the tidal flat surface, aqua blue). The third layer had the smallest variation in temperature (82 cm below the tidal flat surface, sky blue). This is related to the more frequent heat exchange between the tidal flat surface, water, and the atmosphere. Peaks (low or high) in temperature occasionally occurred. The lowest temperature was approximately 15°C, and the highest temperature was approximately 35°C. The highest and lowest temperatures occurred near the surface of the tidal flat (red lines in Figure 8 ). The temperature peaks might have been caused by the air temperature because the sensor was located near the seabed surface. Simultaneously, the sensor temperature was also affected by the tidal flat. The specific heat capacity of the tidal flat was less than that of the water body; therefore, the tidal flat temperature increased faster than the sea water temperature.

As seen from the time series of the temperature data and tide level data ( Figure 9 ), the temperature fluctuation was small (mostly within 1°C) during the neap tides and large (>4°C) during the spring tides. The water temperature fluctuated around 29−30°C during the field period. Temperature peaks/troughs occurred during low slack waters during the spring tides, with negative troughs occurring during the first spring tidal period (September 1−6) and positive peaks occurring during the second spring tidal period (September 9−14). In the first spring tidal period, temperature troughs occurred in the low slack waters during the same periods of exposure of the tidal flat. In the second spring tidal period, temperature peaks occurred in low-slack waters. During the neap tides, temperature peaks and troughs occurred during low slack water periods and high slack water periods, respectively.

Water temperature was affected by the combined effects of the ocean and air. During neap tides, the water depth changes were small, and the water temperature was less affected by the atmosphere. During the spring tides, the water depth was shallow in the spring low-slack waters, and the water temperature was sensitive to the atmosphere and tidal flat. This indicates that the temperature fluctuated with the ocean and air temperatures. In the low slack waters of the first spring tidal period, the sensor was exposed to air and showed air temperature instead of water temperature. When the water depth is shallow, the seawater temperature is sensitive to the atmosphere. As Typhoon Lingling impacted the bay around September 6, 2019 (during the first spring tidal periods), the air temperature was low (22°C on September 2), compared with that during the second spring tidal cycle (32°C on September 3).

Figure 10 shows the time variation of water flux per unit width at the two stations. At station S1 ( Figure 10A ), the water fluxes were mainly along the shore in the west-east and south-north directions, with large eastward and southward fluxes, as indicated by the large values in the x direction and large negative values in the y direction. In addition, during the typhoon, the water fluxes at station S1 shifted from west to east. At station S2 ( Figure 10B ), the water fluxes were large in the east and north directions during the field work. The net water fluxes at station S1 were mainly seawards. The net water flux at station S2 was mainly landward.

Comparing the heat fluxes in the unit widths of the two stations ( Figure 11 ), the changes in the heat fluxes at the two stations were similar. Large or small variations occurred during the spring/neap tides. These changes may be due to the shallow tidal level and large temperature changes at that time. During the typhoon, the change in seawater temperature decreased, and the change in heat flux was small.

The sea water temperature fluctuated at stations S1 and S2. To study the variation and mechanism of sea temperature during the observation period, we performed the empirical mode decomposition (EMD) analysis on the temperature data ( Figure 12 ). EMD is a signal processing method (Huang et al., 1998; Huang et al., 1999) that decomposes a signal according to the time scale characteristics of the data itself. It has advantages in dealing with nonstationary and nonlinear data, and decomposes a complex signal into a limited number of eigenmode functions (IMFs).

The power spectra of the seawater temperature are shown in Figure 12A . Seawater temperature was decomposed using the EMD method. The instantaneous frequency time series ( Figure 12B, C ) and power spectrum analysis of the IMFs ( Figures 12B', C' ) were obtained by calculation. Small-scale features (>0.5 Hz) (e.g., turbulence) can be detected at instantaneous frequencies of IMF 1 and IMF 2. IMF 3 represents relatively low-frequency motion. The peak frequencies of IMF 1 at both stations S1 and S2 were mostly in the turbulent frequency range (>0.5 Hz). IMF 2 had substantial turbulence information at its peak frequencies, which were in the wave frequency range (0.04–0.5 Hz). IMF 1 and IMF 2 were therefore classified as the turbulent components.

Each IMF component represents the change in different scales, and each IMF component is a narrow-band signal ( Figure 12D ). IMF 1 had the largest amplitude and highest frequency. The frequency decreased in IMF 2 and 3, and the amplitude also decreased. The temperature fluctuation alternated between large and small fluctuations, which was the same trend as that observed for the spring-neap tidal cycle. When decomposed into five IMF components, the overall trend of temperature increased, as shown by the residual ( Figure 12D ). This indicates that the seawater temperature at both stations increased slightly during the field work ( Figure 12D ).

The tidal flat temperature of each sensor was decomposed into a series of IMFs using Hilbert transform (Huang et al., 1998). The instantaneous frequency time series and power spectrum analysis of the IMFs decomposed using EMD are shown in Figure 13 . Sensor 1 in the surface layer exhibited the largest variation in amplitude and showed a high degree of agreement with the original data ( Figure 13A ). Sensor 3 exhibited the smallest amplitude variations. Sensors 2 and 3 had similar trends and energies in the turbulent frequency band (>0.5 Hz), as both sensors were in the subsurface layer of the tidal flat.

Each IMF component had a relatively stable quasi-period, and the variance contribution was used to express the effects of the amplitude and frequency on the raw signal data. The IMFs of the seawater temperature and surface tidal flat temperature were similar, but the IMFs of the middle (20 cm below the seabed) and bottom (82 cm below the seabed) layers in the tidal flat were less correlated. This is because both the sea water temperature and surface tidal flat temperature are subject to similar sea-air-tidal flat interactions.

The variance contribution of the IMF 3 component was large, and the average period of the IMF 3 component was close to 24 h, which was mainly due to the daily variation in weather and influenced by the heat flux cycle ( Figure 13E ). The average period of IMF 2 was close to 12.42 h, which was mainly due to the tidal effect caused by the irregular semi-diurnal tide characteristics in the bay. Therefore, the temperature fluctuations of seawater and tidal flats were mainly affected by air temperature and tides. Thermal discharges spread with tide, with higher temperatures during ebb tides and lower temperatures during flood tides ( Figure 9 ), as the power plant was located upstream of the field stations ( Figure 3B ). Therefore, the second most influential factor of the temperature change (12.42 h) may be caused by thermal discharges. The average period of IMF 1 was less than 4 h, which may be related to the inundation of the tidal flats or exposure. The peak frequencies of IMF 1 for tidal flat temperature decomposition were mostly in the turbulent frequency range (>0.5 Hz), which may thus be subject to waves. The overall trend of the tidal flat temperature decreased during field work, as shown by the residuals ( Figure 13E ).