We use the CROCO model (version 1.0, https://www.croco-ocean.org), which has proved useful for simulating internal tides (Guo et al., 2020a; Guo et al., 2020b). It is a split-explicit, free-surface ocean model based on the Boussinesq approximation. In our application, the wavelength of internal tides is much larger than the water depth (Vitousek and Fringer, 2011), and so following previous studies, we use the hydrostatic configuration (Jan et al., 2007; Jan et al., 2008; Rainville et al., 2010; Osborne et al., 2011; Powell et al., 2012; Kerry et al., 2013; Liu et al., 2019).

Our model has 40 vertical non-uniform sigma layers and a horizontal resolution of 1/30°. In our simulation domain ( Figure 1 ), topography is complex, with multiple islands, ridges, continental slopes and deep-sea basins. The topography data is from GEBCO2019 (general bathymetric chart of the oceans, https://www.gebco.net/), then further smoothed to avoid the pressure-gradient errors arising from sigma coordinates (Niwa and Hibiya, 2004). The minimum and maximum depths are 20 and 6000 m. To represent vertical mixing, we use the K-profile parameterization (KPP) mixing scheme (Large et al., 1994). For horizontal mixing, we use the Laplacian method. The explicit lateral viscosity is set to be 10 m² s^-1 and the bilaplacian background diffusivity is 20 m² s^-1 in the implicit diffusion of advection scheme (Cambon et al., 2018).

Table 1 summarizes the setup of numerical experiments. All of the experiments exclude both wind forcing and heat flux. The model initial and boundary conditions for temperature, salinity, velocity and sea level, for all of the experiments are from the HYCOM 10-year (2002-2011) monthly-averaged reanalysis dataset at a horizontal resolution of 1/12° (GOFS 3.0, https://www.hycom.org/dataserver/gofs-3pt0/reanalysis). For the winter cases in Table 1 , the 10-year monthly-averaged HYCOM outputs in January are used as the initial and boundary conditions; whereas for the summer cases, the July HYCOM outputs are used. In the control run (CtrlS and CtrlW in Table 1 ), tidal forcing is not specified at the open boundaries and thus the internal tides do not exist. For other experiments (cases with or without the Kuroshio in Table 1 ), they are driven by the barotropic semidiurnal tide (M₂) or the barotropic diurnal tide (K₁) at the open boundaries. The use of a single tide component here allows us to avoid the complex nonlinear interactions between internal tides of various components (Xie et al., 2010; Liu et al., 2015; Cao et al., 2018; Guo et al., 2020a; Song and Chen, 2020). The water levels and barotropic currents of each tidal component for the open boundary are obtained from the Oregon State University ocean topography experiment TOPEX/Poseidon Global Inverse Solution (TPXO9; Egbert and Erofeeva, 2002) tidal model with a horizontal resolution of 1/30°. The Flather condition (Flather and RA, 1976) is used for barotropic currents. Specifically, a 10-cell wide sponge layer is set at each lateral boundary to absorb the baroclinic energy and avoid energy reflection.

For the cases with the Kuroshio (M₂KS, M₂KW, K₁KS, K₁KW in Table 1 ), our model setup is the same as the control run, except that the tidal forcing is now included. As the initial and boundary conditions for temperature, salinity, velocity and sea level from the HYCOM outputs are realistic, the initial density and stratification are horizontally non-uniform, and the initial velocity field contains the Kuroshio information ( Figure S1 in the supplementary material ). In contrast, for the cases without the Kuroshio (M₂nKS, M₂nKW, K₁nKS, K₁nKW in Table 1 ), the initial and boundary conditions for temperature and salinity is obtained by spatial averaging the HYCOM outputs over the LS region. Therefore, both the density and stratification are horizontally uniform ( Figure S1 in the supplementary material ). Such uniform stratification has been used in previous studies on internal tides (e.g., Kerry et al., 2013; Xu et al., 2016; Wang et al., 2018; Chang et al., 2019). It is noted that the inflow/outflow boundary condition is specified in the with-Kuroshio cases and excluded in the without-Kuroshio cases. Based on the thermal-wind relation, a horizontally uniform hydrography corresponds to an absence of vertical shear, i.e, an absence of baroclinic flow. In addition, in these model runs (M₂nKS, M₂nKW, K₁nKS, K₁nKW in Table 1 ), the initial velocity is set to zero, so the Kuroshio flow is absent and thus its effects on internal tides are effectively removed.

Each case in Table 1 is run for 15 days for the following reasons. First, a 15-day period is long enough to allow internal tides to propagate through the model domain and this integration time was also used in previous studies involving the SCS or the LS (Jan et al., 2008; Xu et al., 2016; Guo et al., 2020b). Second, the domain-averaged sea surface height (SSH) and kinetic energy time series indicate that our model solution has reached equilibrium in 15 days (not shown). An even longer time of model integration may lead to large deviations between the simulation and the initial states, in particular, for the experiments without wind forcing and data assimilation. A relatively short model run can help our model solution maintain a climatological steady state that is close to that in the HYCOM reanalysis. Such a realistic Kuroshio is important for evaluating the Kuroshio’s effect on the internal tides. Third, the Kuroshio intrusion into the SCS has different paths ( Figure 1 ) and generally shows a seasonal pattern (Nan et al., 2011; Nan et al., 2015). However, the timescale of the dominant Kuroshio variability is O(100) days, and thus much longer than the timescale of internal tides (Zhang et al., 2001; Johns et al., 2001; Yin et al., 2017; Chang et al., 2019). Therefore, for the with-Kuroshio cases, we assume that the Kuroshio is quasi-steady and include the Kuroshio information in only the initial and boundary conditions. Using the quasi-steady setup can help us focusing on fundamental dynamics, making it easier to interpret the Kuroshio’s influence.

Our model setup has also been found useful in relevant previous studies. For example, to investigate the interaction between the M₂ internal tides and subtidal flow in Hawaii Ridge and the northeast of Taiwan, Zaron and Egbert (2014) and Chang et al. (2019) both used Simple Ocean Data Assimilation (SODA) ocean analysis data as the initial background field, and the models run for 14 and 30 days, respectively. Hsin et al. (2012) assessed the relative importance of open sea inflow/outflow, wind stress and surface heat flux in regulating Luzon Strait transport (LST) and its seasonality through several elimination model experiments. Their inflow/outflow at the open boundaries of the experiment NO was derived from the 10-year mean from North Pacific Ocean (NPO) model outputs, which is similar to our model settings of the with Kuroshio cases.

The change of Kuroshio invasion is related to the wind in the vicinity of the LS, such as the wind stress curl off southwest Taiwan (Wu and Hsin, 2012) and the East Asian monsoon (Hsin et al., 2012).But following Hsin et al. (2012); Jan et al. (2012) and Chang et al. (2019), wind forcing is not included in our experiments. The reasons are as follows. Winds can generate internal waves (Kitade and Matsuyama, 2000). Wind forcing can also influence the distribution and dissipation of internal tide energy (Hall and Davies, 2007; Masunaga et al., 2019). We avoid such added complexity here. Our current model setup without wind reduces the experiment to the simplest system, which allows us to focus on basic physics of the Kuroshio effect on internal tides. For the experiment NO in Hsin et al. (2012), the model runs are more stable in the absence of wind stress and heat flux forcing while fixing the open ocean inflows/outflows at the 10 year mean from NPO model.

To assess the effects of background stratification, and thus the Kuroshio flow shear, on internal tides, we use the energy equation of the baroclinic tides from Song and Chen (2020). In this equation, the variable a is decomposed into a background mean state a^m and a perturbation component a':

For example, a can represent pressure p, density ρ, or velocity vector u . The perturbation component can be further decomposed into barotropic a^bt and baroclinic modes a^bc ,

The variables η and H represent the ocean surface elevation and water depth, respectively. With this approach, the background field is described by the background mean-state variables (a^m ) and assumed constant.

Using the above variable decomposition, one can obtain the time-and-depth-integrated nonlinear baroclinic tide energy equations as

where Tran^bc denotes the transfer of energy from the mean flow to the baroclinc tidal flow:

Here ρ_c is the constant density and u_h means horizontal velocity vector. The term Tran^bc arises from the nonlinear interaction between the background shear and internal tides, and thus it is hereafter termed I_m-bc (Song and Chen, 2020). The next term Conv represents the total conversion rate of energy from barotropic to baroclinic tidal flow, hereafter called ‘total conversion rate’,

where I_bt-bc arises from the nonlinearity of the system and is zero in the linear case (Song and Chen, 2020). The term Conv_linear represents the linear conversion rate,

Conv_linear has three factors: the amplitude of the pressure perturbation at the bottom p'_θA (-H), the vertical component of the barotropic flow w_btθA (-H), and the Greenwich phase difference cos(θp’-θ_Wbt )), where subscript A indicates amplitude. For the derivation of Eq. (6), see Zilberman et al. (2011), Kerry et al. (2013) and Kerry et al. (2014). The term Div^bc from Eq. (3) quantifies the divergence of the energy flux:

where F_non ^bc represents the horizontal advection of local mechanical energy of internal tides (ke^bc + ape) by the background circulation. Finally, ε^bc the right side of Eq. (3), denotes the dissipation rate of internal tides. Compared to the energy equations in the linear system, the nonlinear energy equations have three additional terms: I_m-bc [Eq. (4)], I_bt-bc [Eq. (5)], and the horizontal divergence of F_non ^bc which is part of Div^bc [Eq. (7)]. For details of the derivation, see Song and Chen (2020).

To diagnose the above energy terms, the horizontal tidal velocity u'_h , the density perturbation ρ', and the pressure perturbation ρ' are first filtered with a specific band to extract the semidiurnal (1.83-2.20 cpd) and diurnal (0.80-1.15 cpd) motions. Given the relatively low frequencies of the background flow, we use a low-pass filter with a cutoff frequency of 0.3 cpd to extract subtidal motion. As Eqs. (3) - (7) are based on the assumption that the background flow is steady, we assume that the simulated background state is constant over three days, so we use the last three days of model outputs for the energy diagnosis.

In the background state, the simulated Kuroshio flows northward along the eastern side of the Philippines, crossing the Balintang Strait and then entering the SCS ( Figure 2 , top row). It is about 120 km wide, with maximum speeds of ~1.2 m s^-1 in winter and ~1.3 m s^-1 in summer, averaged over the upper 100 m of the 121°E meridional section ( Figure 2 , the bottom row). Part of the Kuroshio flows out along the southern part of Taiwan and then travels northward along the eastern part of Taiwan. Overall, the modeled flow agrees well with that from the SODA reanalysis data. The RMSEs (root mean square error) of the current speed for cases CtrlS and CtrlW are 0.15 m s^-1 and 0.16 m s^-1 in the whole domain, respectively. The probability that the modeled current direction differs from SODA by less than 60° is ~70%, with less than 35° being about 50%, indicating good agreement.

Concerning SSH, we find the results in Figures 2C, D to closely agree with the altimeter data in Nan et al. (2011) (not shown in the figure). Moreover, the difference between the maximum and minimum SSH values, for both the simulation and SODA reanalysis data, exceed 0.50 m, which is much greater than the RMSE values (0.08-0.09 m). Also, the correlation coefficients are greater than 0.84. Thus, the SSH also agrees well with SODA.

Concerning the three main Kuroshio paths into the SCS, the leaking path dominates in winter, the leaping path dominates in summer, and the looping path instead appears just southwest of Taiwan more frequently in winter than in other seasons (Wu and Chiang, 2007; Nan et al., 2015). Similarly, our simulation generally has the leaking path in summer and the looping path in winter ( Figures 2C, D ). Overall, the simulated Kuroshio velocities and paths appear reasonable.

Our model also reasonably captures the background hydrographic fields. For example, the meridional cross-sections of temperature and salinity at 120.75°E agree well with those from WOA18 ( Figure 3 ). The RMSEs of temperature for cases CtrlS and CtrlW are both less than 1, and the RMSEs of salt are both less than 0.2, with the correlation coefficients of temperature exceeding 99%. In addition, for both WOA18 and our model solution, the isotherms in the LS region tilt toward the east (the Philippine Sea) (not shown).

Using the harmonic analysis method from Pawlowicz et al. (2002), we obtained the cotidal charts of M₂ and K₁ tides from the last three days of modelling outputs. In Figure 4 , we compare the simulated results to the TPXO9 dataset. The distribution of simulated surface elevation and cophase lines are generally consistent with TPXO9, except that the simulated amplitude is slightly larger than TPXO9 near the shelf region. Small-scale fluctuations of the cophase lines occur in the LS and SCS regions, and these are mainly caused by the baroclinic tide signals (Niwa and Hibiya, 2004; Xu et al., 2016). Consistent with the altimetric estimates from Zhao (2014), the horizontal wavelength of M₂ and K₁ internal tides in our simulation, inferred from these fluctuations, are about 160 and 320 km, respectively. The RMSEs of the amplitudes are 11.70 and 6.95 cm, respectively for cases M₂KS and K₁KS, and the correlation coefficients of the amplitudes exceed 92% with a confidence level of 95%. Here, we further calculate a new RMSE specifically for tides in another way. Following Cummins and Oey (1997),

where A and ϕ represent the amplitudes and phases of the surface elevations. The subscripts o and m refer the results from TPXO9 and our model, respectively. The new mean RMSEs (calculated using Eq. 8) of the eight cases (M₂KS, M₂KW, M₂nKS, M₂nKW, K₁KS, K₁KW, K₁nKS and K₁nKW) are 4.50, 5.03, 4.58, 5.39, 4.34, 4.35, 4.37, and 4.38 cm, respectively.

To assess our model fidelity on internal tides, we ran additional experiments that include semidiurnal tides (M₂ and S₂) and diurnal tides (K₁ and O₁) at the open boundaries. Results in Figures 4E, F show that the simulations well capture the propagation pattern of internal tides. For both diurnal and semidiurnal tides, the simulated energy-flux directions agree with observations from Alford et al. (2011). The energy-flux magnitudes also agree with observations, particularly in the region south of 20°N. We provide further statistics in Table 2 . In the semidiurnal case, the simulated maximum energy flux is 55.2 kW m^-1 in the LS region (120°E-122.5°E, 19°N-21.5°N), which is much greater than the energy flux RMSE of 9.1 kW m^-1, and the average energy flux is 15.3 kW m^-1 ( Table 2 ). The RMSE of the energy-flux direction is no more than 23.5°. In the diurnal case, the RMSE of energy flux is 7.1 kW m^-1, and that of the energy-flux direction is 32.5°. By including all semidiurnal tidal components in our simulation, smaller deviations would likely be obtained.

In the simulations with the Kuroshio, the internal tide energy flux and conversion rate are spatially heterogeneous. Figure 5 shows that the dominant generation site of internal tides is in the LS region. Here, the M₂ tide generates energy of ~17.1 GW (case M₂KS), whereas K₁ generates ~10.6 GW (case K₁KS). As marked by the light blue arrows in Figure 5A , there are three M₂ energy beams radiating from the LS: one points to the northeastern shelf of the SCS (northern beam), another directs southwest to the deep-sea basin (southern beam), and the third one east to the Pacific Ocean. This three-beam pattern is generally consistent with previous observational and modelling studies (Kerry et al., 2013; Zhao, 2014; Xu et al., 2016). The northern beam heads northwestward, converging near the Dongsha Islands, and then dissipating at the continental slope. The southern beam, much weaker than the northern beam, propagates southwestward, crossing the deep-sea basin, and then traveling for more than 1000 km to reach the Vietnam coast or the Nansha Islands due to its small dissipation rate. For the K₁ internal tides, two main beams radiate from the LS: one points southwest to the deep-sea basin of the SCS and the other points east to the Pacific ( Figure 5B ). The complex wave forms and propagation paths of internal tides radiating from the LS into SCS may be due to the impact of rotation, complex topography, and the mutual interference with locally generated internal tides.

Besides the LS region, local generation sites of internal tides occur within the SCS. The strength of generation of these internal tides is controlled by both the topographic slope and the strength of barotropic tides. Also, the strength distribution of the K₁ barotropic tides differs from that of M₂ ( Figure 4 ). Consequently, the generation sites of M₂ internal tides differ from those of K₁ within the SCS. In particular, the generation of M₂ lies mainly at the northeastern continental slope (box “1” in Figure 5A ) and the southwestern continental slope (box “2”), whereas those for K₁ are similarly marked in Figure 5B , being around the Xisha Islands and the northwestern continental slope (“1”), and around the Nansha Islands and the southwestern continental slope (“2”). As shown in Zu et al. (2008), this difference is due to the fact that the energy flux intensities of the two barotropic tidal components vary in different regions within the SCS.

The seasonality of internal tides can be evaluated by comparing the summer cases in the top row of Figure 5 to the winter cases in the bottom row. In this realistic stratification scenario, the source and propagation direction of the tides in the SCS vary little from season to season. However, as shown in the next section, the energy magnitude undergoes a significant seasonal cycle.

In the nonlinear case, there is energy exchange between the background shear and internal tides. The corresponding energy exchange rate I_m-bc can be quantified using Eq. (4). Overall, the magnitude of I_m-bc is only about 1-10% of the magnitude of the total conversion rate (Conv). Considering the spatial pattern of I_m-bc in summer, both its peak magnitude and extent are larger in the with-Kuroshio case than those in the without-Kuroshio case ( Figure 7 ). The tidal current can be generated by tidal forcing from the open boundary, which produces relatively weak background flow in the absence of Kuroshio. This weak background flow contributes a relatively small energy exchange rate between the background state and the internal tides ( Figure S2 in the supplementary material ).

Consider the horizontal shear Shear_H and vertical shear Shear_v :

The distributions of Shear_H and Shear_V for cases M₂KS and M₂nKS are shown in Figure 8 . In the with-Kuroshio case, the Shear_H at 100 m depth along the Kuroshio path typically exceeds 1 × 10^-5 s^-1 ( Figure 8A ). In the without-Kuroshio case, the Shear_H is generally less than 1 × 10^-5 s^-1 over the whole region ( Figure 8B ). For the Shear_v , the large values are also mainly distributed along the path of the Kuroshio, and the magnitude for the without-Kuroshio case is much smaller over the whole region ( Figures 8C, D ). In general, the notable background shear is mainly distributed along the Kuroshio path.

The energy exchange I_m-bc in Figure 7 is large in the regions where both the background shear ( Figure 8 ) and internal tide energy ( Figure 5 ) are large. For the M₂KS case ( Figure 7A ), the magnitude of I_m-bc in the LS region is about 10^-2 - 10⁰ W m^-2, which is 1-2 orders of magnitude larger than that in the open ocean (about 10^-2 W m^-2). Along the Kuroshio path (e.g., Hengchun Ridge and northeastern Taiwan), both the background shear and internal tide generation are intense, corresponding to a large I_m-bc . For case K₁KS ( Figure 7C ), elevated values of I_m-bc mainly occur in the LS region. As the continental slope near northeast Taiwan is not the main generation region of the K₁ internal tides, the magnitude of I_m-bc is small.

Besides the nonlinear interaction between the background shear and internal tides, there is also nonlinear interaction between barotropic and baroclinic tides (Section 2.2). The term I_bt-bc from Eq. (5) represents the energy exchange rate between barotropic and baroclinic tides. Compared to I_m-bc , large values of I_bt-bc cover a smaller area ( Figures 9A-C ). For case M₂KS, I_bt-bc is large in areas with prominent topography features, namely the East China Sea continental slope northeast of Taiwan and the two ridges of the LS (the Dongsha Plateau and the Ryukyu Ridge) ( Figure 9A ). In case M₂KS, I_bt-bc is also large along the Kuroshio path, especially in the southern part of the Hengchun and Lanyu ridges. For case K₁KS, however, the area with elevated values of I_bt-bc is smaller than that of case M₂KS, mainly along the two ridges of the LS ( Figure 9C ).

For case M₂KS, the magnitude of I_bt-bc is only about 10^-3-10^-1 W m^-2 in the LS region, much smaller than the total conversion rate, and even much smaller than I_m-bc . Therefore, I_bt-bc is a negligible nonlinear energy term in this case. Nevertheless, it can shed light on the source of internal tides and the Kuroshio effect. For both M₂ and K₁ internal tides, the magnitudes of I_bt-bc in the without-Kuroshio case are much smaller than those in the with-Kuroshio case ( Figure 9 ), indicating that the Kuroshio tends to promote the interaction between barotropic and baroclinic tides. This sensitivity of I_bt-bc to the Kuroshio may be related to baroclinic tide generation. In particular, baroclinic internal tides can be generated when the background flow with temporal variability at tidal frequency passes a topography obstacle (e.g., a seamount or ridge). In the without-Kuroshio case, weak background flow shear probably corresponds to weak flow on topography. The weak flow decreases the energy exchange rate between barotropic and baroclinic tides and weakens the baroclinic internal tide.

In the baroclinic tide energy equation, the divergence of energy flux is induced by two processes: the advection of tide energy by background flow and the work done by pressure [Eq. (7)]. The contribution from advection is represented by the nonlinear term F_non ^bc in Eq. (7). Representative results are shown in Figure 10 . For the with-Kuroshio cases in summer (M₂KS and K₁KS, Figures 10A, B ), the F_non ^bc flux is largest in the LS, where it points northwestward and crosses the Hengchun Ridge to the SCS. Then, in winter (M₂KW), the F_non ^bc flux is slightly stronger than in summer, directing northwestward and converging in the Dongsha Plateau ( Figure 10C ). The flux pattern here is similar to the leaking path of the Kuroshio ( Figure 1 ), indicating that the contribution of the Kuroshio to F_non ^bc is important for the with-Kuroshio case. For the M₂nKS case ( Figure 10D ), which does not include the Kuroshio, the maximum energy flux is only ~3 kW m^-1. It is much smaller than the maximum value of ~15 kW m^-1 for the with-Kuroshio cases shown in Figures 10A-C . The direction of the F_non ^bc flux seems random, with little resemblance to the Kuroshio path. The small value of energy flux in the without-Kuroshio case arises from weak background flow near the LS.