1 Introduction

Large-amplitude internal solitary waves (ISWs) are frequently detected by in situ and remote sensing observations in the northern South China Sea (SCS)^[1-4], where they are generated from the eastern Luzon Strait and eventually dissipated on the western continental shelf after hundreds of miles evolution^[5]. ISWs are thought to have significant effects on biological productivity and sediment transport^[6]. Strong currents associated with large-amplitude ISWs make them potentially hazardous to underwater navigation and offshore drilling^[7].

Satellite images have been proved to be very helpful in investigating ISWs in the SCS^[8-12]. By examining the spatial characteristics of 116 internal wavepackets in satellite images from 1995 to 2001 in the SCS, Zhao et al.^[9] identified two types of internal wavepackets. In the deep basin, the wavepackets contain only one solitary pulse with or without oscillatory tails, while on the South China shelf, all the internal wavepackets are composed of a group of rank-ordered ISWs. As part of the nonlinear internal waves initiative (NLIWI), an in situ measurement in the northern SCS with 10 moorings spanning Luzon Strait, through the deep basin to the upper continental slope, Alford et al.^[3] observed the similar dynamics. In spite of satellite and field observations, many issues are still not clear in the SCS, especially disintegration of ISWs and formation of wavepackets.

With the advantage of modeling the wave evolution in the SCS under various conditions, the weakly nonlinear Korteweg-de Vries (KdV) theories^[13] play a primary role in elucidating the essential features of the observations^[14-15]. The bottom topography can cause the dramatic deformation and transformation of ISWs, which has been studied for a long time^[16-17]. If the depth varies very sharply, the ISW passes the transition zone without the change in the amplitude and shape, then it evolves into a finite number of solitary waves with an oscillatory tail. The outcome can be predicted explicitly based on the inverse scattering transformation^[18]. If the depth varies very slowly, the ISW undergoes a stage of adiabatic decay, and the evolution process can also be predicted by theoretical results^[19-20].

However, under the real ocean condition, as for the topography between two extreme cases mentioned above, there have been mainly numerical studies. Grimshaw et al.^[21] investigated the rotation effect on nonlinear internal waves with the Ostrovsky equation^[22] and examined his study with a series of laboratory experiments. The combined effects of the variable topography and background rotation on shoaling oceanic ISWs in the SCS were also studied^[5]. Although their calculations qualitatively agreed with the predictions of the theoretical development for the Ostrovsky equation, the amplitudes of the leading waves were about 150 m 250 m at the water depth of 300 m. These waves were found to be unrealistically large and may resulted from the absence of dissipation as they pointed out in their article. The propagation of an ISW in the ocean is accompanied by energy loss due to bottom friction^[23-25]. Chezy friction, parameterized by a quadratic term, is most commonly used in modeling the turbulent bottom boundary layer (BBL) with a suitable empirical drag coefficient. Of noteworthy is the work of Miles^[26] and El et al.^[27] who both gave a typical value of 0.01 for the oceanic ISWs, Grimshaw et al.^[28] who provided an empirical drag coefficient of 0.002 5, Henderson^[29] who found that the bottom stress can be well represented by a simple Chezy friction model with the coefficient of 0.001 5, and Wuest and Lorke^[30] who drew a conclusion that the empirical coefficient in the lake ranges from 0.001 to 0.003. The drag coefficient is usually assumed to be a constant value or within a certain range. However, there is no specific guideline for selection of the coefficient. Hence, with the aim of providing a method to assist the determination of the drag coefficient, the present work is motivated.

In this paper, the propagation and disintegration of ISWs in the SCS are studied by taking into account weakly nonlinear and dispersive effects, variable topography and bottom friction. The theoretical model is briefly presented in Section 2. In Section 3, the density stratification and topography in the SCS are used to study the friction effect, and the drag coefficient is fully discussed based on the field observation and theoretical analysis. Next, numerical simulations are conducted, and the amplitude of the leading ISW during the evolution process is discussed in Section 4. All results are summarized in Section 5.

2 Model description

The KdV equation has been widely used in describing the ISWs observed in the SCS^[31-33]. For a continuous stratified fluid, it can be written as

(1)

where η is the vertical displacement of the pycnocline, x is the horizontal coordinate, and t is time. The coefficients c, α, and β are the linear long-wave speed, the nonlinear coefficient, and the dispersion parameter, respectively. For simplicity, the Boussinesq approximation and the rigid-lid approximation are adopted here. The wave speed c is determined from the Taylor-Goldstein equation in the long-wave limit^[34],

(2)

where z is the vertical coordinate, U(z) is the horizontal background shear current, N(z) is the buoyancy frequency, and H is the water depth. The modal function ϕ is normalized with its maximum. The nonlinear coefficient α and the dispersive coefficient β are given by^[35]

(3a)

(3b)

(3c)

Equation (1) is valid only for a constant-depth fluid. However, when the water depth varies slowly in the horizontal direction, an additional term should be included^[36]. Therefore, the new equation reads

(4)

where

(5)

The term Q represents the amplification factor in the linear long-wave theory due to the variable topography, and the terms with the subscript 0 are evaluated at some fixed point x₀ (it is set to the initial position of an ISW in this paper).

The bottom friction term is then introduced into Eq. (4) to inhibit the unrealistic growth of wave amplitude on the continental shelf area and can be expressed as^[37-38]

(6)

where C_D is the empirical drag coefficient, and P can be expressed as

(7)

After the transformation of variables^[37],

(8)

where s is a temporal variable, and ζ(x, s) is the wave amplitude function. To the same leading order of approximation, Eq. (6) reduces to

(9)

Introduce the nondimensional variables

(10)

into Eq. (9). Then, the following nondimensional equation can be obtained:

(11)

It can be seen that Eq. (11) has exactly the same form with Eq. (9).

Following Holloway et al.^[23] and Grimshaw^[15], the mass and energy of ISWs can be written as

(12)

(13)

where the integrals are calculated over a wave period. Similar to Miles^[26], the mass-transport and energy-transport equations can be obtained from Eq. (11),

(14)

(15)

The relations (14) and (15) are often used to validate the numerical algorithm of the generalized KdV equation.

3 Drag coefficient 3.1 Theoretical analysis

The amplitudes of ISWs observed in the deep basin in the SCS vary from less than 50 m^[39] to more than 200 m^[40]. Although there are large differences in the amplitude and energy in the deep basin, they follow similar traveling paths from the deep sea to the continental shelf, and thus their propagation distances are approximately equal. It implies that ISWs are dissipated over the same distance in the SCS, and this process is independent of the initial states. The nondimensional energy E_N can be defined as

(16)

where E₀ is the initial energy of the ISW.

The relation of the drag coefficient and the initial amplitude of ISW is studied for the flat bottom. When the fluid is of constant depth, Eq. (11) can be rewritten as

(17)

where ₁, ₁, ₁, and Q₁ are constants during the evolution process. We assume that the amplitude of the ISW varies slowly with . Then, the leading-order approximation to can be written as

(18)

where V is the linear phase speed, and () is the corresponding amplitude of the ISW varying with . During the evolution process, Eq. (15) is always satisfied. Substituting Eq. (18) directly into the energy-transport Eq. (15) with the initial condition

(19)

yields

(20)

with ± depending on whether the ISW is depression or elevation. By substituting Eqs. (18) and (20) into Eq. (16), the nondimensional energy can be expressed as

(21)

We propose the hypothesis that different waves have the same nondimensional energy after propagating for a certain distance, and it follows that E_N is independent of its initial amplitude ₀, namely, C_D0 should be a constant, which indicates that the drag coefficient is inversely proportional to the initial amplitude of an ISW. This hypothesis will be verified in the subsequent section by comparing with the field observations.

3.2 Field observation and numerical simulation

Alford et al.^[3] reported a field measurement conducted in the SCS by using the NLIWI. 10 physical oceanographic moorings were deployed from Luzon Strait to the continental slope (see Fig. 1(a)). During the NLIWI experiment, 8 ISWs were detected, and their amplitudes at MP1 and LR2 are shown in Fig. 2. In the present section, the simulations of ISWs propagating from mooring MP1 to LR2 are carried out, and the results are compared with the in situ observation.

The averaged topography profile along the mooring track is shown in Fig. 1(b), which is taken from ETOPO1^[41]. The typical topography presented in Fig. 1(b) can be split into four parts. The segment from A2 to A3 is the deep ocean part characterized with the large water depth and the slowly varying topography. The shelf slope region from A3 to LR1 features a drastic change of the water depth about 2 000 m in a distance of 100 km. The shelf break from LR1 to LR2 is a transition region between the eastern shelf slope and the western continental shelf. The continental shelf from LR2 to MP2 is a shallow gently sloping region. The buoyancy frequency shown in Fig. 1(c) is derived from the summer climatological World Ocean Database 2013^[42]. The background stratification is taken to be spatially uniform, and background currents are not included. Therefore, the coefficients α, β, c, Q, P, and I can be calculated based on the given topography and background stratification.

Since the sign of the nonlinear coefficient does not change along the track line, the polarity reversal is not discussed. In a similar research region in the SCS, Grimshaw et al.^[5] pointed out that the cubic nonlinear coefficient is quite small compared with the quadratic nonlinear coefficient. Besides, over the sufficiently long propagation path, the rotation plays an important role in the dynamics of shoaling ISWs. The propagation of ISWs in the SCS over a length of 300 km, which is approximately the same distance in this paper, is studied by Xie et al.^[33]. They concluded that the rotation plays a minor role in comparison with nonlinearity. Therefore, the cubic nonlinearity and background rotation are not considered in the present paper.

Equation (9) is discretized with Fourier modes^[43] in s and integrated with the standard split-step method^[44] in the evolution variable x, and the 2/3-dealiasing strategy is used for the nonlinear terms. The resolutions of s and x were chosen with care to give appropriately resolved solutions, and in most calculations, Δs=100 s and Δx=10 m were used in dimensional variables.

Simulations in the present section are initialized with a single, first-mode ISW propagating westward from the slope shelf MP1, namely, the initial condition of a well-developed ISW is given by

(22)

where ζ₀ is the amplitude of ISWs at MP1, obtained from the in situ observation shown in Fig. 2. The coefficients α₀, β₀, Q₀, and c₀ are the localized parameters at mooring MP1.

The coefficients c, α, β, , and along the track from MP1 to LR2, shown in Fig. 3, are calculated based on the bottom topography and background stratification. From MP1 to LR2, the water depth changes from 1 500 m to 420 m in a distance of 74 km. The linear long wave speed c decreases to half of its initial value, while the nonlinear coefficient α is quite stable along the track. Compared with other parameters, the dispersion coefficient β drops most (it is about one twentieth of its initial value at LR2). The amplification factor increases monotonically from 1 to 1.5 at the location x=74 km.

The bottom friction is introduced in the evolution model to obtain more realistic amplitudes of ISWs in the SCS. A first-mode ISW with the amplitude η₀ detected at mooring MP1 is chosen as the initial condition. Then, the drag coefficient C_D is determined when the numerical solutions at LR2 match the field observation in the amplitude. The amplitudes of ISWs at MP1 and the corresponding drag coefficients are shown in Fig. 4. There is an approximately inverse relationship between the drag coefficient C_D and the initial amplitude, which can be expressed as

(23)

where ₀=η₀/H₀, and H₀ is the depth at the initial position. There is a dot far from the fitting curve. We suppose that this is probably an inaccurate dot. Without this dot, the constant on the right side of Eq. (23) is -7.7×10^-4. Anyway, Eq. (23) supports the assertion made in the last section, and furthermore, it indicates that although the fact that C_D0 is a constant is obtained in a constant depth, this relation is also satisfied in the SCS when an ISW experiences a drastic change of water depth from MP1 to LR2 as shown in Fig. 4.

4 Evolution of ISW in realistic oceanic environment 4.1 Friction effect

The observation of the generation and evolution of ISWs in the SCS suggests that the first-mode internal waves were first detected at A2^[3]. Hence, in this section, we deal with a single ISW evolving from the deep basin A2 to the middle of the upper continental slope MP2. As shown in Fig. 5, this region covers the deep basin (A2-A3), shelf slope (A3-LR1), shelf break (LR1-LR2), and continental shelf (LR2-MP2), and the bottom depth decreases approximately from 3 000 m to 300 m monotonically in a distance of 308 km.

Equation (9) with the coefficients calculated by Eqs. (2), (3), (5), and (7) shown in Fig. 5 is used in studying the propagation and disintegration of ISWs in the SCS. Along the track line from A2 to MP2, the linear long wave speed c decreases from 3 m/s to 1.1 m/s, and the dispersion coefficient β decreases approximately from 10⁶ m²/s to 10⁴ m²/s. For a KdV-type soliton, the characteristic width of an ISW can be expressed as , which means that the decreasing of the dispersion coefficient is a dominant factor in reducing the width of the ISW from thousands of meters in the deep basin to hundreds of meters on the continental shelf. The spatial variation of the nonlinear coefficient α is small in comparison with the variation of the dispersion coefficient. The amplification factor Q mainly changes at the shelf slope, while it is almost constant in the deep basin and in the continental shelf region.

The evolution of a depression KdV soliton is numerically investigated with and without the bottom friction. In Fig. 6(a), an initial ISW with the amplitude of 80 m is released at A2 with t = 0 h at the wave trough. The drag coefficient C_D is set to zero in this simulation. Figure 6(a) shows the time series of wave displacements at selected locations with the gap of 15 km. It is not until MP1 that the topography shows a strong impact on the incident solitary wave. The amplitude of the leading solitary wave is doubled as it propagates to LR1. At the same time, the solitary wave undergoes a disintegration process, with a shelf forming behind. As it further evolves to MP2, the amplitude of the leading wave is more than 300 m, which is nearly four times the initial amplitude. The behind shelf separates from the leading solitary wave and evolves into a trailing wavepacket. Similar simulations without friction were conducted by Grimshaw et al.^[5] in the SCS. They observed an ISW whose amplitude evolves from 64 m initially in the deep basin to approximately 250 m on the continental shelf, and claimed that the unrealistically large wave was due to the absence of dissipation. Figure 6(b) shows the simulation results conducted with the drag coefficient determined by Eq. (23). There is no obvious difference between the cases with and without friction until MP1, which means that the friction has few impacts on the solitary wave in the deep basin. From MP1 to LR2, the amplitude of the leading solitary wave increases by approximately 100%. After that, the amplitude decreases significantly, being half of its initial value as the wave arrives at MP2. The ISW also endures a disintegration process where a trailing wave tail is formed due to the interaction with the bottom topography. This disintegration phenomenon is consistent with the satellite observations of Zhao et al.^[9] and Jackson^[10]. It is found that the bottom friction is very important in preventing wave amplitudes from growing unrealistically large.

4.2 Amplitude variation

The amplitude is the most important quantity in characterizing an ISW. From the numerical simulation, the ISW in the SCS, under the combined effects of topography and friction, undergoes a sophisticated evolution and a disintegration process and finally forms a wavepacket on the continental shelf. The horizontal current induced by the ISW, in the uniform depth, could be expressed as^[24]

(24)

where c, η, and ϕ are defined in Section 2. Therefore, ISWs with the larger amplitude can induce stronger currents and have larger potential threats on offshore structures. Besides, the leading solitary wave has the largest amplitude in the wavepacket. Therefore, it is useful to study the amplitude variation of the leading solitary wave. The water depth is decreasing monotonically from A2 to MP2, and the amplitude of the solitary wave is increasing at the same time as shown in Fig. 6(a). Therefore, the shoaling effect leads to the increasing of the wave amplitude. On the contrary, the effect of bottom friction tends to suppress the growth of the wave amplitude. Similar to the definition of the Ostrovsky number characterizing the ratio between nonlinearity and rotation^[45] and the Ursell number characterizing the ratio between nonlinearity and nonhydrostatic dispersion^[39], a nondimensional number characterizing the ratio between the shoaling effect and the friction effect can be defined as

(25)

For a solitary wave with T>1, shoaling inhibits dissipation, and the amplitude of the leading wave increases, whereas for T < 1 the friction effect dominates, and the amplitude of the leading wave decreases. A special example is the flat bottom, where T=0 leads to the decreasing of the amplitude of the solitary wave during the evolution process. After Eq. (23), C_D|η| can be approximately assumed a constant along the propagation path. The coefficients c, I, P, and Q along the propagation path are only dependent on environment variables according to Eqs. (2), (3), (5), and (7). Therefore, T is only dependent on environment variables and can be calculated just based on the bottom topography and buoyancy frequency.

Along the propagation path as shown in Fig. 7(b), three numerical simulations with different initial amplitudes and the corresponding drag coefficients are carried out. Figure 7(c) shows the variation of the nondimensional parameter T along the path. The area with T bigger than 1 which is displayed with shadows indicates that the wave amplitude increases in this region, while decreases in the rest part. The wave amplitude A (normalized by its initial value a₀) along the traveling path in three experiments is shown in Fig. 7(a), which displays a similar trend in variation. Five regimes are divided based on the nondimensional parameter T as shown in the same figure. Before L1, the amplitude is stable with a slight increase of 2%. From L1 to L2, the amplitude increases, then from L2 to L3 it decreases for a short distance. From L3 to L4, the amplitude increases again and reaches its maximum. Finally, after L4, the amplitude experiences a continuous reduction process. Before L1 and after L4, three curves almost coincide, since the bottom topography is almost flat in the two regions, and the relationship between the drag coefficient and the initial amplitude is quite satisfied. However, there are differences among three curves from L1 to L4. This may result from the drastic changes of topography and the disintegration of ISWs in this region. The similar law was reported by Alford et al.^[3] about the amplitude variation of ISWs in the SCS, who documented that the wave amplitude increases from MP1 to LR1 and reaches it maximum near LR2 and finally decreases in LR3. By comparing Fig. 7(a) with Fig. 7(c), a conclusion can be made that the nondimensional parameter T is useful in predicting the variation of the leading solitary wave.

5 Summary and conclusions

While the topography effect on the generation and propagation of ISWs in the SCS has been extensively studied and documented^{[3, 5, 33]}, and the effect of bottom friction on solitary waves has been separately studied^{[29, 38, 46]}, their combined effects have not previously been investigated in detail in the SCS. The generalized KdV equation is presented as a basis to develop the theoretical analysis which is testified by the observation of ISWs in the SCS by Alford et al.^[3].

Based on the real topography and typical buoyancy frequency, the present simulations show that, if we consider the topography variation, the shoaling effect will cause the unrealistically large amplitude ISWs on the continental shelf, which was also pointed out by Grimshaw et al.^[5]. The Chezy friction has been most commonly used and proved to be useful in modeling the dissipation process in the real ocean. Its coefficient C_D is an empirical number and often assumed to be a constant. However, through the theoretical analysis of a case of constant depth, an inverse relationship between the initial amplitude and the drag coefficient is found. Moreover, numerical simulations based on the realistic oceanic environment also suggest that this relationship can be approximately satisfied in the SCS. Moreover, a quantitative formula between the amplitude of ISW and the drag coefficient is given.

The shoaling effect leads to the increase in the amplitude of ISWs in the SCS while the friction effect tends to inhibit their growth. A nondimensional parameter T defining the ratio between the shoaling effect and friction effect is put forward to take into account both effects. Numerical results show that T is a good indicator in determining the trend in variation of the amplitude of ISWs in the SCS. In conclusion, the topography in the SCS can be divided into three parts.

(Ⅰ) The deep basin region Both the topography and friction effects are weak in the deep basin, and the ISW can keep an almost permanent waveform.

(Ⅱ) The slope shelf and shelf break region The topography varies drastically on the shelf slope. Thus, the shoaling effect is the dominant factor in determining the wave amplitude. As the wave amplitude increases, and the topographic slope decreases, the friction effect is strengthened and gradually balances with the shoaling effect. Consequently, the ISW reaches its maximum amplitude at the shelf break.

(Ⅲ) The continental shelf region The interaction between the ISW and bottom topography is strong due to the largeness of the amplitude and shallowness of the water. Besides, the topography varies slowly in this region which weakens the shoaling effect. Therefore, the bottom friction effect is the dominant factor affecting the wave amplitude on the continental shelf.

The discussion in the context relies on the KdV-based interpretation of physics. However, the weakly nonlinear theory may have some restrictions in interpreting extremely large-amplitude ISWs, for example, the observation of Lien et al.^[47-48]. This restriction can be avoided through some fully nonlinear equations based on the Euler equation for steady solitary waves ^[49-50]. The turbulent BBL is the source of Chezy friction. It is customary to represent the overall dissipation caused by the turbulent BBL with a Chezy friction term in the KdV-type equation^{[23, 27]}. There are also other types of dissipation, like shear and convective instabilities in the interior of ISWs^[51], and this shall be further investigated in the future.