Encountering shoaling internal waves on the dispersal pathway of the pearl river plume in summer

Fundamentally, river plume dynamics are controlled by the buoyancy due to river effluent and mixing induced by local forcing such as winds and tides. Rarely the influence of far-field internal waves on the river plume dynamics is documented. Our 5-day fix-point measurements and underway acoustic profiling identified hydrodynamic processes on the dispersal pathway of the Pearl River plume. The river plume dispersal was driven by the SW monsoon winds that induced the intrusion of cold water near the bottom. The river effluent occupied the surface water, creating strong stratification and showing on-offshore variability due to tidal fluctuations. However, intermittent disruptions weakened stratification due to wind mixing and perturbations by nonlinear internal waves (NIWs) from the northern South China Sea (NSCS). During events of NIW encounter, significant drawdowns of the river plume up to 20 m occurred. The EOF deciphers and ranks the contributions of abovementioned processes: (1) the stratification/mixing coupled by wind-driven plume water and NIWs disruptions (81.7%); (2) the variation caused by tidal modulation (6.9%); and (3) the cold water intrusion induced by summer monsoon winds (5.1%). Our findings further improve the understanding of the Pearl River plume dynamics influenced by the NIWs from the NSCS.


Results
Meteorological variability. Winds measured at both sites were mostly southerly/southwesterly (Supplementary Figs. 1b, 2a), which were occasionally interrupted by northerly/easterly winds coinciding with lower atmospheric pressure. At ZHJ2, a stormy weather event occurred between 18:00 July 25 and 11:00 July 27 in which the gust wind speed reached 13.7 m/s. This suggests the effect of the peripheral circulation of Tropical Storm Mirinae or the passing tropical rain belt 46 . The average wind speed was 4.3 m/s during the experiment at ZHJ2. Spatial variability of the coastal currents. The progressive vector (PV) plots of instantaneous currents at both sites show the net NE displacements ( Supplementary Fig. 1c, 2b,c), which are the typical summer current patterns seaward of the PRE 37,41,45 . At ZHJ1, the non-tidal current (at 16.6 m) was largely flowing north-northeastward, along the Dangan Channel ( Supplementary Fig. 1). The form number (F) indicates the tidal current was mainly composed of the semi-diurnal tide (Table 1). Superimposed by the M 2 and K 1 tidal components on the non-tidal flow ( Supplementary Fig. 1c,d and Table. 1), the instantaneous current was in a zigzag pattern in northwest-southeast directions (on-offshore). However, the tidal energy ratio (ER) was 38.9%, suggesting that the tide was not the primary forcing to drive the instantaneous flow for the short observational period (Table 1). Near the end of the observation, the tidal signal in the instantaneous current became weak, and the current direction remained consistently northeastward. The average speed of the instantaneous current in the measurement was 0.18 m/s. At ZHJ2, the 5-day current observation showed the PV displacements of instantaneous currents rapidly decreased with depth, indicating the slower current speed in the lower water column ( Fig. 2b and Supplementary  Fig. 2). The average instantaneous current speed was 0.39 m/s at the surface, but only 0.12 m/s near the bottom. The differences between the surface and bottom flows were found not only in the current speed, but also in the current direction. Because of the interaction among the large-scale current 47 , bottom frictional effects 48 , and wind-induced Ekman pumping 47 , the current structure in the vertical mainly turned counterclockwise as the water depth increased ( Fig. 2b and Supplementary Fig. 3). The instantaneous flow to ENE in the upper water column, and then the flow veered toward NE near the bottom layer. Additionally, the non-tidal bottom current showed stagnation between 10:00 July 26th and 10:00 July 27th (marked with a red star in Fig. 2b,c) though the surface current flowed consistently northeastward. www.nature.com/scientificreports/ In general, the instantaneous flow at each depth was composed of a northeastward non-tidal component superimposed by a short circular tidal component (Fig. 2b,d). The tidal flow was basically composed of NE-SW orientated diurnal and NW-SE semi-diurnal tidal components (Fig. 2e). The F value indicates the tidal current was mixed with diurnal dominance in the upper water column and mixed with semi-diurnal dominance in the lower water column (Table 1). Overall, the K 1 constituent was dominant at surface, and M 2 became more important in the lower water column where the ER was over 40%. However, the tidal current at sub-surface layers (7 m) showed a significant increase of K 1 strength ( Fig. 2e and Table 1). This variability also showed in the F and ER value as a 'sandwich structure' above 21 m, in which an abrupt increase at 7 m and decrease at 14 m. This might be due to the NIWs perturbation in the current field.
To understand the tidal motions in the upper part of the water column, the stick diagram of vectorially averaged tidal velocity above the 5-m depth is plotted against the fluctuations of water elevation measured by the ADCP on the mooring (Fig. 3a). At the mooring location, because of diurnal inequality, tidal flows were landward during one (or two) consecutive rising tides (from lower low-water to higher highwater, marked "R") and seaward during one (or two) consecutive falling tide (from higher highwater to lower low-water, marked "F"). Each rising and falling cycle included two floods and ebbs. The tidal range was around 2-m, which is defined as a micro-to meso-tidal environment 49,50 . Water-column structures. The hydrographic profiling revealed temporal variation in the water-column structure at the ZHJ2 (Fig. 3b-d). The vertical salinity values ranged between 28.6 psu near the surface and 34.6 psu near the bottom. The vertical temperature structure changed from 30 °C at the surface to 21.9 °C near the bottom. The sigma-t values ranged from 18 to 24.15 kg/m 3 . Because of the salinity and temperature structures, the density field displayed strong vertical density gradients especially at around 10 m and showed the lowest water density at the surface (lowest salinity also) at end of experiment.
The above vertical structures differentiate two types of regimes. The upper water column was occupied by warm, less saline water sourced from the Pearl River effluent. Due to the on-offshore tidal oscillations, the plume regime was modulated with alternating water masses of low salinity, high temperature (during rising tide;  www.nature.com/scientificreports/ marked with the rectangles) or high salinity, low temperature (during falling tide). Furthermore, the surface water was rapidly displaced downward occasionally for short durations (e.g. selected examples marked "1" and "2" in Fig. 3b). In the lower water column, a secular trend of decreasing temperature existed and is perceived by the depth of the 22 °C isotherm from 10 to 20 m above bed. Punctuated by the tidal current, the colder water signals showed pulses (Fig. 3c).
Mixing parameter (S p ). The extent of boundary-related mixing to the water column was evaluated by the mixing parameter S p (Fig. 4a). The wind-induced S p and the current-induced S p were calculated to quantify the contributions of wind and current in the surface layer (below 3 m). When S p > 2, the water column is well stratified. When S p < 1, the water column is well mixed. The values of S p induced by currents were mostly around or greater than 2 during the experiment at ZHJ2. This implies that bottom friction drag on the current-generated mixing was insufficient to mix the entire water column 21 . However, around 79% of the wind-induced S p through air-sea boundary was smaller than 2. This suggests the wind contributed sufficiently to the mixing force in the surface layer especially during the period between 18:00 July 25th and 11:00 July 27th (marked by a red bar). It is noted that this period was the longest when the S p values were continuously lower than 2.
Stratification parameter (E). The Static Stability Index (E) quantified the degree of water-column stratification. The positive value of E signifies water-column stability. Conversely, negative E denotes unstable water column. There was a layer of high positive E values in the upper water column (Fig. 4b) and formed a dynamic barrier preventing vertical mixing as described by previous works 21,51 . However, the depth of high positive E values fluctuated vertically, and E values had occasional gaps. The gaps most frequently occurred in periods of low E values and favorable wind mixing when S p were lower than 2 marked by the red bar (Fig. 4). Additionally, comparisons of E and salinity profiles indicated the stability gaps co-existed with surface plume water drawdowns (representative examples marked "1" and "2"; Figs. 3b, 4b).
Propagation of internal solitary waves from northern South China Sea to ZHJ2. Using the echo sounder, EK60, typical NIWs are recognized along the transect from DS to the ZHJ2 (Fig. 5a-f). The echograms provide a strong evidence that the energy of the NIWs were propagating at depth approximately 25 m from the NSCS into the Guangdong coastal area. During the shoaling process, the NIWs transitioned from depression waves to elevation waves between the depth range of 75 to 100 m (also indicates the critical depth), which were 175 to 225 km away from DS (Fig. 5c,d). The results show the amplitudes of depression wave and elevation wave were around 25 m and 10 m, respectively. This deformation of shoaling NIW was also found at nearshore of Dongsha Atoll by previous studies 10, [52][53][54] .
The thermometers mounted on the mooring at DS displayed that the NIW frequently appeared at the end of July in 2016 ( Fig. 5h-i) were followed by wave packets. According to harmonic analysis results, the NIWs presented in diurnal periods that were mainly composed of K 1 and O 1 tidal components ( Table 1). The ER was less than 50% in depths shallower than 55 m possibly because of the nonlinear effect by the passing NIWs. Since the NIW is considered to be driven by the barotropic tide 4,55 , the harmonic tidal analysis was also used to predict Table 1. The tidal constituents at different depths in in both cruises (June and July) and at DS site. Only the first four tidal components (K 1 , O 1 , M 2 , and S 2 ) are listed in the table. The "x" symbol indicates the amplitude was insignificant or there was lack of data. The ER ratios indicates the tidal-to-total energy ratio, and F means form number. www.nature.com/scientificreports/ www.nature.com/scientificreports/ the arrivals of the NIW at Sta. DS until the end of experiment, and the results indicate the NIWs consistently occur diurnally (Fig. 5h).
At ZHJ2, the echogram shows that the energy of NIWs presented as the wave packets with 6-8 h duration in the upper water column (Fig. 5b,g1,g2, Supplementary Fig. 4). As wave packets arrived, the vertical water movement was perturbed, and the maximum downward displacement was 20 m. Simultaneously, the bottom sediment resuspension occurred. The NIW patterns were also recorded by CTD and current profiles. The surface warm water was drawn down as the convergent current occurred, causing the 'sandwich structure' in the F and ER above 7 m (Figs. 3b-d, 5g3-6, Table 1). This convergent current pattern was similar to the observation offshore New Jersey 12 . The low Ri values (Ri/0.25 < 1) indicated the water became unstable when the surface warm water was drawn down (Fig. 5g7-8). Afterwards, the increasing Ri value indicated the stability returned when Text "R" and "F" above the sea-surface fluctuations indicate the period from lower low-water to higher highwater (rising tide) and from higher highwater to lower low-water (falling tide), respectively. CTD profiling measurements of (b) salinity, (c) temperature, and (d) density are contoured, respectively. The white lines in (c) indicate 21.9 °C and 22 °C isotherms. The black dots above the CTD data denote times of profiling. The periods with onshore current are marked with the black rectangles. Circled numbers "1" and "2" indicate the examples of the NIWs described in the text and in following figures. www.nature.com/scientificreports/ the warm water was restored to the surface layer. However, it should be noted that since the CTD and current profiling were taken at hourly and 10-min intervals, respectively, potential aliasing could exist with capturing higher frequency NIWs. Therefore, the presence of NIWs revealed in the hydrographic and current data here should be used in a qualitative manner.

Discussion
To further decipher the process-response relations among tide, non-tidal forcing (e.g. wind), and water-column responses (salinity and temperature), a multivariate analysis tool, Empirical Orthogonal Function (EOF) analysis, was used. The tidal and non-tidal forcing were represented by the along and across-shore components in order to clarify the water-column responses to the hydrodynamic processes. The EOF analysis decomposed above spatial/ vertical and temporal co-varying variables into independent (orthogonal) modes that bear separate processresponse patterns in the data set (Fig. 6). The 1st eigenmode explains 81.7% of the correlations (standardized co-variance), and its eigenvectors show salinity co-varied inversely to temperature and non-tidal current in the alongshore component (Fig. 6a). This indicates the fresher and warmer river plume water (negative salinity and positive temperature) with the northeastward non-tidal current (positive non-tidal alongshore component) dominated the co-variability in this mode. Since the surface current direction was connected with the wind direction ( Supplementary Fig. 5), this mode describes the transport of the river plume by the northeastward wind-driven flow, which is the most important factor to influence the co-variability of the entire data set.
The contour plot of the eigenweightings of the 1st mode shows the temporal pattern of the river plume signature in the vertical structure (Fig. 6b). The line of zero crossings of the eigenweightings indicates the partition between buoyant river-plume regime at the surface and the coastal water regime underneath 21,23 . The depth of this line varied between 5 and 10 m, which delineates the lower boundary of the river plume influence. The spikes on the partition indicate water mass perturbations that is not resolved by this mode.
The influence of tides was resolved by the 2nd eigenmode (Fig. 6c,d), which explains 6.9% of the correlations. The sign of the tidal components mainly suggests the dominance of NW-SE-ward tidal flows (inversely between alongshore and across-shore tidal components). The temporal variability of vertical eigenweighting structures shows strong periodicities (Fig. 6d). During the early part of the flood tide (from lower low-water to lower highwater; periods marked with "F1" in Supplementary Fig. 6), the strong landward tidal flows dominated www.nature.com/scientificreports/ the entire water column at diurnal intervals. During the successive ebb tide (marked with "E1"), the alongshore diurnal tidal flows dominated the water column, whose strength decreased with depth. Additionally, the eigenweightings showed a phase offset in the tidal flow between the upper and lower water column at the depth around 15 m. The landward tidal component (positive eigenweighting) presents the diurnal period near the surface, but semi-diurnal period in the lower layer. Since the river plume regime accounts for 81.7% of the correlations, the influence of the tide was 'blurred' . Subsequently, the co-variability due to the river plume regime of the 1st eigenmode was removed from the original dataset 56 . EOF was again used on the residual dataset. The result enhanced the tidal influence. The diurnal and the semi-diurnal tides account for 24.1% and 43.6% of the residual correlations, respectively ( Supplementary Fig. 7). The water column shallower than 15 m was dominated by diurnal-dominated mixed tide, which was composed of K 1 tide mainly with NE-SW orientation ( Fig. 2e and Table 1). On the other hand, deeper than 15 m, the tidal current was semidiurnal-dominated mixed tide, which was comprised of M 2 tide with NW-SE current direction. The frequency analysis was also used to provide additional information of the data, and results are shown in Supplementary Fig. 8.
The 3rd eigenmode of the original dataset explains 5.1% of the correlations (Fig. 6e). The eigenvectors show the water temperature were grouped inversely to alongshore and cross-shore nontidal flows. Since alongshore  (g 3-4 ) are the synchronized temperature profiles recorded by CTD for the same periods as (g 1-2 ) which are marked by the circled numbers 1, 2 in previous figures. The black dots above the temperature contours denote the times of profiling. (g 5-6 ) are contoured cross-shore non-tidal currents, in which red color is landward, and blue color is seaward. (g 7-8 ) are contoured Richardson number. The red color indicates the stable condition, and the blue indicates the turbulence regime. (h) The temporal variation of the temperature structure in the water column measured at DS. The white circles with black in portion indicate the period from new moon to full moon. The translucent part in (h) indicates the predicted temperature variation by harmonic analysis. (i) Is the enlarged segment of (h) at the time when the last two NIWs passed before retrieving the mooring at the DS. The time marked with black dashed rectangle corresponds to the black dashed rectangle marked on the spatial echogram (f). www.nature.com/scientificreports/ non-tidal current is regarded as wind-driven (Fig. 2, Supplementary Fig. 5), this mode demonstrates the coldwater intrusion related to the wind, which leads to upwelling near the coast 37,42,43 . The contour plot of eigenweightings indicates the thickening of the cold water in the lower column that gradually extended upward into shallower depths (the red patch in Fig. 6f). Ultimately, the intrusion reached 20 m above bed. The eigenweightings also depict that the water temperature fluctuated at semi-diurnal frequencies between 06:00 July 25 and 12:00 July 27, which corresponded to the across-shore tidal currents (Figs. 3c,6f). The EOF analysis clearly established a hierarchy in the process-response pattern at a fixed point on the pathway of the Pearl River plume. Firstly, as the fresher and warmer Pearl River effluent exited the PRE, the typical southwesterly summer monsoon and northeastward wind-induced current were the primary forcing to cause the river plume to disperse northeastward and seaward ( Supplementary Figs. 1 and 2). Then, the Pearl River effluent is propagating along the Guangdong coast and enters the Taiwan Strait 36,42,57 . On this setting, the presence of the buoyant river plume caused strong stratification in the upper part of the water column along the river plume dispersal pathway (between 5 and 10 m depth defined by EOF; Figs. 3b-d, 4b, 6a,b). The stratification created a density barrier against water exchange between river plume water and ambient seawater beneath (Fig. 4). Above the stratified layer, the mixing was induced by wind because the low S p indicates the prevailing wind-mixed environment. The wind induced mixing could disrupt the in-situ stratification occasionally when the water density gradient was low 21,27 (marked as the red bar in Fig. 4b). www.nature.com/scientificreports/ The disruption in the water stratification was also caused by the vertical water displacement that entrained the buoyant plume water downward abruptly that appeared as the distinctive pulses (marked "1" and "2" in Figs. 3b, 4b, 5g1-8). Eventually, turbulence was generated (Ri < 0.25) and caused temporal gaps in the stratification. Our observations imply that the events of water displacements corresponded to the passing NIWs, which were propagating from the northern SCS basin westward to offshore PRE (Fig. 5). This interpretation is also supported by the satellite images that indicate internal wave fronts propagate shoreward on the shelf in the northern SCS 58,59 ( Supplementary Figs. 9, 10).
The NIW-induced water displacements occurred at diurnal frequencies offshore PRE (Figs. 3b-d, 4b, 5g), because NIWs from the SCS remained in the tidal signals along the propagating pathway 53 . After the NIWinduced mixing energy dissipated, the water-column stability returned. Since the stratification prevents vertical exchange between the river plume water and ambient seawater, the NIW-induced perturbations could contribute to mixing processes among different water masses offshore PRE. To the best of our knowledge, our study is the first occasion to observe the influence of far-field NIWs on the river plume dynamics and local mixing processes. Although other studies pointed out that mechanisms (e.g. the plume-related hydraulic jump and tidal pumping etc.) might trigger the similar phenomenon 31,60-62 , our measurements (e.g. strong across-shore current variation, consistently NE wind; Figs. 2a and 5g5-6) demonstrates that the perturbating energy was from the offshore.
Tidal motions, which were the second most important forcing, resulted in the diurnal and semi-diurnal fluctuations in the signature of the buoyant river plume (Fig. 3). Flooding tide pushed the buoyant plume landward, and ebbing tide pushed the plume water seaward. Thus, noticeable temporal variability in the vertical structures of salinity and temperature occurred, and water stratification fluctuated diurnally. Similar mechanisms have been discussed in river plumes of the Changjiang River and Columbia Rivers 41,63,64 .
Although tidal processes only account for small percentage of the total co-variability in the dataset, it is by no means that tidal forcing was less complex in the hydrodynamic mechanisms. After removing the 1 st eigenmode component from the data, the tidal processes explained 67.7% of the residual co-variability. The results indicated the vertical partition of the tidal influence was around 15-m depth. Above this depth, the diurnal tide, K 1 , dominated the tidal forcing; below this depth, the semi-diurnal tide, M 2 , controlled the tidal forcing (Figs. 2e, 6c,d, and Supplementary Figs. 7, 8).
The tertiary mechanism was the cold-water intrusion in the lower water column related to the wind-driven coastal upwelling (Figs. 3c, 6e,f). Throughout the study period, the vectorially averaged southwesterly winds were upwelling favorable 37 (Fig. 2a). Therefore, the offshore water flowed northwestward near the bottom (Fig. 2b), leading to the intrusion of cold water near bottom (Fig. 3c). The thickness of the bottom water mass with vertically homogeneous temperature structure, defined by the depth of the 22 °C isotherm, which gradually moved upward in the water column (Figs. 3c, 6f). However, when the wind direction changed to non-upwelling favorable conditions due to the influence of the Tropical Storm Mirinae, the bottom current started to stagnate, and the intrusion of the offshore water paused (Figs. 2a,c, 3c). Within this period, diurnal and semi-diurnal tidal currents became the dominant factor to modulate the bottom water-mass movements (Figs. 3c, 6f). The intrusion resumed with the increasing magnitude of northeastward current as the southwesterly wind started blowing again. This implies that the contributions of the wind and tidal forcing to changing the bottom water mass property cannot be neglected and should be considered when evaluating biogeochemical mechanisms for coastal hypoxia at the nearby Mirs Bay 29,65 . Although previous studies have reported that the prevailing southwesterly wind induced a net landward bottom flow by using modelling 37,43,47 , our research provides the temporal variation of bottom current filed and its response to the wind with in-situ measurement.
In addition to the process-response hierarchy, the EOF analysis also revealed vertical partition between adjacent regimes at discrete depths. The boundary between the wind-driven river plume and ambient shelf regimes existed between 5 and 10 m depth. Within the tidal regime, a partition appeared at about 15-m depth in both the semi-diurnal and diurnal frequencies. At this point, the mechanism for this partition is subject to future studies. The last partition is between the cold-water intrusion regime near the seafloor and the overlying ambient shelf regime. During our study, the boundary of this regime moved toward the surface, suggesting increased flux of cold-water mass.
Major findings of this study are illustrated in a schematic diagram in Fig. 7. The northeastward currents driven by upwelling favorable monsoon winds dispersed the Pearl River plume in the NE direction and triggered the intrusion of cold-water mass near the bottom (Fig. 7a,b). The buoyant plume water created strong stratification at the surface, which set the 'stage' for other physical processes. Since tidal currents caused the plume water mass to swing landward in the flood and seaward in the ebb, it influenced the thickness of the plume water mass and stratified layer (Fig. 7b). The wind field also could increase the thickness of the stratified layer and disrupt the water stability occasionally when the density gradient was low in the water column 21,27 (red bar marked in Figs. 4b, 7b). Furthermore, in some instances, the stratification in the river plume regime was disrupted by the shoreward propagating NIWs from the northern SCS (Fig. 7c). The NIWs created perturbations to draw surface water downward (Fig. 7b). Our findings provide a new insight into how the NIWs influence the mixing processes in the far field of the river plume dispersal. This could improve our understanding of the transport and dispersal processes of the terrestrial material sourced from the PRE to the oceanic sink in NSCS shelf and Taiwan Strait.

Methods
Shipboard observations. Two cruises were carried out by using R/V Ocean Researcher III on June 3rd-5th at ZHJ1 and R/V Ocean Researcher I on July 24th-29th at ZHJ2 for shipboard monitoring and mooring deployment and retrieval (Fig. 1a) www.nature.com/scientificreports/ At ZHJ1, the 75 kHz shipboard ADCP was used to obtain the current field with 1-min sampling interval and 0.12 m/s precision. The first bin of shipboard ADCP was at 16.5 m with 8-m bin size. At ZHJ2, hourly hydrographic profiling was conducted by using Seabird SBE 19 plus CTD with 0.2-m bin size. Along a transect from DS to ZHJ2 (insert in Fig. 1a), the acoustic backscatter strength profile was recorded by SIMRAD EK60 echo sounder (38 kHz) with 1-s sampling interval and 0.4-m bin size to capture the propagating NIWs in the northern SCS. Meteorological data including atmospheric pressure, wind speed and direction were recorded at 10-min intervals by the weather station onboard the R/Vs.

Mooring deployments.
A 44.7-m taut-line mooring was deployed nearby ZHJ2 (Fig. 1a,b) on July 24th-29th. The mooring was configured with Teledyne RDI 1200 kHz ADCP (upward-looking at 11 m below the surface) and Teledyne RDI 300 kHz ADCP (downward-looking at 12.7 m). The bin-sizes for two ADCPs were 0.5-m (1200 kHz) and 1-m (300 kHz) with blank distance 0.5-m and 1-m, respectively. The velocity resolution both for 1200 kHz and 300 kHz ADCP was 0.1 cm/s. The precision of the pressure sensor mounted on the ADCP was 0.01 decibar. The sampling intervals for all instruments were synchronized at 10-min intervals.
Another 150-m long taut-line mooring was deployed at Station DS (117.248°E, 21.8432°N; Fig. 1a,c), north of Dongsha Atoll, near the Sta. S7 conducted by Ramp et al. 53 . The mooring was mounted with 4 DST milli-TDs (marked T/P; thermometer with the pressure sensor) and 8 8-bit mini-logs (marked T; thermometer) at 27.8, 46.7, 57.6, 68.8, 91.9, 103.9, 116.2, 128.7, 141.6, and 150.7 m depths, respectively. The sampling rate of each sensor was 6-min to capture the signature of the passing NIWs between June 3rd and July 23rd. The accuracy of pressure measured by DST milli-TD is ± 1.5 m, and the accuracy of temperature recorded by the 8-bit mini-log is ± 0.3 °C.

Hydrodynamic analysis and mixing indicators.
The instantaneous current field is divided into two components, non-tidal current and tidal current, in this study. The tidal components in the current data were  where V p is the predicted tidal current including barotropic and baroclinic tide, m is the number of tidal components used in the analysis; α k is the amplitude of the kth tidal component; ω k and θ k indicate the tidal angular frequency and corresponding phase lag, respectively; ω k = 2π/T k , where T k is the period of the kth tidal component; t is the sampling time. In this study, only tidal components of which the signal to noise ratio (SNR) higher than 2 were used. The detailed information about calculation of SNR is described by Pawlowicz et al. 66 .
The tidal-to-total energy ratio (ER) in the current was calculated to quantify the energy contribution of the synthetic tidal flow to the instantaneous flow 67 and defined as: where n is the length of the data series. V p is the mean tidal current velocity.V is the instantaneous current velocity. When the ER value is lower, it indicates the contribution to fluctuating the current field by the tidal forcing is weaker.
The form number, F, is used to identify the type of the tidal environment. It is defined as the amplitude ratio between major diurnal tidal and semi-diurnal tidal components. F = 0 ~ 0.25: Semi-diurnal tide. F = 0.25 ~ 1.5: Mixed tide with semi-diurnal dominance. F = 1.5 ~ 3: Mixed tide with diurnal tide dominance. F > 3: Diurnal tide. The static stability index (E) was calculated to quantify the water column stability 21,23,51,68 . The water column is stable when E is positive and unstable when E is negative.
where ρ 0 is average water density in the measurement, ∂ρ w and ∂ z is the vertical gradient of the water density and the water depth. The z-axis is positive upward. The buoyancy frequency is also calculated and shown in Supplementary Fig. 11 as the reference.
The stratification parameter, S p , was used to reveal the relative contribution of tidal mixing and wind mixing in the water column 69 . S p is expressed as: when S p > 2, is well stratified, and when S p < 1, the water column is well-mixed. ε stands for the energy dissipation rate per unit mass, and ε induced by wind and current can be expressed as: where ρ a is the air density (kg/m 3 ), C 10 (~ 0.0016) is the drag coefficient for the wind at 10 m height, K is the wind factor (~ 0.004), U 10 is the wind speed at 10 m height, ρ 0 is the average water density in the measurement, and z is the water depth at specific layer, and 3 m was used in this study.C D is the bottom drag coefficient (~ 0.0025), and Ū p is the depth-average tidal velocity (m/s), h is the water depth from bottom to the specific layer (41.7 m used in this study).
The mixing mechanism induced by the passing NIWs was assessed by the Richardson number, Ri. Ri value greater than 0.25 indicates the stable environment; otherwise unstable as the turbulence overcomes stratification. The Ri shown in Fig. 5g7, 8 was divided by 0.25 to normalize the value.
where N is the buoyancy frequency, N is expressed as N = √ gE . g stands for the gravity acceleration. S 2 is the current-induced vertical shear velocity in the horizontal dimension and is expressed as (4) E = -1 ρ 0 ∂ρ w ∂z (5) S p = log 10 1 ε (6) ε wind = 10 4 × ρ a C 10 K|U 10 | 3 ρ 0 z www.nature.com/scientificreports/ and Information System (EOSDIS). We appreciate the help from the staff and students of the Coastal Geology and Processes Laboratory and Yu-Shih Lin's Laboratory in the Department of Oceanography, National Sun Yatsen University. We thank two anonymous reviewers for their helpful comments and suggestions to improve the manuscript.