Oceanic internal solitary waves at the Indonesian submarine wreckage site

On 21 April 2021 local time (20 April UTC), the Indonesian Navy submarine (KRI Nanggala-402) sank near the Lombok Strait, ~100 km north of the Bali Island (see magenta star in Fig. 1a), with 53 crew members dead. On the basis of Moderate Resolution Imaging Spectroradiometer (MODIS) satellite images (Jackson, 2007), NASA demonstrated that powerful “underwater waves” happened in the treacherous region and likely hit the vessel resulting in its disappearance (https://www.npr.org/2021/04/30/992496772/). Besides, NASA reported that the collapse depth of submarine KRI Nanggala-402 was ~200 m, but official reports on the local underwater waves and the voyage depth of the submarine were still lacking. This phenomenon was referred to as oceanic internal solitary waves (hereinafter ISWs). An important behaviour of ISWs is causing large vertical displacements and downward currents in a short time in the ocean interior, and therefore may drag the submarine down to the collapse depth (i.e., ~200 m), where the water pressure exceeds the endurance limit of the submarine.

Figure  1.  Bathymetry around the Indonesia including the Lombok Strait (a), the submarine wreckage location is marked as a magenta star; a satellite image on 25 April 2005 derived from the MODIS (b), in which blue arrows represent wave propagation directions; and numerically-predicted sea surface height gradients ($ \left|\nabla \eta \right| $) at 01:00 UTC on 20 April 2021 (c).

DownLoad: Full-Size Img PowerPoint

There are three conditions generating oceanic ISWs, namely the oscillatory surface (i.e., barotropic) tides, the abrupt topography (e.g., underwater sill and slope), and the stratified water. The Lombok Strait was well-known for all three (Murray and Arief, 1988; Mitnik et al., 2000; Ningsih et al., 2010; Purwandana et al., 2021), so whether the ISWs were the culprit causing the shipwreck in the Lombok Strait was of particular interest. To address this problem, we need to understand the physical dynamics and spatial characteristics of ISWs. Here two approaches were used to investigate ISWs in the Lombok Strait: first using the satellite image to present the spatial variability of ISWs, second numerically reproducing the ISW dynamics on the day of the accident.

Convergence/divergence of the currents induced by ISWs on the ocean surface contributed to strong modulations of sea surface roughness, thereby resulting in distinctive features in the optical true-color images (e.g., Susanto et al., 2005). A snapshot on 25 April 2005 (Fig. 1b) depicted that the Lombok Strait generated ISWs radiated both northward and southward. Three stages of northward-going ISWs were captured in the satellite image (Fig. 1b), namely the generation, propagation and shoaling processes. The ISW was presented as an isolated wave with a short crest length on the stage 1, but converted to a wave packet with a long crest length on the stage 2. Eventually on the stage 3, the wave packet approached the wreckage site and reached the shallow region near the Kangean Island (Fig. 1b).

Based on ~7 000 MODIS images over the past 20 years from 2002 to 2021, we identified wave occurrences in April and estimated wave amplitudes by a theoretical method (KdV theory, Ostrovsky and Stepanyants, 1989).

where the parameters $ {c}_{0} $, $ \alpha $ and $\, \beta $ are linear long-wave phase speed, nonlinear coefficient and dispersion coefficient, respectively. $ \alpha $ and $\, \beta $ are coefficients associated with background continuous stratification and water depth. Here we calculated spatio-temporal varying $ \alpha $ and $ \,\beta $ using monthly climatological density profiles from the WOA18 dataset (World Ocean Atlas 2018) with a horizontal grid resolution of 0.25°. Amplitudes of the ISWs can be extracted in the satellite images (Zheng et al., 2001), following the equation:

where $ D $ is the distance between the center of the bright and dark stripes within the leading wave in the MODIS true-color images, $ l $ is the characteristic half width of an ISW, and $ {\eta }_{0} $ is the estimated wave amplitude.

Regardless of northward- or southward-going ISWs, ISWs were commonly appearing in April (i.e., the month of submarine sinking event), which provided the evidences that ISWs could be the culprit to the shipwreck in April 2021. In terms of northward-going ISWs, the KdV theory was applied to extract wave amplitudes on different stages (Fig. 1a) based on the distance between the bright and dark stripes in the MODIS images. Overall, ISWs had the largest amplitude of ~70 m shortly after generating over the sill (Stage 1), then decreased from ~50 m in the deep water (Stage 2) to ~35 m at the wreckage site (Stage 3).

In comparison with satellite observations, numerical simulations are a more effective approach to characterize the ISW structures in the ocean interior and reproduce wave dynamics in the Lombok Strait (Aiki et al., 2011). Hence, we implemented a three-dimensional primitive equation ocean solver (MIT general circulation model, MITgcm) in the nonhydrostatic mode with realistic conditions during a spring tidal period from 00:00 UTC 17 April 2021 to 00:00 UTC 22 April 2021. Cable News Network (CNN) reported that the submarine started diving at 03:00 on 21 April local time (20:00 UTC 20 April) but lost contact at 04:25 on 21 April local time (21:25 UTC 20 April), so we mainly focused on 20 April (UTC) 2021 (Fig. 2). Model configurations were presented as follows.

Figure  2.  Numerically-predicted sea surface height gradients ($ \left|\nabla \eta \right| $) (a) and numerically-predicted isotherms ($ T $) from 28.5°C to 9°C (bending line in b, d and f from top to bottom) with a interval of 1.5°C and wave-induced velocity ($ {U}_{\rm {bc}} $, color shade) along the transect (red dashed line in panel a) (b) at 05:30 UTC on 20 April 2021; c and d are the same as a and b but at 13:30 UTC on 20 April 2021, and e and f are at 21:30 UTC on 20 April 2021. The wreckage site is marked as the magenta star and ISWs on three stages are marked in blue boxes. z represents water depth.

DownLoad: Full-Size Img PowerPoint

Model bathymetry data was derived from the ETOPO1 global dataset (Amante and Eakins, 2009). To ensure these waves were physically derived rather than products of numerical dispersion (Vitousek and Fringer, 2011), a resolution condition ($ \Delta x<{h}_{\rm p} $) needs to be satisfied, where $ {h}_{\rm p} $ was defined as the depth of the internal interface (pycnocline depth). According to the WOA18 climatology dataset of stratification, $ {h}_{\rm p} $ was approximately equal to 100 m in April, we therefore set horizontal cell ($ \Delta x $) as 100 m in both the longitudinal and latitudinal directions. The entire model domain consisted of 3 000$ \times $2 000 grid cells with 60 vertical layers ranging from 10 m from the surface to 100 m near the sea bed. The initial model stratification was obtained from the climatology dataset WOA18 by spatially averaging the monthly output in April, resulting in horizontally homogeneous temperature initial conditions.Salinity was set as a constant value of 34.5. The model was driven by eight main tidal constituents (M₂, S₂, N₂, K₂, K₁, O₁, P₁, Q₁) on the boundaries with values derived from the Oregon State University TOPEX/Poseidon Solution (TPXO8-atlas data) with (1/30)° resolution. A 10 km wide sponge layer was imposed on each lateral boundary to absorb internal waves and avoid reflection back to the inner region. We started a 5-d simulation from 17 April (UTC) 2021. Quasi-steady conditions occurred after 3 d, so the model results were analyzed over the remaining 2 d (20–21 April). More details of model configurations were shown in Table 1.

Table  1.  Parameters in the numerical experiment

Parameter	Notation	Value
Horizontal cell size	$ \Delta x $	100 m
Vertical cell size	$ \Delta z $	10–100 m
Maximum water depth	$ {H}_{\mathrm{m}\mathrm{a}\mathrm{x}} $	4 400 m
Time step	$ \Delta t $	5 s
Starting model time	$ {T}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}} $	00:00 UTC 17 April 2021
Predicting time	$ {T}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}} $	5 d
Horizontal eddy viscosity coefficient	$ {A}_{{\rm{h}}} $	10–1 m2/s
Vertical eddy viscosity coefficient	$ {A}_{{\rm{v}}} $	10–3 m2/s
Horizontal diffusivity coefficient	$ {K}_{{\rm{h}}} $	10–1 m2/s
Vertical diffusivity coefficient	$ {K}_{{\rm{v}}} $	10–3 m2/s
Bottom stress	$ {C}_{{\rm{d}}} $	2.5×10–3

 | Show Table

DownLoad: CSV

Analogous to the stripe brightness in the satellite images, numerically-predicted surface height gradients can measure the surface roughness induced by ISWs. A model snapshot was selected at 01:00 UTC on 20 April (Fig. 1c). Three stages of northward-going ISWs were clearly identified, as well as a southward-going ISW packet, whose locations were reasonably consistent with those in the satellite image (Fig. 1b). To a certain extent, it verified the model accuracy. Moreover, a transect (red dashed line in Fig. 1c) along the wave propagation directions was selected to demonstrate the vertical structures of ISWs from generation in the Lombok Strait to the shoaling process at the wreckage site (Figs 2 and 3).

Figure  3.  Numerically-predicted sea surface height gradients ($ \left|\nabla \eta \right| $) (a) and numerically-predicted isotherms ($ T $) from 28.5°C to 9°C (bending line in b, d and f from top to bottom) with a interval of 1.5°C and wave-induced velocity ($ {U}_{\rm {bc}} $, color shade) along the transect (red dashed line in panel a) (b) during the period from 19:30 UTC 20 April to 22:30 UTC 20 April with the time interval of 1.5 h. It characterizes the ISW properties right before and after it passing the wreckage site. The wreckage site is marked as the magenta star. z represents water depth.

DownLoad: Full-Size Img PowerPoint

The model results demonstrated that ISWs, generated from the Lombok Strait, presented different characteristics on three stages, i.e., with the amplitude of 90 m, 50 m and 40 m on the generation, propagation and shoaling processes (Figs 2b, d and f). Maximum ISW amplitudes mainly occurred at the water depth between 150 m to 400 m, covering the submarine’s collapse depth of ~200 m, thereby dragging it down to a more dangerous depth. At 21:30 UTC on 20 April 2021, an ISW packet with a leading wave amplitude of 40 m was shoaling and passing the shipwreck site (Fig. 2f). This time is in coincidence with the reported missing time. The ISW packet had a long characteristic width of approximately 50 km and remained fluctuations for over 10 h. Even though the following waves have relatively small amplitudes of 10−30 m, the continuous wave motions are likely to have a sustained impact on the submarine (Fig. 3). It is noteworthy that the future submarine motion should be more careful when passing the Lombok Strait (Stage 1), where the local ISWs might have a reasonably large amplitude of 90 m (Fig. 2b) in the ocean interior. In future, when submarines sail across the Lombok Strait, the voyage depths would be better in the upper 150 m, where the ISW amplitudes are relatively small, so the ISW is unable to drag the vessel down to a dangerous depth.

In summary, intense ISW events in the Lombok Strait have remarkable vertical displacements within a few minutes, thereby significantly affecting the underwater navigation and action of submarine. This numerical study concludes an ISW packet with a maximum amplitude of 40 m at the wreckage site on 20 April (UTC) 2021, which is likely the culprit to the sunk KRI Nanggala-402 submarine. Although satellite observations and numerical modelling have illustrated the crucial role of ISWs in the ocean interior, in situ observations of ISWs are needed to tell a more complete story in the future.

Acknowledgements

We acknowledge the use of MODIS-Aqua imagery from the NASA Worldview application (https://worldview.earthdata.nasa.gov). The numerical simulation is supported by the High Performance Computing Division and HPC managers of Wei Zhou and Dandan Sui in the South China Sea Institute of Oceanology.