The historical hydrographic data of full-depth temperature and salinity profiles collected by the SCS open cruises and several research field campaigns during the past 15 years (2004-2018) are used in this study. These data were obtained from a SBE-911plus conductivity-temperature-depth (CTD) system using frequencies between 8–24 Hz. After pre-and post-cruise calibrations, the accuracies of the CTD sensors are 0.0003 S m^-1 for salinity and 0.0018°C for temperature. To get as complete a vertical profile of the BML as possible without damaging the CTD sensors by colliding with the bottom, acoustic altimeters were used in the more recent cruises to monitor the distance of the sensors to the bottom. In this study, only the downcast data that with maximum observed depth less than 50 m from the bottom are used because those data may more easily observe the structure of the BML.

The raw CTD data quality control applied the following criteria: (i) remove CTD profiles where the original information about the station position and/or water depth is missing or incorrect; (ii) because only deeper locations are considered in this study, stations with depths less than 100 m are excluded; (iii) down sample the raw vertical high-resolution to 1 m and apply a 5 m running mean filter to smooth the data. Applying these criteria and after validation, only 514 CTD profiles remained. Nevertheless, these selected CTD data records covered nearly the entire SCS north of 13°N and provided more than adequate coverage of the northern continental slope region of the SCS ( Figure 1A ). While the majority of the selected CTD profiles were located in depths shallower than 500 m, there were 120 CTD profiles at locations deeper than 1000 m ( Figure 1C ). Most of those profiles were collected from August and September ( Figure 1D ) because the routine SCS open cruises were conducted in those two months.

Two bottom-anchored moorings were particularly deployed at two sites ( Figure 1A ) for obtaining time series to characterize the temporal variations of the BML. One was deployed on the northern continental slope of the SCS (M1) and the other was in the deep basin in the western SCS (M2), in the latter of which a previous study had confirmed the existence of a strong deep western boundary current (Zhou et al., 2020). To obtain simultaneous observations, these two moorings were both deployed in August 2017 and recovered in September 2018, collecting a 14-month long time series for use in this study. The moorings had seven RBR-TDs and seven SBE 37 CTDs measuring temperature and pressure near the bottom at M1 and M2, respectively. The accuracies of the RBR-TDs were ±0.002°C for temperature and ±0.05% over the full-scale range for pressure. The accuracies of the SBE 37 CTDs were ±0.002°C for temperature and ±0.1% over the full-scale range for pressure. The design and configuration of the moorings used in this study are shown in Table 1 . Since there were no salinity measurements on the M1 mooring, the present study only examines the variation of temperature in the near-bottom regions at the two sites. To explore the low frequency variations of temperature, a 72-h low-pass filter was used to remove the inertial, tidal, and other high-frequency signals (Thomson and Emery, 2014). All data were averaged over an hourly interval.

To help analyze the potential influence of tidal currents on the distribution of the H_BML, the barotropic tidal currents computed from harmonic constituents provided by the latest TPXO9‐atlas tidal models (Egbert and Erofeeva, 2002) are used to calculate the tidally driven mixing and dissipations of SCS. The TPXO9‐atlas is a 1/30° resolution global model of ocean tides, which represents optimal least squares fit of the Laplace tidal equation to satellite altimetry data. TPXO9 atlas provides the eight major tidal (M2, S2, N2, K2, K1, O1, P1, Q1), three long period (Mf, Mm, 2N2) and three non-linear (M4, MS4, MN4) harmonic constituents. In this study, only the eight most energetic tidal harmonic constituents of the TPXO9‐atlas solution were used to calculate the yearlong (2017) hourly time series of barotropic tidal currents in the SCS.

According to the vertical mixing parameterization scheme as proposed by St. Laurent et al. (2002), which has been applied to estimate the diapycnal mixing induced by internal tides (St. Laurent et al., 2002; Wang et al., 2016; Tan et al., 2022), the turbulent dissipation rate ϵ and diapycnal diffusivity k_v can be calculated as follows:

where Γ is the mixing efficiency taken to be 0.2 (Osborn, 1980), q = 0.3 is the local tidal dissipation efficiency as suggested by St. Laurent et al. (2002), ρ is density of seawater, N² is the squared buoyancy frequency, and k₀ is the background diffusivity (1×10^-5 m² s^-1). F(z) is the function for the vertical structure of the dissipation, chosen tosatisfy energy conservation within an integrated vertical column, ∫ 0 − H F ( z ) d z = 1 . Because only the barotropic tidal flow is considered in this study, thus the term F(z) is taken to be 1. E (x, y) is the energy flux per unit area transferred from barotropic to baroclinic tides, formulated as

where ρ₀ is the reference density, N_b is the buoyancy frequency at the seafloor, k and h are the wavenumber and amplitude scales for the topographic roughness, respectively. The wavenumber is set to k = 2π/(10 km), we take the horizontal scales of O (10 km) as typical of the roughness. The h² is defined as the variance of bathymetry over a 1/4°×1/4° domain square (using GEBCO_2022 gridded bathymetric dataset). u_bt is the averaged horizontal speed of the barotropic tides over a yearlong time series. In this study, the density of seawater and buoyancy frequency are calculated from the GDEMv3 database (Carnes, 2009).

In addition, a two-dimensional map of the internal tidal dissipation dataset was also used in this study. The dataset consists of global column-integrated maps of internal tide energy sources and sinks with a horizontal resolution of 0.5° × 0.5°. In this dataset, energy sinks are provided for each of M2, S2 and K1 and for “All constituents” (the eight most energetic tidal constituents). The energy sinks are decomposed into five process contributions: (i) dissipation of low modes via wave-wave interactions; (ii) dissipation of low modes scattering by abyssal hills; (iii) dissipation of low modes critical reflection; (iv) dissipation of low modes shoaling; (v) local dissipation of high modes. Units are Watts per square meter. Considering the dissipation of low modes via wave-wave interaction mainly occurs in the stratified water column away from boundary layers (Polzin and McDougall, 2022), thus only the other four processes were used in this study. Detailed information about the dataset and documentation can be found in de Lavergne et al. (2019).

In this study, a relative variance method is used to identify the H_BML in the SCS. This method is based on the ratio between the standard deviation and the maximum variation of the temperature, salinity, or density profiles above the sea bed; the position of the minimum relative variance is defined as the top of BML (Huang et al., 2018b). The relative variance method is an objective method that determines the H_BML that is less dependent on arbitrary criteria (Huang et al., 2018b). Although the relative variance method was first proposed to identify the surface mixed layer, its performance in determining the H_BML is also superior to other available methods (Huang et al., 2018a). A detailed description of the method and implementation can be found in Huang et al. (2018a) and Huang et al. (2018b).

We use the relative variance method separately on the profiles of potential temperature, salinity, and potential density to obtain three estimated values of the H_BML in each CTD cast. The quality index (QI) defined in Lorbacher et al. (2006) was used to evaluate the quality of the estimate of H_BML, and profiles with QI<0.5 were discarded. It should be noted that there may exist real differences between the three H_BML estimated values, and the value with higher QI is adopted as the observed H_BML. In addition to the formal objective analysis of the H_BML, all profiles used in this study were visually inspected to detect possible errors due to contaminated samples or accidental spikes. Figure 2 shows an example of temperature and salinity profiles collected at 117.06°E, 21.45°N near the Dongsha Islands, where a well-mixed layer clearly exists in the near-bottom zone with an H_BML of about 100 m.

We estimate the vertical eddy diffusion coefficient by using the advection-diffusion equation and mooring observations. Assuming that vertical advection is balanced mainly by vertical diffusion (Munk, 1966), the momentum equation of temperature excluding the source and sink terms, becomes

where T is the potential temperature, t is time, z is the vertical coordinate (positive upward), w and A_z are the vertical velocity and vertical eddy diffusion coefficient, respectively. In Eq. (4), the three terms ∂ T / ∂ t , ∂ T / ∂ z , and ∂ 2 T / ∂ z 2 can be estimated from the mooring observations, thus the unknown values of w and A_z can be estimated from a set of linear equations using the least-square fitting method (Thomson and Emery, 2014).

To investigate the distribution of H_BML in the SCS and its sensitivity to ocean topography, two main aspects representing the topographic effects are considered. One is the topographic slope angle (θ) and the other is the topographic ruggedness. The topographic slope angle is defined as the magnitude of the grid topography gradient vector, and use the arc tangent to convert it to an angle. The topographic ruggedness, which measures the degree of irregularity of the topography, is defined by the topographic ruggedness index (TRI), defined by (Riley et al., 1999):

where i and j are the zonal and meridional grid numbers in the specified domain, respectively, and x i j is the elevation of each neighbor cell relative to the center point cell, x 00 . The TRI presents the sum of changes in elevation between a grid cell and its neighboring cells and is equivalent to the standard deviation in two dimensions (Riley et al., 1999). In this study, we calculate the topographic slope and topographic ruggedness on the same 1°×1° spacial domain as that of the bin averaged distribution of H_BML.

The H_BML computed from the 514 CTD profiles range from 4~255 m, with the mean, median, and one standard deviation being 73 m, 56 m, and 55 m, respectively ( Figure 3 ). The mean and median values are much smaller than the mean H_BML (154 m) that estimated from the 201 full-depth hydrographic profiles by Li et al. (2022), but somewhat larger than the median value (47 m) in the global ocean estimated with full-depth CTD data from the World Ocean Circulation Experiment program (Huang et al., 2019). The probability density distribution of the H_BML ( Figure 3 ) demonstrates that 45% of the H_BML are in the range of 20~80 m, and 27% of the H_BML are thicker than 100 m, with a positive skewness (1.25).

The H_BML shows different distribution characteristics with water depth. The thickest averaged H_BML (133 m), with a standard deviation of 72 m, occurs at a depth of ~700 m in the SCS ( Figure 4A ). The correlation between the H_BML and water depth is positive for water depths shallower than 700 m, and negative for water depths deeper than 700 m. These relationships are quite different from the global distribution presented by Huang et al. (2019), where the H_BML increased exponentially for water depths deeper than 1000 m (the statistics for stations shallower than 1000 m were not presented and stations shallower than 500 m were not considered in their study). These differences between the SCS and the open ocean suggest that the distribution of the H_BML in the SCS may be regulated by more local dynamic factors rather than the general factors in the open ocean.

To obtain a quantitative relationship between H_BML and water depth, we calculate the ratio (R_H/D) between the H_BML and the total water depth (D) to evaluate how the H_BML varies as a function of depth ( Figure 4B ). In general, the R_H/D decreases roughly exponentially with water depth following the least-squares best fit curve:

The mean ratio ranges from 10–50% for depths between ~100–700 m, between 5–10% for depths of ~1000 m, and less than 2% for water deeper than 3000 m. All these values are slightly higher than the results estimated in the North Atlantic (Lozovatsky et al., 2008; Lozovatsky and Shapovalov, 2012), suggesting that bottom mixing is more active in the SCS.

To estimate the seasonal variation of H_BML in the SCS, we calculate the median values of the H_BML and its standard deviations in each season. The results show that the H_BML is relatively small in summer, and large in spring and autumn, suggesting the H_BML has obvious seasonal variation in the SCS ( Figure 5 ). Unfortunately, the feature of H_BML in the winter season cannot be described because only few data were collected during the winter.

To investigate the spatial distribution of H_BML in the SCS, we project the geographic scatter data in the 1°×1° bins and calculate the median values of H_BML in each grid cell ( Figure 6 ). The results show that the H_BML is thicker (>100 m) on the northern continental slope of SCS, especially in the regions adjacent to the west Luzon Strait and Dongsha Islands, where previous studies have suggested that bottom mixing is enhanced (Tian et al., 2009; Yang et al., 2016; Lu et al., 2021). In contrast, the H_BML is thin over both the continental shelves (around 30–60 m) and the deep-sea regions (around 10–50 m). Despite the limited data in each bin, the variations in the H_BML (calculated by using at least 5 points in the bin) is relatively large on the northern continental slope (not shown), indicating a relatively unstable H_BML over the continental slope regions.

Two deep moorings were deployed to measure the variability of the near-bottom potential temperature over the continental slope (M1) and deep-sea region (M2). Figure 7 shows the temporal variations of the bottom potential temperature and its corresponding gradient. The potential temperature and its variations at the M1 were much larger than that at M2. Their mean potential temperature with the one standard deviation was 2.31 ± 0.04°C and 2.07 ± 0.01°C, respectively. In particular, the mean potential temperature lapse rates (analogous to the “lapse rate” in the atmosphere) calculated from the two moorings were -5.75×10^-4 °C/m and -3.21×10^-4 °C/m, respectively.

Although the differences in the potential temperature profiles between M1 and M2 are significant, quasi-homogeneous layer structures in the near-bottom zone can be clearly seen. During most of the observation period, the homogeneous layer is much thicker at the M1 mooring than that at the M2 mooring. The height of the homogeneous layer appears roughly between 100~120 m above the bottom at M1 and 40~60 m above the bottom at M2. This confirms the earlier observation that the H_BML over the continental slope is thicker than it is in the deep-sea regions. In addition, the structures of the quasi-homogeneous layer at M1 are more complex than that at M2, suggesting that the quasi-homogeneous layer is much unstable at M1 than that at M2. The result is also consistent with the CTD observations that the H_BML over the continental slope is unstable than that in the deep-sea regions.

In the section, we combine the observational results to analyze the potential formation mechanisms of the BML in the SCS, especially focusing on the BML differences between the northern continental slope (M1) and the deep-sea regions (M2). Although it is unclear whether the BML can be regarded as the classic bottom Ekman Layer (Armi and Millard, 1976; Beaulieu and Baldwin, 1998), the velocity shear within the bottom Ekman layer drives stable mixing that keeps the layer unstratified. In bottom Ekman layer dynamics, the turbulent Ekman layer formed above the ocean’s bottom can be expressed as h ≈ 0.4 u * / f , where the u* is friction velocity and f is the Coriolis parameter. In practice, the friction velocity u*, which correlates with bottom currents and bottom drag coefficient (C_d ), is not clearly determined from the observations. One simplest procedure to reinterpret the Ekman layer height (h) with friction velocity u* depend on the local mean current and constant drag coefficients (Stahr and Sanford, 1999). Because few current measurements are available, the fluctuations of bottom temperature linked with bottom currents are used to analyze the mixing strength and possible energy sources that influence the structure of the BML.

Spectral analysis of the bottom potential temperatures (using hourly data) shows that variability over the continental slope was dominated by the internal tidal and near-inertial signals, while only low-frequency oscillations (~60 days) were significant in the deep ocean ( Figure 8 ). It is suggested that the differences in the H_BML between the continental slope and deep-sea regions may be caused by different dynamical processes. In the SCS, the internal tides are widely distributed on the continental slope (Alford et al., 2015; Wang et al., 2016), with most of those signals emanating from the Luzon Strait (Zhao, 2014; Alford et al., 2015). The near-inertial signals are likely injected into the upper ocean by typhoon processes (Xu et al., 2013; Zhang et al., 2022) or by other upper ocean mechanisms (Alford et al., 2016). The dominance of these two signals suggests that the internal tidal and near-inertial motions may play the primary roles in the strong mixing along continental slopes (Tian et al., 2009; Wang et al., 2016; Lu et al., 2021). For the deep ocean, the low-frequency oscillations (30-120 days) may be induced by the topographic Rossby waves or the deep ocean eddies, which has been shown in recent observational studies (Zhou et al., 2017; Zhou et al., 2020; Zheng et al., 2021).

To confirm whether the tidally driven mixing and dissipation match the BML distribution pattern in the SCS, we calculated the distribution of baroclinic energy conversion (E), turbulent dissipation rate (ϵ), and diapycnal diffusivity (k_v ) in the SCS (see Section 2.3). As expected, the spatial distributions of the depth-integrated E, ϵ, and k_v reveal significant spatial variations ( Figures 9A–C ). The E, ϵ, and k_v near the Luzon straits and the northern continental slopes usually can exceed 1×10^-3 W m^-1, 1×10^-7 W kg^-1, 1×10^-3 m² s^-1, respectively. The regions with these high values correspond to the areas in which the H_BML is thicker in the SCS ( Figure 6 ). Whereas except some regions with rough topography such as islands, seamounts, and shelf breaks the E, ϵ, and k_v in the central SCS are about one or two orders of magnitude weaker, which corresponds to the thinner H_BML in the SCS ( Figure 6 ). The above typical characteristics of mixing and dissipation pattern match well with the internal tide dissipation as estimated by de Lavergne et al. (2019) ( Figures 9D , 10 ), as well as consistent with the previous observational reports (Tian et al., 2009; Yang et al., 2016; Shang et al., 2017; Lu et al., 2021) and model estimations (Wang et al., 2016). It suggests that the tidal energy is a major source of energy to inducing the strong dissipation and diffusivity, thus influencing the distribution of BML in the SCS.

In addition, the above results reveal the diapycnal diffusivity at M1 is higher than that at M2, which can also be confirmed by the mooring observations (see Section 2.5). The calculation shows that the estimated average eddy diffusion coefficient varies between 1.0×10^-6 to 5.3×10^-3 m² s^-1 at M1 and -2.3×10^-7 to 5×10^-3 m² s^-1 at M2. The estimated mean eddy diffusion coefficients were 1.3×10^-3 m s^-1 and 8.4×10^-4 m² s^-1 at M1 and M2, respectively. These values are in agreement with practical observations estimated from the Thorpe-scale method and direct observations in previous studies (Yang et al., 2016; Shang et al., 2017; Lu et al., 2021), suggesting that the bottom vertical mixing over the continental slope (M1) is stronger than that in the deep-sea regions (M2), which may explain why the H_BML over the continental slope is thicker than in the deep-sea regions.

We analyze other potential factors that may contribute to the bottom mixing and their distribution differences in the SCS. CTD and mooring observations are used to consider the roles of topography, internal tidal dissipation, and density stratification in the distribution of H_BML between the continental slope (M1) and deep-sea regions (M2). The results in Figures 11A–C show that the H_BML increases with increasing topographic slope, topographic ruggedness, and internal tidal dissipation. The typical values of the topographic slope, topographic ruggedness, and internal tidal dissipation are all larger at M1 than that at M2, and likely all contribute to the thicker H_BML over the continental slope. Figure 11D shows the H_BML as a function of the buoyancy frequency estimated within the BML. This result shows no monotonic relationship between the H_BML and the buoyancy frequency as suggested by Huang et al. (2019), in which they showed that the H_BML tends to be thinner with stronger stratification. However, our results suggest that the stratification may not be important for determining the BML differences between the continental slope and the deep-sea regions. In other words, the dominant factors controlling the distribution of the H_BML in the SCS are dynamic rather than thermodynamic.

Taken together, we conclude that the dynamic processes over the northern continental slope controlling the BML are the energetic, high-frequency forcing, together with the large slope and steep topography, which combine to cause strong tidal energy dissipation and vertical mixing near the bottom in these regions. As a result, the BML on the northern continental slope is relatively thick. Conversely, in the deep-sea regions, the dynamic processes are low-frequency, and the topographic roughness and slope are relatively smooth and gentle, so that the tidal energy dissipation and bottom vertical mixing are considerably weaker, resulting in a relatively thin BML in the deep-sea regions.