Introduction

Large amplitude nonlinear internal solitary waves (ISWs) are ubiquitous features in the ocean (see www.internalwaveatlas.com). The Northeastern (NE) South China Sea (SCS, Fig. 1) is a region where large amplitude internal waves are frequently observed1. The observed ISWs in NE SCS are generally considered to be generated at the Luzon Strait, the primary water passage connecting the SCS and the Pacific Ocean. They are caused by the interaction between the strong tidal forcing and the abrupt topography located there2,3. These internal waves then propagate westward crossing the deep sea basin region, pass the sharp slope area, entering the shoaling shelf zone and ultimately dissipate on the continental shelf. The propagation and evolution processes and mechanisms of the ISWs as they cross the SCS are rather complicated due to strong interactions between currents, eddies, fronts, internal waves, stratification and topography2,4,5. These large amplitude ISWs induce strong currents which have substantial but not well-studied effects on oceanic mixing, biological processes, sediment transport, offshore oil engineering and submarine navigation1,6,7.

Figure 1
figure 1

Bathymetry of the NE South China Sea.

The purple line shows the seismic section carried out from 5:30 to 21:50 on 27 June 2012 (local time, UTC + 8:00) with heading of 65 degrees from the north. Red dots and the adjacent numbers mark distances from the starting point. Black dashed curves are surface signatures of internal solitary waves from satellite images3. Red lines are the imaged solitary waves in the study. Figure produced with the Generic Mapping Tools (GMT). Bathymetry data is from http://topex.ucsd.edu.

Full size image

Although their importance is well recognized, studies of ISWs are far from being comprehensive. In-situ hydrographic observations of the ISWs are still rare and many studies rely only on remote sensing and numerical modeling1,8.One new technique that shows promise in observing small-scale features over large areas and full depth of the ocean is seismic oceanography, which is capable of capturing finescale structures in the ocean with a spatial resolution of O (10 m) using standard seismic reflection profiling9. In less than 10 years, the seismic reflection method has been used to image numerous interesting oceanic phenomena, such as internal waves10, subsurface eddies11,12,13, horizontal vortices14 and subsurface currents15. Nevertheless, the spatial features of near-surface ISWs remain unexplored by the seismic reflection method due to the large source-receiver distance of the seismic setups and the significant noise of source bubble oscillations.

In order to image the near-surface water column structures in the NE SCS (Fig. 1), an optimized seismic observational system was deployed on 27 June 2012 (Methods). Two prominent near-surface ISWs were clearly seen on the seismic image in tandem for the first time (Fig. 2). They were also captured in a satellite image obtained for the same day. Propagation velocities of the ISWs derived from the seismic, satellite and on site observation data are all consistent. Our results suggest that the seismic reflection profiling, in conjunction with conventional hydrographic observations, could be an excellent method to study near-surface ocean processes such as ISWs.

Figure 2
figure 2

Seismic image of the water column in the NE SCS.

(a) Seismic image of the section in Figure 1 (only the portion of horizontal and vertical ranges of 15–70 km and 0–300 m, respectively, is shown). Black triangles are the locations of XBT stations. Detailed structures in the two white boxes of (a) are expanded in (b) and (c), respectively. Figure produced with Seismic Unix (SU).

Full size image

Results

General features of the water column

A clear seismic image of the upper ocean thermocline structure between ~50 to 300 m in depth was obtained along the continental slope in water depths of 400–1000 m (Figs. 1 and 2). Strong features in the images and large impedance variations estimated from simultaneous XBT profiles are in high correspondence (Supplementary Fig. S1). Further comparisons with hydrographic measurements/derivations show that temperature is the dominant factor of the seismic response (Supplementary Fig. S2). Therefore, it is reasonable to view the seismic reflection surfaces as being parallel to isotherms in high-vertical gradient areas.

The most pervasive feature of the seismic image are coherent reflections from quasi-horizontal linear features of 5–10 km in length, which include undulations with a vertical extent of around 20 m and a horizontal extent of less than 1 km (Fig. 2a). These are likely due to a background field of small internal waves10. These internal waves are more continuous in the deep slope region at the right side of the image than in the shallow slope region on the left side, where the internal waves are intermittently interrupted by patches of blanking zones a few km in horizontal extent and 50–100 in vertical extent, possibly indicating regions with well-mixed water. The measurement from XBT1 within such a patch confirms these blanking zones coincide with regions where there is a lack of strong gradients in the vertical variation of the hydrographic components (Supplementary Figs. S1a and S2a). The interactions among eddy currents, topography, internal waves, or other nonlinear processes are likely responsible for these well-mixed patches.

Internal solitary waves captured by seismic survey and satellite

In addition to the field of small internal waves, two isolated regions showing internal vertical displacements of O(50 m) over a horizontal scale of less than 1 km were also captured during the cruise on 27 June 2012 (Fig. 2). These displacements are interpreted as being the large leading waves in packets of internal waves. Each wave packet is dominated by one large depression followed by smaller waves (Fig. 2b,c). The vertical influence of the initial depression extends from 50 m to deeper than 300 m (Fig. 2b). These are typical morphological features of ISWs as described elsewhere by hydrographic observations or model descriptions6,16,17,18. The amplitudes of the features are smaller than typically observed for ISWs in this region2,7,16,19. However, our observations occurred during neap tides, when ISWs would have smallest amplitude (Supplementary Fig. S3). Hereafter, these two wave packets are named as ISW1 and ISW2, respectively. Neap tides are also a time when low-frequency undular bores have been observed in this region2,20 and such a bore may be responsible for the plunging layers between 150 m and 200 m at the left edge of Figure 2b.

A NASA Terra satellite image (Fig. 3), acquired on 27 June 2012 by the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument, shows a number of dark bands forming curving wavefronts with a lateral extent of up to 100 km in the general vicinity of the seismic track. These appear to come from two packets, whose leading edges cross the seismic track at about −3 and 75 km. By identifying these wave crests with the two large displacement features seen in the seismic section, the space-time differences between the seismic and satellite observations can be used to estimate mean phase velocities for ISW1 and ISW2 to be 1.75 m/s and 2.33 m/s, respectively (Table 1).

Table 1 Parameters observed/derived during velocity estimation. Tsat: satellite acquisition time; Tseis: seismic acquisition time; LT: local time, equals to UTC + 8:00; Dist: ISW distance between satellite and seismic observation; A: angle between ISW propagation and seismic line heading; η0: amplitude of ISW. Vsat: velocity from satellite observation; Vsa: apparent phase velocity from seismic data along the seismic line; Vseis: phase velocity from seismic data along the wave propagation direction; c: mode-one linear phase speed; KdV: speed estimated using continuous stratified ocean system; Vvs: vertical mean velocity from seismic data; Vvx: vertical mean velocity from hydrographic data
Full size table
Figure 3
figure 3

Filtered MODIS image.

The observed ISWs at 11:05 on 27 June 2012 (local time, UTC + 8:00) in a MODIS image (http://oceancolor.gsfc.nasa.gov). Two red squares mark the locations of ISW1 (left) and ISW2 (right) when observed on the seismic image. Figure produced with the Generic Mapping Tools (GMT).

Full size image

Mutual verification among seismic, satellite and site observations

In addition to using the space-time difference between the seismic and satellite observations (named as approach A), phase velocities of ISW1 and ISW2 are also obtained using two other approaches. These include an analysis of the transient velocities from pre-stack seismic data (named as approach B) (Supplementary Figs. S4–S6)14,21 and a theoretical speed derived from hydrographic data based on the Korteweg-de Vries (KdV) equation (named as approach C)6,17,22.

In approach B, the ISW trough is easy to identify in signals from different source-receiver pairs and wave phase speed can be inferred from these observations (Supplementary Fig. S4). By assuming a constant propagation speed and aligning the common offset gathers, an in-plane propagation speed can be estimated (Supplementary Fig. S5). Apparent velocities along the seismic line are 1.88 m/s and 2.67 m/s for ISW1 and ISW2, respectively (Table 1). Then using the angle between the ISW propagation and the seismic line estimated from the satellite image (Table 1; Fig. 3), the apparent velocities in the seismic line direction are projected into the wave propagation direction and the actual wave phase velocities are estimated to be 1.70 ± 0.06 m/s and 2.19 ± 0.34 m/s for ISW1 and ISW2, respectively (Table 1). These are consistent with the estimated wave speeds using approach A.

In approach C, the data from XBT1 and XBT2 are used to derive the analytical parameters of ISW1 and ISW2, respectively, because the distance of each XBT-ISW pair is very close. With no background current and a continuously stratified ocean, the mode-one linear wave phase speeds for ISW1 and ISW2 are estimated as 1.40 and 1.68 m/s, respectively (Table 1). These wave speeds can be corrected for nonlinear effects by assuming that a continuously stratified Korteweg-de Vries equation adequately describes the propagation. The nonlinear phase speeds are 1.63 m/s and 1.92 m/s. Here, the estimated phase speed of ISW1 is comparable to the observed values while the speed of ISW2 is about 12% lower than the seismic observation, although still within the uncertainty of the seismic observation. This discrepancy between the theoretical estimates and the observed phase speeds may be caused by the neglected background currents19,23, although satellite altimetry-derived velocity fields (not shown) suggest that the mean current at least is approximately perpendicular to the direction of propagation. Another explanation for the difference is that ISW2 appears to be transforming its shape and is not a pure, symmetric solitary wave.

Other important parameters during the wave propagation stage such as waveforms and vertical velocities can also be derived from the seismic data and compared with analytical estimates from KdV theory.

While the seismic image directly provides the waveform of ISW by tracking the reflective horizons, the hydrographic data can also predict the wave shape using the theoretical models. In this study, the ratios of the amplitude to the upper-layer depth are around 0.7, so a weakly nonlinear theory16, e.g., the KdV theory may be sufficient. The analytical solitary solution of the classical KdV equation is η(x, t) = η0sech2((x − vt)/Δ), where the nonlinear velocity v and the characteristic width Δ are related to the mode-one linear speed c, the nonlinear parameter α, the dispersion parameter β and the maximum wave amplitude η0 by v = c + αη0/3, Δ2 = 12β/αη0 (Supplementary Table S1). The ratios of the seismic observed full widths at the half wave amplitude to the characteristic widths are 1.71 and 1.73 for ISW1 and ISW2, respectively, while the analytical ratio is 1.76. Therefore, the KdV theory predicted solitary wave shapes match well with the seismic observed waveforms (Fig. 4).

Figure 4
figure 4

Waveform comparison of ISWs.

Two predicted waveforms of ISWs (white lines, derived using the stratified KdV theory from XBT1 and XBT2 profiles) are superimposed on their seismic images. The predicted waveforms have been contracted by simulating the Doppler-Effect for wave-ship movements and then projected to the direction of seismic section. Figure produced with Seismic Unix (SU).

Full size image

Similarly, the vertical velocities of the observed ISWs are also compared with analytical predictions of the KdV theory. From the seismic image, the vertical velocity between the peak and the half peak of the displacements is measured using the relation w = η0/L (Supplementary Fig. S7). Estimated vertical velocities are 13.8 cm/s and 8.7 cm/s for ISW1 and ISW2, respectively. In comparison, KdV theory predicted values are 12.9 cm/s and 7.5 cm/s using the equation (Supplementary Fig. S8).

Discussion

Seismic profiling has previously been used to image oceanographic features at depths of 100 s to 1000 s of meters in the ocean. Less attention has been paid to features closer to the surface, because the conventional geometries of marine seismic survey are unfavorable for imaging regions that close to the surface. However, by optimizing the survey geometry, we have successfully formed a seismic image of the upper 300 m of the water column (but below the mixed layer) in the NE SCS. The image shows that thermocline fine-structure in this area is dominated by typical internal waves with a vertical extent of O(10 m) that can stretch for up to 10 km. These features are weakened or interrupted by patches a few km across where vertical profiles are smoother, verified by in-situ XBT observations. The most prominent features are the two distinct displacements in the section, associated with ISWs that form discrete packets radiating away from their source regions during neap tide when ISWs are rarely observed in this area2,7. The in-plane propagation speeds of these displacements were estimated from the pre-stack seismic data alone and the results are consistent with propagation speeds estimated in other ways. In addition, parameters like wave shapes and vertical velocities are successfully determined from the seismic data and can be comparable to the theoretical predictions.

Thus we have shown that studying ISWs becomes more comprehensive using an optimized seismic array and incorporating conventional observations. Other fine-structure features observed on the seismic image await confirmation. One of these features are the intriguing regions of possibly enhanced mixing occurring immediately after passing of the ISWs, identified by clear disappearance/weakening of the reflections around 200 m depth (Fig. 2), which may be a result of shear instability16. Although weak stratification in this area is more easily destroyed than the strong stratification above, whether such appearances are universal and whether they represent quantitative descriptions of mixing is unknown. Another feature requiring explanation is the flat-step structure of the ISW2 trough which is unable to be modeled by single KdV soliton. Its dynamic implication (wave splitting or merging, or other process) also remains to be assessed.

Therefore, we believe that continued studies of other interesting stages or phenomena of ISWs, such as generation, evolution, diffraction, dissipation, polarity conversion and mixing will benefit from coordinating seismic observations with other conventional techniques.

Methods

Seismic observation and processing

A seismic oceanographic experiment was conducted aboard the R/V SHIYAN2 in the NE South China Sea on 27 June 2012. In order to optimize for imaging near-surface features, a Sercel GI gun was chosen as the seismic source since it has good performance in suppressing bubble oscillations. However, due to the limited source energy from a single Sercel GI gun (with volume of 105/105 cu-inch), a modified acquisition geometry was necessary to improve the data quality. The nearest offset (source-receiver distance) was set to 90 m to reveal the shallow structures, which are overlooked by the conventional geometry in which the nearest offset is around 250 m. A Sercel Seal oil filled streamer with 72 receivers (12.5-m interval) was towed at 5 m depth. The air gun was triggered every 25 m (~10 s) to increase the fold number. The record length was 6 s with sampling interval of 2 ms. Such a geometry ensured a satisfactory ray sampling of the near-surface structure and a high SNR (Signal Noise Ratio) of the stacked seismic image in processing.

Seismic data were processed using the pre-stack depth migration (PSDM) method24,25. The conventional procedure of PSDM consists of six steps: (1) geometry definition, (2) bad trace deletion, (3) direct wave attenuation, (4) band-pass filtering, (5) PSDM with velocity analysis and (6) common image gather stacking. Here, the velocity analysis of the fifth step was omitted. Instead, the acoustic velocity was fixed at 1510 m/s according to our hydrographic observations. Migration tests have shown that it is an acceptable model with fine imaging. According to the seismic geometry and processing parameters in the study, the shallowest reliable seismic reflection is ~50 m depth, which is similar to the mixed layer depth from the XBT data.

Velocity estimation of the ISWs from seismic data

The propagation velocities of the ISWs can be estimated using the pre-stack seismic data because of the redundancy of the seismic acquisition spanning a finite time14,21,26. However, the near-surface reflection signals from ISWs tend to be muted by normal move out (NMO) and suppressed by direct wave attenuation. Thus the size of the common mid-point gathers (CMPs) or folds of the datasets decreased from 18 to around 7 traces (i.e., only the first 27 near-source traces of the common offset gathers (COGs) are left over). Therefore, the time and distance spans are too short to estimate the movement of the reflections precisely and intuitively according to Sheen's approach using the CMPs14. To overcome this limitation, we resort to the remained COGs which can be used to track the ISWs clearly and improve the approach accordingly. Integrity of the waveforms on the COGs helps us to derive the crest/trough points used for movement estimation.

The principle of the in-plane (apparent) velocity measurement is illustrated in Supplementary Fig. S4 with the assumption of uniform vessel speed. For the given coordinate system, horizontal locations are, cmp2 = −(g2*dg + ofs0)/2, cmp1 = −(g1*dg + ofs0)/2 − (s2 −s1)*ds, where (s1, g1) and (s1, g2) are the numbers of the shot-receiver pairs that capture the crest/trough of the moving internal wave at the midpoints cmp1 and cmp2, respectively; dg, ds, ofs0 are the receiver interval, shot interval and nearest offset, respectively. Thus, the horizontal velocity of the wave crest is measured by v = (cmp2 − cmp1)/T = −[(g2 − g1)*dg/2 − (s2 − s1)*ds]/[(s2 − s1)*dt], where T is the shot time duration between s1 and s2 and dt is the shot time interval. Taking ISW1 as an example, we first auto-tracked the waveforms from 27 COGs and then interpolated to find the troughs of the ISW and the corresponding source-receiver pairs. Then, with the known receiver span (g2 − g1) of 26, the only unknown shot span (s2 − s1) was estimated using the slope of the regression based on the RMS criterion (Supplementary Fig. S5). Finally, the apparent horizontal velocity v was derived, as well as the actual phase velocity on the direction of wave propagation. For the real data, the source range, time duration and ship heading were estimated precisely using the navigation data of the cruise.

To confirm the utility of this velocity estimation method, we estimated the velocity of a static seamount, whose spatial scale was comparable to ISWs, to perform the whole procedure to measure its velocity (Supplementary Fig. S6). The estimated velocity for the static seamount was only -0.04 m/s, suggesting any bias is small. A larger source of uncertainty is from the slope regressions and 95% confidence intervals estimated from the regression were ±0.06 m/s and ±0.34 m/s for ISW1 and ISW2, respectively. The primary factor of the large error of ISW2 is that the troughs of the ISW2 on COGs are too gentle to allow the trough to be accurately located.

Velocity estimation of the ISWs using XBT data

We used the in situ observations from XBT1 and XBT2 to calculate the nonlinear phase velocities of ISW1 and ISW2 correspondingly. Here, the salinities of the profiles were estimated using the historical CTD data in the study region from the World Ocean Database 200927. Thus, the density, velocity, impedance, buoyancy frequency were derived using the temperature/salinity/depth relationships. For a continuously stratified fluid, the first-mode linear velocities c were solved from the Sturm-Liouville equation22. Then, the parameters α and β of the one-dimensional KdV equation were calculated for the same stratification23, as well as the characteristic width Δ (Supplementary Table S1)6. Finally, the nonlinear phase velocities were estimated via v = c + αη0/3.