SCHISM is a sophisticated three-dimensional finite-element model for simulating fluid and ecological dynamics within the realm of oceanography (Zhang et al., 2015, 2016). Operating on unstructured grids (UG), the model boasts several unique features: (1) Utilization of unstructured hybrid triangular/quadrilateral grids in the horizontal dimension and option between hybrid SZ coordinates (hybrid S and shaved z coordinates) or the innovative LSC2 layering scheme in the vertical dimension; (2) integration of a cutting-edge high-order implicit convection transport scheme (TVD2), ensuring both mass conservation and monotonicity; (3) high-order transport capabilities in the horizontal dimension, employing third-order weighted essentially non-oscillatory (WENO) formalism (Ye et al., 2019), alongside a novel momentum equation convection scheme that includes optional high-order Kriging and explicit local adaptive dissipation filters; (4) incorporation of new horizontal viscosity schemes, including biharmonic viscosity, adept at effectively filtering out spurious inertial modes without inducing excessive dissipation; (5) commendable tolerance for grid quality in non-tidal areas, eliminating the necessity for depth smoothing. In essence, this model represents a fusion of efficient and accurate semi-implicit finite-element and finite-volume methods. The model also incorporates the generic length scale (GLS) proposed by Umlauf and Burchard (2003) to evaluate the eddy coefficient. The advantage of the GLS is that it includes most terms from the closure schemes of two-equation models. Augmented by the Euler–Lagrange algorithm, SCHISM addresses the hydrostatic form of the Navier–Stokes equations. This comprehensive approach enables the model to simulate a wide spectrum of physical and biochemical processes within water bodies, making it a versatile tool for advancing oceanographic research.

The model domain coverage is shown in Figure 1A , where the grids in the west encompass the entire Beibu Gulf and extend south to the basin of the northern South China Sea (NSCS). The eastern boundary extends to the Taiwan and Luzon Straits and encompasses the Pearl River Estuary. The horizontal grid comprises 161,277 nodes and 316,763 triangular elements. The coastline data utilised in the model are derived from the global self-consistent, hierarchical, high-resolution geography database (GSHHG, https://www.ngdc.noaa.gov/mgg/shorelines/shorelines.html), with specific modifications introduced using electronic nautical charts for the NSCS. Regarding bathymetry data, ETOPO 2022 is employed with a resolution of 15 arcsec (approximately 460 m, https://www.ncei.noaa.gov/products/etopo-global-relief-model). The maximum water depth is approximately 4,800 m, in the south-eastern part of the computational domain. The grid resolution increases from the south-east boundary towards Beibu Gulf, maintaining a resolution of approximately 18-km in the open sea. Local refinement is implemented with a resolution of approximately 8 km to address significant depth variations on the north-western slope of the SCS. Within Beibu Gulf, the resolution narrows to 2–4 km and further increases to 500 m near the QS. In the vertical direction, this study adopts a flexible vertical layering approach (local coordinates) previously utilised by Zhang et al. (2015) and Yu et al. (2017). The main vertical grid is configured with an average of 47 layers ( Figure 1B ), including 30 layers of master grids in regions shallower than 100 m and a maximum of 77 layers, providing a comprehensive representation of the vertical structure of the water column.

The tidal current and elevation forcing at the open boundary are generated using FES2014 (https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes.html), encompassing 14 tidal components (M2, K1, O1, S2, Q1, P1, N2, K2, M4, M6, MS4, SSA, MM, and MF) driven by harmonic constants. The model produces hourly wide-area water level results, and harmonic analysis using the T_Tide software package (Pawlowicz et al., 2002) extracts the cotidal charts for diurnal tides, represented by K1, and semi-diurnal tides, represented by M2. The amplitude and phase distributions of these tidal constituents are shown in Figure 2 . The tidal constituents in FES2014—a numerical outcome derived from a global ocean tidal finite element model that exhibited significantly enhanced accuracy compared with its predecessors, namely FES2004 and FES2012—are used as reference for this analysis. The amplitude and phase distributions of both the diurnal and semidiurnal tides in the model are consistent with those of FES2014 and align with the findings of previous studies (Xu et al., 2010; Minh et al., 2014; Ding et al., 2017). Specifically, the M2 tide reveals a degenerate counterclockwise amphidromic point in the north-west of Beibu Gulf, with an amplitude of approximately 40 cm, and the amplitude increases towards the Leizhou Peninsula in the eastern part of the QS. The K1 tide features an amphidromic point near Hue (Vietnam) and achieves an amplitude of 90 cm.

The residual current and water level at the open boundary are generated using the HYCOM + NCODA global analysis data (https://www.hycom.org/). These data not only establish the initial conditions for the model but also establish a ‘hot start’, incorporating sea surface height, temperature, salinity, and flow velocity data. Atmospheric forcing at the sea surface is based on ERA5 reanalysis data (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=form) from the European Centre for Medium-Range Weather Forecasts (ECMWF). Specific details on the HYCOM and ERA5 data are presented in Table 1 , which provides a comprehensive overview of the atmospheric and oceanic input parameters used to drive and initialise SCHISM for accurate simulations in the QS.

Surface elevations of 16 tidal gauge stations are selected for analysis ( Figure 1A ) using data from the Global Tidal Data Platform (http://global-tide.nmdis.org.cn/). The validation period spans from January to March 2021, and the results, as illustrated in Figure 3 , reveal robust fitting outcomes at all 16 stations, with errors being mostly less than 10 cm. Notably, Shantou station at the model boundary exhibits larger errors, which are attributed to boundary effects, whereas station 7 in Haikou may present differences due to depth interpolation or coastline resolution issues (Ding et al., 2017). Overall, the model establishes its credibility by accurately simulating water level variations.

Verification of the surface flow velocity and direction in the QS was conducted by comparison with data obtained from the forecasted values of the Global Tidal Data Platform for a selected time period in March 2021. The distribution of the monitoring stations is shown in Figure 1A , and the validation results are shown in Figure 4 . The modelled flow velocity generally aligns with the forecasted values both in magnitude and trend. However, some discrepancies in details are evident; these are primarily attributed to the limited selection of the two characteristic angles in the forecasted values. This limitation can result in significant angular differences during flow reversal. Despite these discrepancies, the good fit of the long time series data highlights the robustness of the model and its accurate description of the flow field in the QS.

The distributions of tidal amplitudes and phase lags in the QS are shown in Figure 5 , where the diurnal tide is characterised by K1, and the semidiurnal tide is represented by M2. The semidiurnal tidal waves exhibit a general west-to-east propagation through the strait, featuring a quasi-amphidromic point at the south-western corner of the Leizhou Peninsula. The amplitude gradually increases along the direction of propagation, ranging from approximately 10 cm at the western entrance to 50 cm at the eastern entrance. The tidal wave propagation speed can be estimated from the density of the co-phase lines. Notably, the speed of semi-diurnal tides rapidly increases after passing through the QS, whereas diurnal tides exhibit a different propagation direction. Co-phase lines for diurnal tides indicate propagation from east to west across the QS, with a noticeable decrease in speed upon entering the strait. A quasi-amphidromic point emerges at the north-eastern corner of Hainan Island. Post-modulation, the distribution of the co-phase lines becomes nearly parallel to the northern and southern coasts. Another diurnal tide enters from the southern mouth of Beibu Gulf, reaches the west side of the strait, and propagates from west to east at a speed similar to that inside the strait. Convergence of these two diurnal tides is observed north-west of Hainan Island.

The amplitudes of K1 increase with the phase lag, with minimum values at the eastern mouth ranging from 10 to 20 cm and maximum values at the western mouth reaching 70–80 cm. These overall characteristics align with the results of the numerical simulations of Shi et al. (2002). However, the numerical values and intervals of equal-amplitude lines and equal-phase lag lines within the strait and on the eastern side may vary owing to changes in topography and the influence of model/interpolation accuracy.

This study also considered the computational methods of Wu et al. (2020) and Dai et al. (2017) to provide the tidal flux distribution characteristics in the QS in the form of tidal constituents, as shown in Figure 6 . The calculation formula is as Equation 1:

where ( Φ x , Φ y ) represent the components of tidal energy flux density in the ( x , y ) directions, t is time, T is the tidal period, h is water depth, ρ is the density of seawater, which is assumed to be 1,025 kg/m³ in this study, g is the acceleration due to gravity, which is taken as 9.8 m/s² in this study, ζ ˜ represents the elevation of the sea level (in meters), ( u ˜ , v ˜ ) represent the vertical average flow velocities in the ( x , y ) directions (m/s), Z and θ represent the amplitude (m/s) and phase lag (°) of the tide, respectively, ( U , V ) are the amplitudes of the tidal flow components in the ( x , y ) directions (m/s), and ( ξ , η ) represent the corresponding phase lags (°).

The overall characteristics of the tides in the QS presented in this study align with those described by Wang and Wang (2009) for the diurnal and semi-diurnal tides (represented by K1 and M2, respectively) and Zhao et al. (2010). However, our study provides a more detailed distribution in the QS. The tidal waves at the eastern and western entrances of the strait predominantly originate from the NSCS, propagating westward and bifurcating into two main branches that enter Beibu Gulf and the western coast of Guangdong Province (Shi et al., 2002).

In the semi-diurnal tide, the M2 tidal fluxes are relatively large, entering the strait from the western entrance and reaching their maximum value in the central part of the strait, i.e. up to 2×10⁴ kW/m. The K1 tidal fluxes further illustrate the characteristics of tidal waves propagating in opposite directions on both sides of the strait. Energy is transmitted from both the east and west sides into the interior of the strait with magnitudes close to those of the semi-diurnal tides. The transition zone of the interaction is located approximately along the Mulan Bay–Beila Port line, as observed from the distribution of tidal amplitudes and phase lags in Figure 5 .

After passing through this transition zone, the propagation speed of the tidal waves rapidly decreases. This suggests that although the tidal waves from the east and west sides intersect north-west of Hainan Island based on the lines of equal phase lags, this may already be a result of the superposition of waves from both directions. From the perspective of tidal energy propagation, it is likely that the tidal waves on both sides of the strait meet at its eastern entrance. The current numerical model results further support the existence of local tidal wave resonance in the QS, encompassing interactions among tidal waves as well as between tidal waves and terrain (Shi et al., 2011; Ding et al., 2013; Zhu et al., 2014; Cheng et al., 2017; Li, 2022).

The distribution and variations of the background residual flow field in the QS were examined through numerical simulations, offering insights into the seasonal and interannual variations of the water transport. The results shown in Figures 7 , 8 were obtained separately for 2021, 2022, and 2023, serving as reference experiments for the exploration of the residual flow field and subsequent mechanism. Monthly averages were used to derive the residual flow distribution for each month. Unlike longer-term seasonal averaging, which might smooth the overall patterns, this study opted for monthly averages to capture more pronounced variations and achieve a more effective comparative analysis. The results for February represent the winter residual flow, whereas those for July represent the summer residual flow.

Figure 7 reveals that the winter surface residual current in the QS generally flows westward, reaching maximum flow velocities of 0.4–0.5 m/s. The distribution of flow velocity correlates with water depth. Near the prominent dividing ridge on the west side of the strait, there is a high-velocity zone that increases with water depth, which is particularly evident in winter 2022 ( Figure 7B ). Additionally, irregular flow velocity distributions or eddy-producing areas are evident near the eastern entrance of the strait because of the complex terrain there. This is consistent with the findings of Shi et al. (2002) and Chen et al. (2009). Importantly, these areas do not significantly impact the residual currents within the strait and tend to disappear with the smoothing of the terrain. The distribution of the flow vectors reveals that the background flow enters from the north-eastern entrance of the strait and crosses the channel along the Leizhou Peninsula into Beibu Gulf. The flow velocity in the middle of the strait exhibits a north-high and south-low distribution; this pattern occurs also in summer ( Figure 8 ).

The summer surface residual current in the QS ( Figure 8 ) generally exhibits weaker velocities compared to those of the winter surface residual current, with the maximum residual current velocities reaching 0.2 m/s in certain areas near specific topographical features. The high-velocity zone in the central–western region of the strait resembles the distribution of the K1 tidal energy flux ( Figure 6A ). The influence of topography on current distribution is more pronounced. Numerical simulation results for 2021–2023 indicate that the QS experiences westward transport during summer, which is consistent with previous studies (Shi et al., 2002; Zu et al., 2005; Dasen, 2006; Chen et al., 2009; Zhu et al., 2014; Zheng, 2015; Ding et al., 2017).

Nevertheless, interesting features appear especially in summer 2022 ( Figure 8B ). During this period, a distinct eastward flow is evident for the surface residual current in the central part of the strait. This contrasts with the significant westward flow observed in summer 2021. Monthly eastward transport from the northern part of Haikou to the northern part of Mulan Bay is evident, accompanied by a noticeable weakening of the strong westward flow along the northern coast of the Leizhou Peninsula. This implies the existence of another form of eastward transport in the QS during summer, besides the diurnal oscillation characteristics proposed by Zhang et al. (2009), which manifest as subtidal oscillations in this study.

The observed eastward transport in summer 2022 may indicate interannual variability or influence of larger-scale factors, offering a new perspective on the contradictory findings of previous studies. This highlights the dynamic nature of the residual currents in QS, which are subject to varying conditions and influencing factors over time.

In summer 2021, the large-scale residual flow results presented in Figures 8 , 9 suggest a prevailing westward flow in the QS from the surface to the bottom. In the profile ( Figure 9A ), there are two westward transport flow nuclei in the northern and southern parts, whereas the central region exhibits weak residual flow. In contrast, the results of summer 2022 indicate an eastward flow trend in the surface layer of the QS, with strong eastward transport observed in the water column. The eastward flow is concentrated in the middle part of the strait (near 20.16° N), reaching a maximum intensity in the surface layer in excess of 0.05–0.10 m/s and extending to a depth of 20 m, gradually decreasing with increasing water depth, while the northern and southern parts of the QS continue to exhibit smaller westward transport flow nuclei. This pattern is in line with the overall structure revealed by the measurements of Wang et al. (2014) along 110.18° E. The results for 2023 are almost the same as the average ( Figures 8C , 9C, D ), showing a transitional pattern between 2021 and 2022. In the middle part of the QS, the surface water exhibits weak eastward flow with an intensity of less than 0.05 m/s and a depth of about 10 m.

Comparatively, the simulated results for 2022 reveal more pronounced eastward flow features in both the horizontal and vertical profiles. Hence, 2022 is considered as a “significant eastward summer flow year” of the period of 2021–2023. Conversely, 2021 is considered as a “significant westward summer flow year.” Finally, 2023 is considered an intermediate year for comparison.

To scrutinise the intricacies of water transport processes in the QS, this study examines local factors. Three characteristic points, namely st1–st3, are selected to represent the western, eastern, and central parts of the strait, corresponding to areas of medium outflow (st1), strong inflow (st3), and weak flow (st2), respectively ( Figure 7A ). Momentum balance analysis is performed in both the x-direction (parallel to the banks) and y-direction (perpendicular to the banks) within the interior of the QS, while considering the motion equation. The methodologies of Shi et al. (2002); Jiang (2018); Wang et al. (2021), and Hench et al. (2002) were comprehensively utilised, leading to the derivation of horizontal momentum equations in both directions (as shown in Equation 2 and Equation 3):

In the equation, f is the coefficient of Coriolis parameter, f u , f v are the Coriolis force terms, ρ 0 is the density of seawater, τ is the wind stress generated by the 10-m wind at the sea surface, p is the pressure of seawater, υ is the kinematic viscosity coefficient, C D is the bottom friction coefficient, Δ = ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 represents the horizontal Laplacian operator, H ( x , y , t ) = h ( x , y ) + ζ ( x , y , t ) is the water depth at a certain moment, and ζ ( x , y , t ) is the sea surface height at a certain moment. When seawater is approximately homogeneous in density, the pressure gradient terms 1 ρ ∂ p ∂ x and 1 ρ ∂ p ∂ y and the related variations in sea surface height g ∂ ζ ∂ x and g ∂ ζ ∂ y can be approximated. The results are shown in Figure 10 , where A C C , H A D V , W i n d s t r , V I S C and B S T R represent the local acceleration, non-linear advection, wind stress, turbulent viscosity, and bottom friction terms, respectively. All terms were filtered for 48 h. Because this study mainly focuses on the east–west flow, in the following, only the equation with the term f u (i.e. across the direction of the QS) is shown, with results corresponding to various characteristic points from 2021 to 2023. When f u (red line in Figure 10 ) is positive, it represents an eastward flow at the characteristic point, whereas negative values indicate a westward flow. However, the horizontal pressure gradient in the east–west direction also affects the u-component flow. In this study, it is represented by the difference in sea level in the area adjacent to the QS, as shown in Figure 1A , by calculating the difference between the mean values of Area1 and Area2: Δ E l e v ( t ) = E l e v a r e a 1 ( t ) − E l e v a r e a 2 ( t ) .

In Figure 10 , the Coriolis force and pressure gradient terms at st1–st3 are significantly larger in magnitude than the other factors, exhibiting quasi-synchronous variations. This suggests that the QS is in quasi-geostrophic balance, in accordance with the conclusions of Shi et al. (2002). Furthermore, the non-linear advection term is prominently enhanced only in winter and is generally much smaller than the Coriolis force term in spring and summer.

Analysing the interannual variations at the three characteristic points in 2021 (i.e. significant westward summer flow year; Figures 10A–C ) shows that all points exhibit strong westward transport in winter, followed by a robust eastward oscillatory flow in spring (Zhang et al., 2009), and a return to westward flow in summer. In 2022 (i.e. significant eastward summer flow year; Figures 10D–F ), there is almost no eastward oscillatory flow in spring, and the strait maintains westward transport from January onwards. However, it gradually becomes positive and shows eastward flow in the form of subtidal oscillations starting in May and is most pronounced from late June to early August. As mentioned earlier, although tidal effects may not have been completely filtered out in the results, studies have suggested that tidal effects in the QS drive the westward flow (Shi et al., 2002). Therefore, if tidal effects are entirely removed, the eastward flow may appear in a more stable form. In 2023 (intermediate year; Figures 10G–I ), the results resemble those in 2021; however, the eastward oscillatory flow in spring 2023 appears later than that in 2021, and despite the eastward transport in summer 2023, its intensity is much weaker than that in summer 2022.

Although the Coriolis force and pressure gradient terms display quasi-synchronous variations, the current experiment does not definitively establish that the pressure gradient perpendicular to the direction of the QS determines the east–west flow. Instead, these two factors appear to complement each other. Compared to other local variables, the pressure gradient ( Δ E l e v ) in the east–west direction may have a more direct impact. In 2021, the water level difference between the western and eastern sides of the strait gradually increased with continuous westward transport, maintaining a west-high–east-low water level state in summer, with an average difference of less than 0.1 m. In 2022, water level accumulation on the west side of the QS occurred later, and the accumulation rate increased rapidly in summer, with an average difference of approximately 0.15 m. The water level difference peaked in late June and was accompanied by the strongest eastward transport. The water level difference in 2023 was average, with values similar to those in 2021, albeit with more obvious fluctuations.

Based on a comparison of several momentum balance results, it can be generally concluded that the average water level difference in the east–west direction of the QS has a significant impact on its flow direction; however, different behaviours may occur in different seasons. In late winter and early spring, a small, stable Δ E l e v (stable only at 5 cm) may cause subtidal oscillations and a compensatory eastward flow. In summer, a larger Δ E l e v is required to trigger a stable eastward flow mechanism. In the experiments conducted thus far, both in spring and summer, the strait has a positive water level difference (i.e. a higher water level on the west side); however, the flow direction of the strait remains inconsistent. This usually means that the critical value for the east–west flow in summer is not near zero. According to the previous discussion, this value is estimated to be approximately 10 cm (although it may vary at different locations within the QS). This phenomenon may be caused by density gradient variations and may be related to the influence of the current west of Guangdong Province and other NSCS shelf-slope current systems on the east side of the QS (Xue et al., 2004; Fang et al., 2012; Shu et al., 2018).

The interannual variations mentioned in Section 4.1 prompt the investigation of external dynamics that transform Δ E l e v into an “eastward transport promoting state.” Besides the potential influence of the NSCS current system, Δ E l e v is likely driven also by interactions at the air–sea interface. Previous studies have shown a close relationship between the eastward flow in the QS and the south-western monsoon during summer. However, momentum balance analysis revealed that the local north–south wind magnitude is only about one-tenth of the Coriolis force term. Therefore, this study further emphasises its potential role by adopting the vorticity balance method of Gao et al. (2015) for Beibu Gulf. The monthly average wind stress vorticity from February to July in 2021–2023 was calculated in the vicinity of the strait, as shown in Figures 11 , 12 , where the vectors represent the monthly average wind field.

From late winter to spring, the wind stress curl along the QS exhibits a basic symmetric distribution. In February 2022 ( Figure 11B ), there is a strong positive wind stress curl on the western side of the strait, corresponding to a positive input of vorticity (Gao et al., 2015). Furthermore, a positive vorticity causes surface divergence, whereas a negative vorticity curl produces the opposite effect. This leads to a decrease in sea surface divergence and intensifies the westward transport during winter. The wind stress curl causes a certain degree of sea level decrease, which may explain why Δ E l e v in 2022 and 2023 exhibits a decrease or even a negative phase during late winter and early spring.

Conversely, Figure 13 shows the relationship between the Coriolis force term ( f u ) at different locations and the variation in the east–west-wind (U-wind) speed. During spring, there is little change in the U-wind direction, and intermittent eastward flow occurs more frequently during the weakening process or around zero values of the long-term easterly wind (as indicated by the orange box), with the situation being more pronounced in 2022 and 2023. The subtidal-period eastward flow tends to be compensated by the relaxation process of the wind, and there is a short lag in the flow direction change relative to the wind direction change. When the easterly wind strengthens again, f u quickly trends towards a negative phase, in accordance with the results of Wang et al. (2021) for the Taiwan Strait. The period of the subtidal oscillations is likely related to the tidal period, as determined by comparing it with the local tidal cycles ( Figure 3 ). However, it is difficult to explain the complete phase reversal between f u and the U-wind direction in summer, as well as the variation in the strength of the eastward flow (indicated by the purple box), indicating that the wind-driven influence in summer is dominated by other forms.

In Figure 12 , the difference in wind stress curl between the two sides of the strait reappears. A positive wind stress curl in the east causes surface divergence, whereas a negative wind stress curl in the west causes surface convergence, leading to a further increase in Δ E l e v . In June and July, the region transitions to the south-western monsoon, while near the QS, the wind is mainly southerly, which has a limited impact on the east–west flow. Wind mainly causes water level changes through the wind stress curl, and the interannual differences in the wind field also have an impact on the eastward flow, mainly concentrated in the eastern region of the strait. According to Gao et al. (2015), the effect of wind stress curl is more pronounced in areas with weaker exchange with the open ocean. Therefore, the positive wind stress curl area in the eastern part of the strait in 2021 ( Figure 12D ) being much larger than those during the same periods in 2022 and 2023 ( Figures 12E, F ) does not imply stronger seawater divergence. In contrast, the broad input of positive wind stress curl may further intensify the current on the west coast of Guangdong, leading to water accumulation on the east side. This also explains why, although the difference in wind stress curl on both sides of the strait is similar, there is no eastward flow in summer 2021, which has the widest positive wind stress curl area. In summer 2022 ( Figures 12E, H ), there is a negative vorticity input in the eastern region, thereby weakening the south-western intrusion of the coastal flow to the west of Guangdong. The non-linear superposition of various conditions contributed to the generation of eastward transport during summer.