To investigate the existence of internal tide turning depths in the SCS, the hydrographic data collected by the open cruise of the NSFC in the SCS in March of 2012 is analyzed. During the cruise, a total of 48 sets of hydrographic profiles were obtained, and 35 of these profiles have observational depths exceeding 3 000 m (see locations in Fig. 1). The observational depths of all the profiles exceed 200 m. The extensive horizontal spatial coverage and vertical depth range of the conductivity-temperature-depth (CTD) profiles allow an accurate analysis of the existence of internal tide turning depths in the SCS. These CTD profiles were recorded by the Seabird SBE-911 and have a temperature accuracy of 1 mK and a conductivity accuracy of 0.3 mS per meter. The vertical resolution of the data is 1×10⁴ Pa.

To supplement the NSFC hydrographic measurements and uncover more observational evidence of internal tide turning depths in the SCS, WOD18 data, which represent the most comprehensive collection of ocean hydrographic data currently available, were also analyzed. The WOD18 data include quality-controlled observational data for more than 15.7 million temperature profiles and 8.5 million salinity profiles (Boyer et al., 2018) and are available from the National Oceanographic Data Center (https://www.nodc.noaa.gov/OC5/WOD/pr_wod.html). Given the main focus of this study, only WOD18 high-resolution CTD hydrographic profiles within 4°–24°N and 105°–121°E with the deepest measurement exceeding 200 m are considered. After screening, a total of 2 077 WOD18 profiles from the years 1773 to 2019 met these requirements, within which 23 profiles have observational depths greater than 3 000 m (15 deep-sea measurements in the SCS and 8 deep-sea measurements in the Sulu Sea). All the analyzed WOD18 profiles are at the observation depths.

The WOD18 stations cover nearly the entire SCS (Fig. 1); they have better spatial coverage than the NSFC observations due to the longer period of data accumulation, but the NSFC observations have more deep-sea hydrographic measurements. To obtain more abyssal data, gridded data products such as the US Navy Generalized Digital Environment Model (GDEM) is also analyzed, which are available at 78 standard depths (ranging from the surface down to 6 600 m; Carnes, 2009). However, temperature-salinity plots indicate that the abyssal GDEM data show obvious data interpolation (not shown here). Calculating the buoyancy frequencies in weakly stable waters utilizing original profiles is more reliable than using a gridded dataset. The final dataset used in this study consisted of 48 NSFC and 2 077 WOD18 temperature/salinity profiles, with 58 profiles having observational depths exceeding 3 000 m. Using this dataset, this study aimed to obtain a spatial map of internal tide turning depths in the SCS.

Quality control was applied to the raw temperature and salinity profiles from each station; duplicate profiles and abnormal values were removed. The original hydrographic data have a “wiggly” nature due to the observation accuracy of the instruments (Fig. 2). For instance, in the WOD18 data, the accuracies of the temperature and conductivity sensors are 1–5 mK and 0.3–2 mS per meter, respectively (Boyer et al., 2018). Fluctuations in temperature and salinity data can lead to large fluctuations in the buoyancy frequency. Therefore, the raw hydrographic temperature and salinity data ( ${T}_{\rm{raw}}$ and ${S}_{\rm{raw}}$, respectively) were smoothed before calculating the buoyancy frequency. The entire water column was vertically divided into several bins with a constant depth, and ${T}_{\rm{raw}}$ and ${S}_{\rm{raw}}$ in each bin were averaged (Fig. 2).

The ocean stability is weak near the surface and at great depths; thus, the determination of turning depths requires the most precise buoyancy frequency calculation possible. Van Haren and Millot (2006) sought to determine the buoyancy frequency in weakly stable waters and concluded that CTD observations can be used to estimate water stability with an error of ~0.8 f (with an average over 1×10⁶ Pa, where f is the inertial frequency). In the SCS (~10°–20°N), 0.8 f ranges from 2.53×10^–5 rad/s to 4.98×10^–5 rad/s, which is less than the tidal frequency (1.41×10^–4 rad/s and 7.29×10^–5 rad/s for semidiurnal and diurnal frequencies, respectively); thus, the observed temperature and salinity profiles of the CTD data can indicate the existence of internal tide turning depths in the SCS.

The buoyancy frequency indicates the natural frequency of vertically oscillating fluid parcels in the continuously stratified water and can be derived in terms of the vertical gradients of potential temperature $ T $, and absolute salinity $ {S}_{{\rm{a}}} $, such that

where $ g $ is the gravitational acceleration, and $ \alpha $ and $ \beta $ are the depth-dependent thermal expansion coefficient and the saline contraction coefficient, respectively. $ {N}^{2}\left(z\right) $ is not calculated by the vertical potential density gradients because the density of water is not a directly observed variable, and calculating the buoyancy frequency utilizing the potential density expression may underestimate the strength of the stability (Lynn and Reid, 1968). In this study, the Gibbs SeaWater toolbox was adopted to calculate the buoyancy frequency (http://www.teos-10.org). This toolbox is based on the new thermodynamic seawater equation TEOS-10, and the ocean state is depicted more accurately by considering the spatially dependent $ {S}_{{\rm{a}}} $ (Li et al., 2013).

As suggested by King et al. (2012), the squared buoyancy frequency ( $ {N}^{2}\left(z\right) $), rather than the buoyancy frequency was calculated to determine the existence of a turning depth because the squared version was used in the internal wave governing equation. The squared buoyancy frequency can also be negative, indicating convective instability. Vertical profiles of the bin-averaged squared buoyancy frequency ( ${N}_{\rm{bin}}^{2}\left(z\right)$) and the associated standard deviation ( ${SD}_{N{^2_{{\rm{bin}}}}}$) were obtained by utilizing a Monte Carlo method following the calculation process of King et al. (2012). First, the bin-averaged temperature and salinity ( ${T}_{\rm{bin}}$ and ${S}_{\rm{bin}}$) were spline-fitted to the entire water column, and the standard deviations of temperature and salinity ( ${SD}_{T_{\rm{bin}}}$ and ${SD}_{S_{\rm{bin}}}$) and the number of samples in each bin ( ${M}_{\rm{bin}}$) were obtained. Then, a set of synthetic data samples was devised based on sample functions (2) and (3); the number of samples was set to 500 following King et al. (2012).

Finally, $ {N}^{2}\left(z\right) $ was calculated based on the synthetic data, and ${N}_{\rm{bin}}^{2}\left(z\right)$ and ${SD}_{{N^2_{\rm{bin}}}}$ were obtained. Generally, $ {N}^{2}\left(z\right) $ is large at the main pycnocline (a gradual interface separating two fluids of different densities) and decreases exponentially from the pycnocline to both the surface and bottom boundaries. In the SCS, the main pycnocline depth is approximately 120 m, which was obtained from satellite altimetry estimation (Chen et al., 2017). Thus, the distance from the main pycnocline to the ITLTD is far more than that to the ITUTD, and the bin widths are different for ITLTD and ITUTD identification. The sensitivity of the calculated results to different bin widths is shown in Section 3.1.

The buoyancy frequency commonly peaks at the depth of the pycnocline and decreases from the pycnocline to both the surface and bottom boundaries. The existence of ITLTDs was first investigated. An example of the deep-sea $ {N}^{2}\left(z\right) $ calculated at one station (14.1°N, 117.0°E) is shown in Fig. 3; obviously, the depth-dependent $ {N}^{2}\left(z\right) $ crosses the M₂ and K₁ tidal frequencies at depths of approximately 3 660 m and 3 870 m, respectively. The K₁ ITLTDs are greater than those of M₂ at the same locations due to the lower tidal frequency. The distance between the ITLTDs and the seafloor is approximately 380–590 m, which corresponds to a thick, deep, “homogeneous” layer of weak stability. The calculated ITLTDs are basically consistent under different bin values of 50 m, 100 m and 200 m (Fig. 3). These findings confirm that ITLTDs exist in deep marginal seas in addition to the abyssal open ocean.

Notably, an optimum bin size must be selected to calculate the turning depth. The squared buoyancy frequency profiles may swing drastically due to small bin values (Fig. 3a), which could lead to an incorrect supposition of the existence of several turning depths. Buoyancy frequency profiles become smoother with a larger bin size and are more representative of the vertical stratification structure. However, at the other extreme, buoyancy frequency profiles may not be able to cross tidal frequencies due to large bin sizes because large bin values correspond to a reduction in the number of deep samplings (Fig. 3c). Thus, in a subsequent analysis, a 100 m bin size as suggested by King et al. (2012) was selected to avoid multiple crossings between the squared buoyancy frequency profile and the tidal frequencies and to simultaneously ensure the reliability of the calculated ITLTD. However, the 100 m bin size could also result in multiple crossings, and the ITLTDs at these stations were obtained by calculating the mean of the crossing depths.

In contrast to the very limited abyssal observations, a stronger data foundation was available for the investigation of the ITUTDs. Figure 4a shows a sample of upper ocean $ {N}^{2}\left(z\right) $ values calculated at one station (20.25°N, 117.21°E), in which the thickness of the upper mixed layer (within which the density discrepancies are less than 0.01 kg/m³) is approximately 25 m. A closer view of the vertical distribution of $ {N}^{2}\left(z\right) $ values within the mixed layer (Fig. 4b) shows that the water stability in the upper mixed layer is very weak, and the calculated $ {N}^{2}\left(z\right) $ profiles with small bin values cross the tidal frequencies many times. Similar to the ITLTD calculation, the bin values were increased to smooth the $ {N}^{2}\left(z\right) $ profile. Considering that the pycnocline is much closer to the ocean surface than the ocean bottom, $ {N}^{2}\left(z\right) $ decreases upward from the maximum value at the pycnocline to tidal frequencies within a very short distance if ITUTDs exist; the bin values in the ITUTD calculation are much less than those in the ITLTD calculation. When the bin sizes are greater than 5 m, the $ {N}^{2}\left(z\right) $ profiles are smooth and are still able to cross tidal frequencies, and this result confirms the existence of ITUTDs and indicates the fact that propagating internal tides can reflect at the subsurface. However, the estimated ITUTDs are different under different bin values and it is difficult to determine an optimal bin value convincingly. One possible reason is that, in contrast to the relatively quiet abyssal ocean, the upper boundary is influenced by winds and surface buoyancy fluxes, and the $ {N}^{2}\left(z\right) $ profiles may not exhibit a monotonous decreasing trend from the pycnocline to the sea surface. The coexistence of multiple ITUTDs at one single profile may be closer to the real oceanic scenario. Thus, ITUTD estimations include more uncertainty than ITLTD estimations. In the following section, a 5 m bin size is selected to obtain a reference value for ITUTDs in the SCS. If there are multiple crossings between the $ {N}^{2}\left(z\right) $ profiles and the tidal frequencies, ITUTDs are determined by calculating the mean of the crossing depths.

Based on the $ {N}^{2}\left(z\right) $ profiles calculated from both the NSFC and WOD18 observations, a comprehensive picture of turning depths in the SCS can be obtained. Among the 58 abyssal ocean observations, ITLTDs exist at 16 (and 12) stations for M₂ (and K₁) internal tides (Figs 5 and 6). At certain stations, ITLTDs may also exist but are not detected because of the incomplete observational coverage of the full water depth. Because K₁ internal tides have lower frequencies, in rare cases, the minimum $ {N}^{2}\left(z\right) $ may cross only M₂ tidal frequencies but fail to reach K₁ tidal frequencies. Horizontally, ITLTDs are distributed over a wide spatial area in the SCS (Figs 5a and 6a) but are concentrated mainly in the southern part of the deep basin. This north–south asymmetry may result from the fact that the water is deeper in the south. The maximum ITLTD exceeds 4 500 m in the vicinity of the Manila Trench. ITLTDs also exist in the Sulu Sea according to limited observations. Vertically, the distances between ITLTDs and the seafloor range from 270 m to more than 1 200 m (Figs 5b and 6b), which implies the existence of a broad evanescent region at multiple locations at the bottom. The maximal distance between the turning points and the seafloor occurs near the southern boundary of the deep basin. The residence time of the deep water in the SCS deep basin is approximately 30–100 years (Broecker et al., 1986; Qu et al., 2006), and $ {N}^{2}\left(z\right) $ profiles for depths greater than 1 000 m can exhibit very limited variations over several years (King et al., 2012). Therefore, the map presented in this study may shed light on the basic spatial features of ITLTDs in the SCS. More high-quality deep-sea hydrographic measurements are needed in the future to complete the details on the map.

Figure 7a shows the spatial distribution of M₂ ITUTDs in the SCS. At each station, the calculated ITUTDs are nearly the same for M₂ and K₁ frequencies (Fig. 4c); therefore, the map of K₁ ITUTDs is not shown here. Among the 2 125 ocean observations, upper turning depths exist at 210 stations for both M₂ and K₁ internal tides. Most of the stations with the presence of ITUTDs are located in the northeastern SCS, where the most intensive observations are located. The maximum calculated ITUTD is 93.7 m; this indicates that the internal wave evanescent region near the surface is not as spacious as that in the abyssal bottom. Moreover, the calculated ITUTDs are different among neighboring stations, which may be caused by the difference in observation times. In contrast to the relatively quiet abyssal ocean, the upper ocean is stirred by winds and surface buoyancy fluxes. Thus, ITUTDs can be highly variable both temporally and spatially. To quantify the uncertainties of the calculated ITUTDs, the standard deviations of the obtained ITUTDs under different bin values of 2–8 m were calculated (Fig. 7b). Note that the standard deviations were calculated only at stations where ITUTDs could be detected under all of the abovementioned bin values. The standard deviations are generally less than 10 m at most stations; at a few stations, the standard deviations can exceed 15 m. Compared to the ITUTDs in the 5 m bin, the percentage of uncertainty ranges widely from 0.7% to 80.7%. Further studies with the support of long-term mooring data may be more convincing in understanding the uncertainties and temporal variations in the ITUTDs.

Considering their relatively shallow water depths, oceanic motions generally have very large horizontal scales, rendering the cosine terms of the Coriolis force insignificant such that they are typically neglected in traditional approximations. In the previous sections, turning depths were estimated based on the traditional Coriolis definition. However, when approaching a turning depth, the internal tidal beams are nearly perpendicular to the horizontal plane (Paoletti and Swinney, 2012); thus, the cosine term of the Coriolis force can no longer be neglected for motions with manifestly vertical characteristics (Gerkema et al., 2008). As the cosine terms of the Coriolis force act only in the zonal direction, the horizontal propagation direction of internal tides should be taken into account. For internal tides propagating at an angle $ {d} $ with respect to the east, the critical $ {N}^{2}\left(z\right) $ can be derived following King et al. (2012):

where $ {f_s} = 2\Omega \sin \varphi $ (in which $ \Omega $ is the Earth’s angular velocity and $ \varphi $ is the latitude in the northern hemisphere) and $ {f_c} = 2\Omega \cos \varphi $ are the sine and cosine terms of the Coriolis parameter, respectively. For zonally propagating internal tides, $ d $ equals zero (or 180°), and $N_{\rm{crit}}^2 = {\omega ^2}$, which is the same as the traditional approximation.

The latitudinal dependence of the critical value of $N_{\rm{crit}}^2/{\omega ^2}$, considering the full Coriolis force and a given internal tide propagation direction, is shown in Fig. 8. For a fixed propagation direction, the discrepancy of the critical $ {N}^{2}\left(z\right) $ value between the traditional approximation and the full definition increases poleward. The biggest difference occurs in the vicinity of the critical latitude (where the inertial frequency equals the tidal frequency, i.e., 74.5° for the M₂ constituent and 30° for the K₁ constituent). The SCS deep basin extends at relatively low latitudes from approximately 10°N to 20°N (the light gray areas in Figs 8a and b), but $N_{\rm{crit}}^2$ can still be significantly affected by the propagation direction. When internal tides propagate at a small angle relative to the zonal direction (less than the critical latitude), the discrepancy between the traditional approximation and the full definition increases with the angle (Figs 8a and b). By contrast, when the angle between the propagation direction and the zonal direction is greater than the critical latitude (the dark gray areas in Figs 8c and d), the $N_{\rm{crit}}^2$ term in Eq. (4) is always negative when it is equatorward of the critical latitude, and internal tides can propagate freely without any turning depth limitation. M₂ internal tides radiate from the Luzon Strait to the SCS along a stronger north-westward beam and a weaker south-westward beam (Zhao, 2014); the southwestward beam directs to the deep basin with an approximately fixed propagation direction of approximately 225° (a 45° angle relative to the zonal direction). The critical value of $N_{\rm{crit}}^2$ for M₂ internal tides under the full definition is approximately half of that under the traditional approximation, indicating a greater turning depth. K₁ internal tides radiate from the Luzon Strait to the SCS deep basin along one main beam with a propagation direction varying between approximately 190° and 240° (a 10°–60° angle relative to the zonal direction) under the effect of the Earth’s rotation (Zhao, 2014). The situation becomes slightly more complicated for the K₁ internal tide. During the first half of propagation (with a propagation direction from 190° to 210°), the discrepancy of the critical value of $N_{\rm{crit}}^2$ between the traditional approximation and the full definition increases significantly. Moreover, during the second half of propagation (with a propagation direction from 210° to 240°), K₁ internal tides are independent of the turning depths. These results are remarkably different from those under the traditional approximation. In addition, the critical latitudes can be shifted by several degrees under the impact of the background relative vorticity (Yang et al., 2018), which suggests that an accurate prediction of the propagation of internal tides in the SCS is still very challenging.