*2.2. ISW Underwater Structure Reconstruction*

ISW amplitude *η*<sup>0</sup> can be derived from satellite images using the following equation [28]:

$$\eta\_0 = \frac{12\beta}{a l^2} = 1.32^2 \cdot \frac{12\beta}{a D^2} \,\mathrm{}^{\prime}\tag{1}$$

where *D* is the distance between the center of light and that of dark stripes, and *l* is the half-wavelength of the ISW. *α* and *β* are nonlinear coefficient and dispersion coefficient, in the Korteweg–de Vries (KdV) equation, and they are calculated from *α* = 3 R 0 −*H c*0 h *∂φ*(*z*) *∂z* i3 *dz*/2 R 0 −*H* h *∂φ*(*z*) *∂z* i2 *dz* and *β* = *c*<sup>0</sup> R 0 −*H* [*φ*(*z*)]<sup>2</sup> *dz*/2 R 0 −*H* h *∂φ*(*z*) *∂z* i2 *dz*, respectively. Here, *φ*(*z*) is the first-mode vertical eigenfunction of the wave and governed by Sturm–Liouville equation:

$$\frac{d^2\phi}{dz^2} + \frac{N^2(z)}{\mathbb{C}\_0 r^2}\phi = 0,\tag{2}$$

with the boundary conditions *φ*(0) = *φ*(−*H*) = 0, where *C*<sup>0</sup> is the eigenspeed of the Equation (2) and *N*<sup>2</sup> (*z*) represents the background stratification obtained from monthly climatological density profiles from the WOA18 (World Ocean Atlas 2018) dataset. This method has been proved to be reliable by comparing the inversion results of ISW amplitude from MODIS images in the northern South China Sea with the wave amplitude measurements from mooring observations [26].

**Figure 1.** Optical remote sensing images over the BS acquired on 12, 14, and 19 April 2021 (local time) by MODIS.

**Figure 2.** Optical remote sensing images over the BS acquired on 15, 16, 18 and 21 April 2021 (local time) by VIIRS (**a**–**c**) and by MODIS (**d**).

The waveshape of the ISW with amplitude *η*<sup>0</sup> can be obtained based on the solution to the KdV equation:

$$
\eta(\mathbf{x}, z) = \eta\_0 \mathrm{sech}^2 \mathbf{x} \cdot \phi(z). \tag{3}
$$

As an ISW arrives, the isopycnal surface, initially at depth *z*, is depressed to the depth *z* 0 = *z* + *η*. After rotating the horizontal axis of the coordinates to the ISW propagation direction, the continuity equation is written as:

$$\frac{\partial u(\mathbf{x}, z')}{\partial \mathbf{x}} + \frac{\partial w(\mathbf{x}, z')}{\partial z'} = \mathbf{0},\tag{4}$$

where *u* and *w* are the horizontal current along the wave propagation direction and the vertical current, respectively. The vertical velocity is regarded as the partial derivative of isopycnal displacement with respect to time, and thus *w*(*x*, *z* 0 ) = *∂η*(*x*, *z*)/*∂t* and the first derivative of isopycnal depth in the vertical direction is given as *∂z* 0/*∂z* = 1 + *∂η*/*∂z*. Therefore, local along-isopycnal horizontal current can be calculated from:

$$u(\mathbf{x}, \mathbf{z}') = -c \frac{\partial \eta(\mathbf{x}, \mathbf{z})}{\partial \mathbf{z}'} = -\frac{c \frac{\partial \eta(\mathbf{x}, \mathbf{z})}{\partial \mathbf{z}}}{1 + \frac{\partial \eta(\mathbf{x}, \mathbf{z})}{\partial \mathbf{z}}}.\tag{5}$$
