## 2 Data and Methods
### Methods  
The process of extracting internal tides from satellite data involves two main steps: data preprocessing and plane-wave fitting [18]. The purpose of the preprocessing step is to remove potential outliers and barotropic signals from the original data and calculate the parameters that will be used in the next step. Initially, all values along one track greater than 3.5-times the standard deviation and those measured at water depths less than 50 m were tagged as outliers and thus removed from the dataset. To extract the barotclinic signals from the original satellite data, the barotropic component predicted by the TOPEX/POSEIDON global tidal model (TPXO9) was subtracted. The wavelengths of the \(M_{2}\) barotropic and baroclinic tides are about 1000 and 100 km, respectively. Therefore, they can be distinguished using a spatial bandpass filter designed as a fourth-order Butterworth filter, considering their wavelengths. The thresholds were selected based on the wavelengths of the mode-1 \(M_{2}\) internal tides in the Arabian Sea as 50 km and 200 km. The wavelengths were computed by Equation (1) [23]:  
\[\omega^{2}=k^{2}c^{2}+f^{2}. \tag{1}\]  
Figure 1: (**a**) Spatial distribution of satellites’ sampling. The position of the selected window is shown on the top right corner, indicated by the red box. Blue dots show the satellite observation points, and the “line segments” actually are built up by dense dots. The red point is the central point of this window. (**b**) Topography of the Arabian Sea. The main locations involved in this study are indicated. This equation describes the dispersion relation of waves under rotation, where \(\omega\) represents the tidal frequency, \(k\) represents the wavelength, \(c\) represents the phase speed, and \(f\) represents the Coriolis frequency (\(f\equiv 2\Omega\sin(\varphi)\), with \(\varphi\) denoting the latitude). Phase speed \(c\) was calculated using the World Ocean Atlas (WOA) 2018 data and the following relation [23; 36]:  
\[\frac{d^{2}\Phi(z)}{dz^{2}}+\frac{N^{2}(z)}{c^{2}}\Phi(z)=0, \tag{2}\]  
where the vertical structures of displacement and vertical velocity are described by \(\Phi(z)\), which can be calculated by the eigenvalue of Equation (2), subject to the boundary conditions \(\Phi(0)=\Phi(-H)=0\) (H expresses the water depth), and \(N(z)\) expresses the buoyancy frequency profile. The Brunt-Vaisala frequency is defined as \(N^{2}=-(g/\rho_{0})d\rho/dz\), where \(\rho\) represents density and is calculated from the WOA dataset. Following these preprocessing steps, the data were used in the below plane-wave-fitting analysis. This research used the plane-wave-fit method proposed by [PERSON] and [PERSON] [23] to analyze internal tides. The main idea of this method is to separate the waves propagating to different directions at a given location (e.g., the red point in Figure 1a). After estimating \(k\) using Equation (1), the plane-wave-fit method [18] was applied to extract the features of the internal tides at the central point of one fitting window (e.g., Figure 1a), including their amplitudes, phases, and propagation directions. This method is similar to harmonic analysis, and both are based on the least-squares approach. Compared with harmonic analysis, which only uses a series of data at a given location, all data points within a \(1.6^{\circ}\times 1.6^{\circ}\) window were included in the plane-wave-fit calculation:  
\[SSH(x,y,t)=\sum_{i=1}^{3}a_{i}cos(\omega t+\varphi_{i}-kxcos\theta_{i}-kysin \theta_{i}), \tag{3}\]  
where \(SSH(x,y,t)\) represents the contribution of internal tides to the sea surface height (SSH). Here, \(\omega\) and \(k\) represent the frequency and wavenumber of \(M_{2}\) internal tides, respectively, and \(k\) can be obtained from Equation (1). Further, \(a_{i}\) and \(\varphi_{i}\) represent the amplitude and phase, respectively, which are the unknown variables in Equation (3), and \(x\) and \(y\) represent the east and north coordinates, respectively. Time is represented by \(t\).