A Characteristics Set Computation Model for Internal Wavenumber Spectra and Its Validation with MODIS Retrieved Parameters in the Sulu Sea and Celebes Sea

: The quasi-linear or nonlinear interactions among different ocean motions dominate the system internal structure and appearance feature presented in spatial and temporal evolution. However, deficiency of the characteristics set computation model for internal wavenumber spectra proves to be a serious barrier to derive interaction mechanisms of internal waves with large or small scale ocean motions. In this study, a characteristics set computation model for internal wavenumber spectra is proposed for complicated offshore environments. The refraction of current shear instability, bottom topography and the reflection at surface and bottom are attentively considered in the complicated characteristics inlaid scheme. Model results are validated with MODIS retrieved internal wave parameters in the Sulu Sea and Celebes Sea. This original characteristics set computation model for internal wavenumber spectra can be used widely and can further improve the understandings of generation, dissipation, nonlinear wave-wave interaction and mixing process of internal


Introduction
The quasi-linear or nonlinear interactions in the ocean dynamic system dominate the physics of internal gravity waves [1][2][3]. The energy spectrum method is a practical and efficient way to interpret the internal wave composition structure and appearance properties. However, deficiency of the characteristics set computation model for internal wavenumber spectra proves to be a serious barrier to analyze interaction mechanisms of internal waves with large or small scale ocean motions. In fact, the characteristics/ray theory has been studied and successfully applied to the spectral propagation of surface waves [4][5][6][7]. It is probable to deeply recognize the horizontal anisotropy, high nonlinearity, spectrum denaturation of internal waves in the offshore if we follow the quite similar treatment of surface waves through the internal wave characteristics set modeling, which can interpret further the propagation physics of internal waves compared to the traditional models based on the primitive Navier-Stokes or KdV equations. Furthermore, satellite imagery can detect packets of internal waves which are common features propagating in the seas [1,8]. These satellite retrieved internal wave parameters (e.g., wavelength, phase velocity, etc.) are useful particularly and expediently to validate the field distributions obtained from the wave characteristics set modeling.
The internal waves appear in the satellite imagery as alternating bands of light and dark strips that result from sea surface roughness variations, which are due to the creation of convergent (rough) and divergent (smooth) zones set up by the internal wave currents that move across the surface in phase with wave's subsurface crests and troughs [9,10]. MODIS (Moderate Resolution Imaging Spectroradiometer), SAR (Synthetic Aperture Radar) and other sensors on board satellites have the capabilities for large area detection of internal waves in their wide swath modes. The internal wave parameters can be obtained directly by using pairs of satellite images separated in time by only a few minutes to a few hours, or retrieved through the TMI (Tidal Period Images) method for application to the internal wave field studies [11]. In this study, we will apply the MODIS retrieved internal wave parameters to verify the proposed characteristics set computation model. This paper is organized as follows. A derivation of the characteristics set computation model is given in Section 2.1. Ocean data related to the MODIS sensors and varying environmental modeled data from a wave-tide-circulation coupled model are described in Section 2.2. The applied results of the characteristics set computation model are presented in Section 3. Finally, a discussion and conclusions are presented in Sections 4 and 5.

Model Derivation and Set-Up
denote the kinetic and potential internal wave energy and SM  denotes the Reynolds average in wave motion.
Otherwise the notation is standard. The first term on the left-hand side of Equation (1) is related to the local mechanical energy variation and the second and third ones denote the energy flux transferred by internal waves and background currents. The first and second terms on the right-hand side of Equation (1) are related to the modulation by larger scale motions through shear instability generations. The third term is related to the energy input through thermal radiation, the fourth and fifth ones are related to the modulation by smaller scale motions through ocean mixing and the last two terms are related to the energy loss rate due to internal viscosity.
, , E k E k k k   designate an internal wavenumber energy spectrum for unit mass water, i.e., and the corresponding wave energy balance equation can be written as: Correspondingly, There are also other boundary input source functions, such as the wind and surface gravity wave energy input terms not included here. In this study, we are concerned entirely with the propagation properties of wave-number spectrum, so we let further numerical modeling performed below.

3-Dimensional Complicated Characteristics Set Equations
In the preceding derivations, we implicitly assume the horizontal invariance of density, topography, etc., as in Yuan et al. [2]; the complex frequency-wavenumber relation at the sea surface can be simplified as: where f denotes the Coriolis parameter; ˆN H represents the depth where (0)= ( ) , Equation (6) can also be practically expressed as: where 1, 2,... i  denote the vertical mode numbers of internal gravity waves. The vertical wavenumber 3 k satisfies: So the group velocities can be obtained as: The ray method or characteristics method is valid for a study of both linear and nonlinear wave packets propagating both in homogeneous and inhomogeneous media [4,5,16,17]. Here we present the practical characteristics equations for further interpreting the refraction induced by topography and current. The characteristics equation describes the propagation law of ocean waves in physical space and it can be written as: By using the motion equations of waves, with î i k U      , after some manipulation, the variation laws of the modulus and azimuth angles of wavenumber can be derived as: Differentiated by Ĥ in both sides of Equation (9), we obtain: where (13) and (15)-(17)

Reflection at the Surface and Bottom
The characteristics or ray theory predicts that internal waves reflect at surface or at the reflection level ˆN H which is supposed as ˆN H H  below. The boundary condition of equations of internal wave motion at bottom is [2]: ( , , ), ( , , ) k k k k l l l l     and the reflection coefficient A satisfy the conservation laws of kinematics [18] and we obtain: The reflected wave satisfies 0 up C n g     at bottom. By employing Equation (12), it yields: So that in numerical modeling, Formula (25) can be seemed as a criterion if the discrete wavenumber is a reflected one or not. Moreover, at surface, Formula (25) is simplified as 3 0 l  . Consequently, the discrete reflected wavenumber spectrum roughly satisfies: where the incident wavenumber may be not the discrete one and its spectrum   1 2 3 , , E k k k is 3-dimensionally and linearly interpolated among the adjacent eight discrete wavenumbers in the phase space.

Grid Distribution and Complicated Characteristics Inlaid Scheme in Physical and Phase Spaces
In physical space, we take the horizontal grid parallel to the longitude and latitude lines in the study domain with the resolution of 1°/3 by 1°/3 or 1°/6 by 1°/6 and the uniform vertical grid spacing of 20 m. The maximum water depth is set to 3000 m. Figure 1 is our model computational domain and bathymetry.
In wavenumber space, we take polar coordinate grids according to the proposed negative exponential form of empirical internal wavenumber spectrum E(k)~k , i.e., k = k e ( )•∆ , i = 1,2, … , K + 1, where ∆k = log ( ) The horizontal and vertical angles of wave direction are discretized uniformly below: The above model settings can be modified according to different requirements. The complicated characteristics inlaid scheme is designed as Figure 2.  Wave statistics based on the computational spectrum have a spatial resolution between 250 m and 1 km, dependent on the collection wavelength [19] and with a geolocation accuracy of at least 60 m (1σ) [20]. Full details of the characteristics of MODIS and its associated data products are given on the NASA MODIS website (modis.gsfc.nasa.gov, accessed on 30 November 2021). Previous studies have been done on the inversion mechanisms of internal waves from satellite data [8][9][10][21][22][23]. Here in this study, the Aqua-MODIS images with a spatial resolution of 250 m are used to retrieve the internal wavelengths, here defined as the distance between two successive wave packets and the corresponding apparent phase velocities.

Modeled Data
The modeled data were presented by using the Key Laboratory of Marine Science and Numerical Modeling (MASNUM) wave-tide-circulation coupled model [24], which employs the MASNUM wavenumber spectral model [6,7] and the POM circulation model [25]. The key mixing role induced by surface waves in the formation of upper mixed layer [15,[26][27][28] was considered in the coupled model. The coupled model covers the computational domain (10° S-30°N, 90-135°E) with the horizontal resolution of 1°/30 by 1°/30 and the vertical 51 sigma layers. The model was integrated for 10 years and the gridded products of SSH anomaly, circulation pattern, subsurface temperature, etc., were validated to the satellite altimetry and in situ observations [24]. The simulated climatological monthly-mean temperature, salinity, etc., of year 10 are used in this study.

Simplified Wave Frequency Dispersion Relations
The vertical mean wavenumber can be obtained through a non-trivial iteration scheme according to Equation (9) and the following transformation of Equation (10): The vertical integration of Equations (9) and (30) constitute the simplified wave frequency dispersion relations. As the result of introducing the vertical mean group velocity, according to Equation (12), it can be expressed as follows: Obviously, the phase velocity   ago, with 5 groups of wave packets propagating across the Sulu Sea to the northwest and 4 groups propagating across the Celebes Sea to the southeast. The apparent phase velocity can be retrieved through the TMI (Tidal Period Images) method [11] based on the distance between two successive wave packets, which can be obtained directly from the image. The mean velocities are 2.26, 2.88 m/s in the Sulu Sea and Celebes Sea, respectively. The corresponding retrieved parameters to wave packets are listed in Table 1.

Celebes Sea
Sulu Sea length is homogenous in the Celebes Sea with the quantity of approximately 140 km, but varies much in the Sulu Sea from 70 km to 140 km. It is evident that this responds to the varied terrain, especially in the shallow water area. Figure 5 shows the modeled phase velocity ap H C of mode 2 in the Sulu Sea and Celebes Sea. Similarly, it varies from 1.8 m/s to 3.0 m/s in the northwestern Sulu Sea and reaches up to 3.4 m in the central and eastern region, which is a little larger than the MODIS retrieved ones listed in Table 1. However, in most of the Celebes Sea, it has a uniform velocity of approximately 3.0 m/s, which agrees well with the MODIS retrieved ones. The modeled wavelength and phase velocity of mode 2 are also listed in Table 1 near the corresponding locations of MODIS observations. Both are relatively consistent except for Packet 1 in the Sulu Sea, which is mainly due to the varied terrain in the offshore area with coarse resolution of our model.    Figure 6, which has similar properties as the phase velocity interpreted above. Its characteristic scale is about 2.0-3.5 m/s for semidiurnal internal tide in deep water area. Three-dimensional group velocity plays a key role in the propagation of internal wavenumber spectrum, which will be discussed below.

Analysis of Spectral Propagation
We wish to explore the internal wave spatial propagation behavior; the integral of internal wavenumber spectrum at any depth layer is presented for further conveniently interpreting, i.e., In order to examine the applicability of the characteristics set computation model stated above, the initial GM75 spectrum [29] is set radially in the Sulu Archipelago with two major directions, northwest to the Sulu Sea and southeast to the Celebes Sea. The spectral propagation can be considered and reduced to the following two successive stages: P1 and P2: During the first stage P1, the wavenumber spectrum mainly propagates vertically with the low horizontal group velocity near the Sulu Archipelago, but reflects at surface and bottom frequently, so the mean wave direction appears as much scattering (Figure 7). During the second stage P2, the wavenumber spectrum propagates thoroughly across the deep water area in the Sulu Sea and Celebes Sea; besides reflecting at surface and bottom, the mean wave direction concentrates in an orderly manner (Figure 8). Furthermore, the refraction of bottom topography considered in Equations (15)- (17) and (19) plays a perceptible role in the western and eastern offshore areas of the Sulu Sea and the western offshore area of the Celebes Sea. Figure 8b indicates that internal waves also propagate affluently to the northeast in the Sulu Sea, as has previously been found from satellite observations [30,31]. However, the western propagating waves detected by [8] are not simulated here in the Celebes Sea, which needs to be further studied in the future.

Discussion
Fine resolution, large swath MODIS or SAR imagery has the ability to survey the properties of internal wave field evolution and applicability to interpret and verify the inherent dynamic mechanisms of wavenumber spectral propagation. The complicated characteristics inlaid scheme is an efficient way to simulate the internal wave propagation. In this study, a characteristics set computation model for internal wavenumber spectra is proposed for complicated offshore environments. The refraction of current shear instability, bottom topography and the reflection at surface and bottom were tentatively considered in the computational scheme. A rather cautious derivation was employed for the interaction dynamics which dominates the wave internal and external properties. However, there are still other processes, such as the enhancement near the turning depth or latitude [32][33][34] which are not considered here, which will be studied later.
Model grid discretization and computational configuration were designed for numerical experiments in the Sulu Sea and Celebes Sea. The modeled horizontal wavelength ap H  and phase velocity ap H C for mode 2 of semidiurnal internal tide agree well with the MODIS retrieved ones. Spatial distributions indicate that these parameters vary much in the Sulu Sea but homogenously in the Celebes Sea, which is largely due to the former varied terrain in the northern offshore. In situ survey data are needed for further comparison with the modeled results.
Numerical spectral propagation experiments indicate that after the first distracting period near the Sulu Archipelago, the waves propagate concentrately in an orderly manner in the deep water area, in which the reflection at surface and bottom plays an important role in long-range propagating and the refraction of bottom topography plays a perceptible role in the offshore areas. The internal waves also propagate affluently to

Conclusions
Packets of internal waves propagating in the seas are common features detected by large amounts of satellite imagery. A characteristics set computation model for internal wavenumber spectra is proposed for complicated offshore environments to interpret further the propagation physics of internal waves. The refraction of current shear instability, bottom topography and the reflection at surface and bottom are attentively considered in the complicated characteristics inlaid scheme. Theoretical analysis and numerical modeling were implemented in the Sulu Sea and Celebes Sea and evaluated with the MODIS retrieved internal wave parameters. The results suggest that the original characteristics set computation model for internal wave spectra can be used widely and will be helpful in improving the understandings of generation, dissipation, nonlinear wave-wave interaction and mixing process of internal waves.

Data Availability Statement:
The data that support the findings of this study are not publicly available. However, data are available from the authors upon reasonable request.