The level-2 sea surface temperature (SST) product measured by the Visible and Infrared Imager/Radiometer Suite (VIIRS) sensor of the Suomi-NPP satellite from June to July 2020 is used to validate the numerical model. It has a spatial resolution of 0.75 km and a temporal resolution of 1 h. The access website is https://ncc.nesdis.noaa.gov/NOAA-20/NOAA20VIIRS.php.

Another SST product is obtained from the Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) system, which has a 5-km spatial resolution and daily temporal resolution (https://ghrsst-pp.metoffice.gov.uk/ostia-website/index.html). The products are merged with microwave and infrared satellite SSTs provided by the Group for the High-Resolution Sea Surface Temperature project (Martin et al., 2012).

Daily sea-level anomaly (SLA) with a spatial resolution of 0.25° is obtained from the Archiving, Validation, and Interpretation of Satellite Oceanographic data (AVISO) project, which integrates JASON-1, TOPEX/Poseidon, and ERS-1/ERS-2 products (https://www.aviso.altimetry.fr/en/data/data-access.html).

The level-2 chlorophyll-a concentration products are obtained from the VIIRS onboard the JPSS-1. The spatial resolution of these data is about 0.75 km, and the temporal resolution is several hours.

We set up the Massachusetts Institute of Technology General Circulation Model (MITgcm) to investigate submesoscale structures. To reduce the influence of the open boundaries in the simulation area as much as possible and to improve the computational efficiency, we run MITgcm using a rectangular one-way nested grid with two levels. MITgcm adopts the z-coordinate and static approximation because, for submesoscale flow of R o ≈ 1 , the vertical scale in the upper ocean is much smaller than the horizontal scale, that is, H / L  ≪ 1 .

The bathymetric chart is derived from ETOPO1 provided by the National Oceanic and Atmospheric Administration, which is an arc-minute global relief model of the Earth’s surface. As shown in Figure 1 , the parent grid named M1 covers a spatial range of 100°E–135°E and 0°N–30°N with a horizontal resolution of 0.25°. The initial conditions of M1 are based on the January climatological data of WOA2018 (https://www.nodc.noaa.gov/cgi-bin/OC5/woa18/woa18.pl). The boundary conditions, sea surface relaxation, and wind field are all interpolated using SODA3 (http://apdrc.soest.hawaii.edu/erddap/search/index.html?page=1&itemsPerPage=1000&searchFor=soda). After 15 years of spin-up, the temperature, salinity, and major currents of the SCS reach statistically steady states.

Based on the results of the M1 climatological state, the nested zone M2 covering 108°E–118°E and 13°N–21°N in the northern SCS is selected for one-way nesting with a horizontal resolution of 0.05°. The M2 boundary conditions are from HYCOM gofs3.1 daily data (https://www.hycom.org/dataserver/gofs-3pt1/analysis), and the initial conditions are from the climatological data of M1. The wind field and other meteorological data for calculating the sea surface heat flux are interpolated from the ERA5 reanalysis data provided by the European Center for Medium-Range Weather Forecasts (https://www.ecmwf.int/en/forecasts/datasets). The simulation period is from 1 October 2017 to 31 July 2020.

The horizontal resolution of M3 (15.5°N–20.5°N, 109.5°E–114.5°E) is refined to 0.01° (about 1 km). The boundary conditions of M3 are provided by the M2 output data, including temperature, salinity, velocity, and sea surface height. The K-profile parameterization (Large et al., 1994) is adopted for parameterizing vertical mixing. The M3 simulation spans from 1 August 2019 to 31 July 2020, with 6-hourly outputs of sea surface height, temperature, salinity, and velocity. The time step is 100 s.

It is necessary to examine the accuracy of the model data before using them. The simulated annual average SST image is compared with OSTIA SST and the simulated annual average velocity image with AVISO products. As shown in Figures 2A, B , the model reproduces the cyclonic circulation in the north of the SCS, and the velocity and range of the SCSWBC have a similar spatial distribution with AVISO. However, the simulated velocity at the southwest boundary and the southeast boundary is different from that of AVISO, which may be due to the influence of the flow correction at the boundary of the model on the simulation results.

The simulated annual average SST is higher in the southeastern corner of the analysis region, which is consistent with the spatial distribution of OSTIA SST ( Figures 2C, D ). A submesoscale front (~10 km) observed by VIIRS satellite data on 6 February 2020 is also captured by MITgcm ( Figures 2E, F ). In general, it is found that the output of MITgcm is reliable and meets the research needs of submesoscale motions.

We also validate the seasonal variations of modeled temperatures and velocities. Figure 3 shows the seasonal variation of temperatures at 5, 55, and 110m. The spatial distribution of temperatures is consistent between SODA3 and MITgcm in all seasons. The fine structures of MITgcm are much clearer than those of SODA3. The mean biases of temperature in spring, summer, autumn, and winter are −0.15°C, 0.10°C, 0.19°C, and 0.18°C, respectively. Figure 4 shows the seasonal variation of the velocities of SODA3 and MITgcm. Although the modeled magnitudes of velocity are lower than the observed value in the four seasons, the main current, i.e., the western boundary current, is significant in both the SODA3 and MITgcm. It indicates our model result is reliable.

The horizontal isotropic wavenumber spectrum is an effective means of estimating the spatial critical values of submesoscale and mesoscale processes (Cao et al., 2019) and is commonly used in the processing of model data (Rocha et al., 2016; Li et al., 2019). Under the assumption of horizontal motion uniformity and isotropy, within the wavenumber range dominated by quasi-geostrophic motion, the slope of the kinetic energy spectrum follows a k − 3 power law (Charney, 1971). In the submesoscale range, the slope of the kinetic energy spectrum is relatively flat and follows a k − 2 power law. Many observational and simulation results have confirmed this conclusion (Callies and Ferrari, 2013; Rocha et al., 2016; Qiu et al., 2017; Cao et al., 2019; Zhang et al., 2020a; Cao et al., 2023). This suggests that by determining the wavenumbers corresponding to the transition between k − 3 and k − 2 in the horizontal wavenumber spectrum, the spatial upper limit of submesoscale processes can be estimated.

The Python program provided by Rocha et al. (2016) is used to estimate and plot the horizontal isotropic wavenumber spectrum. The program first removes the linear trend of the raw data, then performs a two-dimensional discrete Fourier transform on the data and multiplies it by a two-dimensional Hanning window. Subsequently, the power spectral density of the meridional and zonal velocities is obtained separately through spectral estimation methods. Finally, the two-dimensional kinetic energy spectrum is synthesized based on the following formula,

where spec 2 D , spec u and spec v represent the two-dimensional power spectral density and the power spectral density of the meridional and zonal velocities, respectively.

This study area includes both shallow coastal waters and open deep waters, with large differences in topography and water depth. This leads to differences in the Rossby radius of deformation between the shelf and open sea areas, and a single spatial critical value cannot be used to extract submesoscale processes in both areas, which must be discussed separately. This study uses the location of the 210-m isobath as the boundary between the shelf and open water areas.

We estimate the horizontal isotropic wavenumber spectrum of the S1 region on 2 January 2020, and the B1 region on 5 February 2020 ( Figures 5 , 6 ). The slope of the kinetic energy spectrum in the shelf area is close to k − 3 in the mesoscale range but gradually flattens with decreasing wavelength in the range of wavelengths less than 14 km, exhibiting a k − 2 power law ( Figure 6A ), indicating that submesoscale motion is active in this wavelength range. The slope of the horizontal isotropic wavenumber spectrum in the open water area changes near a wavelength of 30 km and is close to k − 2 in the submesoscale range (<30 km) ( Figure 6B ). Therefore, this study selects 14 and 30 km as the spatial critical values for mesoscale and submesoscale processes in the shelf and open water areas, respectively.

Extracting submesoscale processes from full-spectrum data is a key step in submesoscale analysis (Wang et al., 2022). The Eulerian-based filtering method has been widely used for separating submesoscale processes, which includes temporal filtering (Barkan et al., 2017; Zhang et al., 2020a), spatial filtering (Dong et al., 2022; Song et al., 2022; Su et al., 2020b), and spatiotemporal filtering (Capet et al., 2008b; Qiu et al., 2018; Torres et al., 2018). In addition to the commonly used Eulerian filtering methods, many scholars have proposed different methods to improve the accuracy of submesoscale decomposition. For example, Wang et al. (2022) used the Lagrangian filtering method to remove the influence of mesoscale background flow, focusing only on the submesoscale VHT under the mesoscale eddy background. Wang et al. (2023a) proposed a dynamic decomposition method for ocean multiscale motion and successfully applied it to multiscale motion decomposition in the central SCS (Wang et al., 2023b).

This work employs a spatiotemporal filtering method to extract submesoscale processes. We adopt the Reynolds decomposition method to decompose the original physical quantities into annual average terms and perturbation terms. The perturbation terms are then decomposed into mesoscale and submesoscale components. The calculation formula is as follows:

where A denotes the physical quantity to be decomposed (such as w , T , U , V , etc.); A ¯ represents the annual mean result of the physical quantity; A m and A s represent the mesoscale and submesoscale components of the physical quantity, respectively; and A ˜ is the spatiotemporal mean of the physical quantity A . In this study, 7 days are selected as the time-critical value for mesoscale and submesoscale processes (Torres et al., 2018; Wang et al., 2022), and 30 and 14 km are selected as the spatial critical values for open water and shelf area, respectively. A ˜ is obtained by spatially averaging the physical quantity A within a specified range (0.3° × 0.3° for open water, and 0.14° × 0.14° for shelf area) after 7-day averaging. As an example, the decomposition results of the vertical velocity w at 12:00 on 1 January 2020 are shown in Figure 7 .

The submesoscale VHT is defined using the following formula:

where C p is the specific heat capacity of seawater, taken as 4,096 J · kg − 1 · ° C − 1 ; ρ is the density of seawater, taken as 1,024 kg · m − 3 ; and w s   and  T s are the submesoscale components of vertical velocity and water temperature, respectively. Similarly, the mesoscale vertical heat transport can also be obtained from the following formula:

where w m and T m are the mesoscale components of vertical velocity and temperature, respectively.