Numerical Simulation for the Evolution of Internal Solitary Waves Propagating over Slope Topography

: In this study, the propagation and evolution characteristics of internal solitary waves on slope topography in stratiﬁed ﬂuids were investigated. A numerical model of internal solitary wave propagation based on the nonlinear potential ﬂow theory using the multi-domain boundary element method was developed and validated. The numerical model was used to calculate the propagation process of internal solitary waves on the topography with different slope parameters, including height and angle, and the inﬂuence of slope parameters, initial amplitude, and densities jump of two-layer ﬂuid on the evolution of internal solitary waves is discussed. It was found that the wave amplitude ﬁrst increased while climbing the slope and then decreased after passing over the slope shoulder based on the calculation results, and the wave amplitude reached a maximum at the shoulder of the slope. A larger height and angle of the slope can induce larger maximum wave amplitude and more obvious tail wave characteristics. The wave amplitude gradually decreased, and a periodic tail wave was generated when propagating on the plateau after passing the slope. Both frequency and height of the tail wave were affected by the geometric parameters of the slope bottom; however, the initial amplitude of the internal solitary wave only affects the tail wave height, but not the frequency of the tail wave.


Introduction
Internal solitary waves are observed in many domains through synthetic aperture radar images and the measurement data of acoustic Doppler current profiler in oceans worldwide [1][2][3][4][5][6]. ISWs can propagate over long distances in the ocean owing to their large amplitude and concentrated energy [7], and they shoal with the influence of topography when propagating to the continental shelf area. The shoaling of ISWs upon sloping in oceans plays an important role in controlling the stratification and vertical distribution of biogeochemical matter in the water column [8][9][10]. In particular, the interaction between ISWs and the seabed is considered an important factor that promotes sediment suspension on the seabed and instability [11,12]. Meanwhile, strong ISW-induced shear flow during the shoaling process is a non-negligible threat to offshore engineering structures, such as oil drilling rigs, marine risers, and subsea pipelines [13,14].
The evolution characteristics of internal solitary waves on slope topography are very complicated, including various forms, such as breaking, fission, and polarity reversal. The evolution characteristics of internal solitary waves on slope topography have been studied extensively using experiments and numerical simulations aimed at changing the wave profile and flow field characteristics [15][16][17][18][19][20][21][22]. The breaking types of ISWs are studied using experiments and numerical simulations aimed at changing the wave profile and flow field characteristics [23][24][25]. However, ISWs can propagate over mild slopes without breaking, such as fission, polarity reversal, and the formation of trailing waves. Polarity inversion, wave fission and the forming of trailing waves are typical phenomenon while ISWs propagating over a gentle slope bottom. Most ISWs are depression waves due to stratification properties in the ocean, but the depression ISWs can change to elevation waves on the slope when the depth of upper layer is larger than that of the lower layer. Helfrich et al. conducted numerical simulations supplemented by inverse-scattering theory to investigate the change in polarity of the incident waves as they pass through the turning point of approximately equal layer depths, and the separation of the leading waves was observed [26]. Orr and Mignerey observed the conversion of nonlinear internal depression solitons to elevation internal waves through high-frequency acoustic flow visualization [27]. Meanwhile, the trailing waves following the initial waves can appear due to the dispersion effect while ISWs passing the slope topography. Shroyer et al. observed nonlinear internal waves propagating to the New Jersey coast and found the trailing face of the initial depression wave [28]. Zhi et al. found a trailing wave train when ISWs transforms over the slope using asymptotic theory and numerical simulation [29]. Lamb and Xiao conducted high-resolution two-dimensional numerical simulations of incompressible Euler equations to study the evolution of internal solitary waves shoaling onto a shelf, and the shoaling internal solitary waves generally fission into several waves [30]. The previous numerical studies on the evolution of ISWs propagating over slope bottom were mainly based on CFD method, and the boundary element method was rarely used to calculate these topics based on potential flow theory.
As an efficient numerical method, the boundary element method (BEM) is widely used in the calculation of water wave problems based on potential flow theory, which has excellent adaptability to complex boundaries due to the fact that only the boundary of the computational domain needs to be discretized. Longuet-Higgins and Cokelet introduced the boundary element method for the simulation of steep free surface waves [31]. Since then, BEM has been widely used for simulating the surface wave [32][33][34][35][36]. Meanwhile, comparisons of two-dimensional numerical results with laboratory experiments have shown that the potential theory can predict the characteristics of wave shoaling over slopes accurately [37][38][39]. Koo simulated the interaction of internal linear periodic waves with linear free surface waves using the BEM [40]. Gou et al. studied wave diffraction in a two-layer fluid using the BEM with linear boundary conditions [41]. Evans and Ford studied steady ISWs using an integral equation approach, assuming the rigid-lid condition at the top surface [42]. The existing investigations on internal waves using the BEM are mainly for linear periodic waves or steady solitary waves.
In this study, the evolution characteristics of ISWs propagating over slope bottom in a stratified fluid were investigated using a nonlinear numerical model based on the potential flow theory. The numerical simulations were conducted using multi-domain boundary element method which is an effective method to calculate evolution of water waves. The evolution of ISWs on variable bottom was rarely investigated using boundary element method in the existing investigation. A series of numerical cases were calculated with different slope parameters and initial wave amplitudes to study the relationship between the wave evolution characteristics and slope topography. The variations in wave amplitude, wave profile, and tail waves over different slope topographies were discussed.

Physical Model
A system of two-layer fluids whose densities are denoted by ρ k and thicknesses by h k , with k = 1 for the upper fluid layer and k = 2 for the lower layer, is shown in Figure 1. The bottom of the fluid system is a rigid slope seabed with different geometric parameters. Cartesian coordinates are defined with the x-axis at the mean position of the interface between the two liquid layers. The surface at the top of the upper layer is denoted as z = h 1 , and the interface between the two liquid layers is z = η 2 (x, t).
A system of two-layer fluids whose densities are denoted by k  and thicknesses by k h , with k = 1 for the upper fluid layer and k = 2 for the lower layer, is shown in Figure 1. The bottom of the fluid system is a rigid slope seabed with different geometric parameters. Cartesian coordinates are defined with the x-axis at the mean position of the interface between the two liquid layers. The surface at the top of the upper layer is denoted as 1 zh = , and the interface between the two liquid layers is (1) The kinematic boundary conditions at the interface are that the normal velocity is continuous at the interface.
( ) 12 2 , where n is the unit-outward normal of the domain considered. The dynamic boundary condition is that the pressure is continuous at the interface. Using the Bernoulli equation, this can be expressed as: The impenetrable boundary condition at the rigid seabed requires: The flows in the two layers are assumed to be irrotational and incompressible. The fluid velocities u k for the two fluid layers can, thus, be expressed as the gradient of the potentials ϕ k , that is, u k = ∇ϕ k . The velocity potentials satisfy the Laplace equation: The kinematic boundary conditions at the interface are that the normal velocity is continuous at the interface.
where n is the unit-outward normal of the domain considered. The dynamic boundary condition is that the pressure is continuous at the interface. Using the Bernoulli equation, this can be expressed as: The impenetrable boundary condition at the rigid seabed requires: The top of the upper liquid layer is set as a rigid lid: Using Green's second identity, the velocity potentials in the upper and lower fluid layers can be expressed as the line integrals along their boundaries Γ 1 and Γ 2 , respectively.
where r is the field point, q is the source point, G (r, q) = −1/(2π)ln|r − q| is the Green function, and c(r) is the solid angle of the boundary at r. Equation (6) can be further denoted as: where ψ k (q) = ∂ϕ k (q) ∂n .

Numerical Model
The boundary integral equation is solved using the boundary element method. Both boundaries Γ 1 and Γ 2 are discretized into n constant line elements Γ ij from point q j to q j+1 . Equation (7) are discretized are discretized using the constant elements as follows: With the constant elements used, we consider r j = 1 2 q j + q j+1 , r = r i in Equation (8), and we have c(r i ) = 0.5. In addition, we denote ϕ kj = ϕ k r j .
The two fluid layers have different densities such that the multi-domain boundary element method [43] should be adopted. For the two-layer computation domain, the boundary elements can be divided into two types, which are self nodes and common nodes, as shown in Figure 2. The boundary integral Equation (9) can be expressed as Equation (10).
where the superscript s and c denote self nodes and common nodes, respectively. The boundary conditions of the left and right boundaries Γ f are that fluid flux is zero as there is enough distance from solitary waves to the two boundaries. For the rigid-lid surface Γ t and bottom Γ b , impenetrable conditions are applied so the normal velocity is zero. Based on the pressure continuity and normal velocity continuity conditions at the interface, the two regions can be connected to solve the problem, and the multi-domain boundary element method is introduced to develop numerical model for the problem of internal solitary waves propagation and evolution. A specific calculation procedure can be found in the literature [44].
where the superscript s and c denote self nodes and common nodes, respectively. The boundary conditions of the left and right boundaries Γ f are that fluid flux is zero as there is enough distance from solitary waves to the two boundaries. For the rigid-lid surface Γ t and bottom Γ b , impenetrable conditions are applied so the normal velocity is zero. Based on the pressure continuity and normal velocity continuity conditions at the interface, the two regions can be connected to solve the problem, and the multi-domain boundary element method is introduced to develop numerical model for the problem of internal solitary waves propagation and evolution. A specific calculation procedure can be found in the literature [44].
In fact, the normal velocity is continuous, but velocity potential of two-layer fluid is not continuous at the interface. Therefore, the interface between the two liquid layers is updated by Equation (11) because the interface displacement is mainly produced by the flow of the upper layer according to the literature [45].
For the interface, is introduced to update the velocity potential expressed explicitly. Therefore, the time-stepping scheme of the velocity potential at the interface can be derived from Equation (3) and is expressed as: In fact, the normal velocity is continuous, but velocity potential of two-layer fluid is not continuous at the interface. Therefore, the interface between the two liquid layers is updated by Equation (11) because the interface displacement is mainly produced by the flow of the upper layer according to the literature [45].
For the interface, ϕ = ϕ 2 − γϕ 1 is introduced to update the velocity potential expressed explicitly. Therefore, the time-stepping scheme of the velocity potential at the interface can be derived from Equation (3) and is expressed as: The interface and their potential distribution are updated using the fourth-order Runge-Kutta method. The initial solitary waves are given based on analytical theory, and the propagation processes are then calculated using fully nonlinear numerical models. For small-amplitude solitary waves, the initial condition is provided by the KdV theory solutions for a two-layer fluid [46]. The numerical simulations are performed using calculation codes developed based on MATLAB, and the flowchart which numerical simulations followed in the present study is shown in Figure 3.
The interface and their potential distribution are updated using the fourth-order Runge-Kutta method. The initial solitary waves are given based on analytical theory, and the propagation processes are then calculated using fully nonlinear numerical models. For small-amplitude solitary waves, the initial condition is provided by the KdV theory solutions for a two-layer fluid [46]. The numerical simulations are performed using calculation codes developed based on MATLAB, and the flowchart which numerical simulations followed in the present study is shown in Figure 3.

Numerical Model Validation
The verification of convergence in grid size and time step was conducted, and a comparison of the interfacial profiles simulated by the multi-domain boundary element method with various grid sizes is shown in Figure 4. It was found that the computational results converge when the grid size is 0.02 x = along the length of the computational

Numerical Model Validation
The verification of convergence in grid size and time step was conducted, and a comparison of the interfacial profiles simulated by the multi-domain boundary element method with various grid sizes is shown in Figure 4. It was found that the computational results converge when the grid size is ∆x = 0.02 along the length of the computational domain and the time step is ∆t = 0.01, grid size of ∆x = 0.02, and time step ∆t = 0.01 are adopted in the present numerical simulations. To evaluate the BEM model, we compared the time history of the displacement of the interface between the two liquid layers obtained from the computations with that of the experiments [47], as shown in Figure 5. Wave profiles are non-dimensionalized using the top layer thickness h 1 and the calculated wave speed c. The computational results were obtained using BEM with a rigid lid boundary condition at the free surface, and the computational results of the BEM are in good agreement with the experimental data. adopted in the present numerical simulations. To evaluate the BEM model, we compared the time history of the displacement of the interface between the two liquid layers obtained from the computations with that of the experiments [47], as shown in Figure 5. Wave profiles are non-dimensionalized using the top layer thickness 1 h and the calculated wave speed c . The computational results were obtained using BEM with a rigid lid boundary condition at the free surface, and the computational results of the BEM are in good agreement with the experimental data.   In order to further verify the accuracy of the numerical model, the fluid velocity in the inner region was calculated using the present method. The velocity potential in the inner region can be expressed as integral Equation (13).
The fluid velocity in the inner region can be expressed as Equation (14) by differentiating Equation (13).
where, Q(r, q) = r 2πr 2 , Q n (r, q) = 1 2πr 2 n − 2r r 2 (r · n) . In order to further verify the accuracy of the numerical model, the fluid velocity in the inner region was calculated using the present method. The velocity potential in the inner region can be expressed as integral Equation (13).
The fluid velocity in the inner region can be expressed as Equation (14) by differentiating Equation (13).
where,  (14) using boundary element method, and the results agree well with the data in literature [48], as shown in Figure 6. The velocity of fluid in two layers region is obtained by solving Equation (14) using boundary element method, and the results agree well with the data in literature [48], as shown in Figure 6. To investigate the effect of a slope on the evolution of ISWs, a series of numerical simulation cases for different slope bottoms and wave parameters were conducted, as shown in Table 1   c 0 is the long wave speed defined by c 2 0 = gh 1 h 2 (ρ 2 −ρ 1 ) To investigate the effect of a slope on the evolution of ISWs, a series of numerical simulation cases for different slope bottoms and wave parameters were conducted, as shown in Table 1. The densities of two layer fluids are set as ρ 1 = 856, ρ 2 = 996, and the depths are set as h 1 = 0.05, h 2 = 0.25.

Wave Evolution
The evolution of internal solitary waves was conducted, and the validation of amplitude variation while ISWs propagating over slope bottom in regard to mesh sizes is exhibited in Figure 7. It can be found that the convergence of present method to calculate the wave evolution can be ensured for the mesh scale ∆x = 0.02. The wave amplitude variations for different slope heights and initial amplitudes when slope angle β = 20 • are plotted in Figure 8. In the first stage (horizontal location x from −1 m to 0), the amplitude is nearly constant when the ISWs propagates above the flat bottom before slope topography, which is the natural attribute of solitary waves. In the second stage (from the slope start point x = 0 to the slope shoulder), the amplitude increases sharply for the up-slope topography and reaches a maximum at the shoulder of the slope, which can explain why the ISWs are observed around the offshore and coastal areas frequently. In the third stage (after the slope shoulder), the amplitude decreases significantly when the ISWs passes the shoulder and propagates over the plateau. The maximum amplitude increases with an increase in the slope height and can reach a nearly 20% increase from the original amplitude for the slope height h s = 0.18 m and initial dimensionless amplitude a 0 /h 1 = 0.4. The amplitude decreases more sharply and reaches a smaller value after the location nearly x =1.25 m (Figure 4b) for a larger slope height owing to the energy reassignment while shoaling along the plateau because the wave suffers more energy rearrangement on the large height slope bottom. slope start point x = 0 to the slope shoulder), the amplitude increases sharply for the up-slope topography and reaches a maximum at the shoulder of the slope, which can explain why the ISWs are observed around the offshore and coastal areas frequently. In the third stage (after the slope shoulder), the amplitude decreases significantly when the ISWs passes the shoulder and propagates over the plateau. The maximum amplitude increases with an increase in the slope height and can reach a nearly 20% increase from the original amplitude for the slope height hs = 0.18m and initial dimensionless amplitude a0/h1 = 0.4. The amplitude decreases more sharply and reaches a smaller value after the location nearly x =1.25 m (Figure 4b) for a larger slope height owing to the energy reassignment while shoaling along the plateau because the wave suffers more energy rearrangement on the large height slope bottom.     The wave amplitude variations at different initial amplitudes are shown in Figure 9. It can be seen from the figure that the dimensionless maximum amplitude with respect to a0 increases with an increase in the initial amplitude, and a more notable decrease can be The wave amplitude variations at different initial amplitudes are shown in Figure 9. It can be seen from the figure that the dimensionless maximum amplitude with respect to a 0 increases with an increase in the initial amplitude, and a more notable decrease can be found for large wave heights above the slope bottom. Meanwhile, internal solitary waves with different initial amplitudes reach the maximum displacement at a similar location (the slope shoulder nearly horizontal location x = 0.5 m). Figure 10 shows the change of internal solitary wave amplitude under different slope angles (the slope angle β = 5 • , 10 • , 15 • , 20 • , respectively). The maximum dimensionless amplitudes increase with the increasing of slope angles, and the attenuation of the internal solitary wave is more obvious for a steeper slope. found for large wave heights above the slope bottom. Meanwhile, internal solitary waves with different initial amplitudes reach the maximum displacement at a similar location (the slope shoulder nearly horizontal location x = 0.5 m). Figure 10 shows the change of internal solitary wave amplitude under different slope angles (the slope angle β = 5°, 10°, 15°, 20°, respectively). The maximum dimensionless amplitudes increase with the increasing of slope angles, and the attenuation of the internal solitary wave is more obvious for a steeper slope.    The variation of ISWs amplitudes versus different densities jump of two-layer fluids is shown in Figure 11, and it can be found that ISWs reach a larger value of amplitude for the larger densities jump. Meanwhile, the degree of attenuation is not affected by densities jump obviously, which is different with the effect of slope parameters and initial amplitude. Figure 12 shows the variation in the dimensionless maximum wave amplitude a m /a 0 for different slope heights, which indicates that the amplitude maximum increases with an increase in the slope height owing to the strong interaction between the wave and the bottom of the slope with a larger height. The increase in the maximum amplitude becomes sharper for the higher slope, which also indicates the apparent effect of slope topography on the evolution of internal solitary waves.
The variation of ISWs amplitudes versus different densities jump of two-layer fluids is shown in Figure 11, and it can be found that ISWs reach a larger value of amplitude for the larger densities jump. Meanwhile, the degree of attenuation is not affected by densities jump obviously, which is different with the effect of slope parameters and initial amplitude. Figure 12 shows the variation in the dimensionless maximum wave amplitude am/a0 for different slope heights, which indicates that the amplitude maximum increases with an increase in the slope height owing to the strong interaction between the wave and the bottom of the slope with a larger height. The increase in the maximum amplitude becomes sharper for the higher slope, which also indicates the apparent effect of slope topography on the evolution of internal solitary waves.

Tail Waves
In addition to the change in the wave height, the waveform of the internal solitary wave changes when it propagates on a slope terrain at different moments at t = 2 s, 4 s, 5 s, 7 s, as shown in Figure 13. In the process of climbing the slope, the lower-layer water continues to become shallower, which decreases the particle velocity of the leading face and results in a steeper waveform. When it passes the top of the slope, the leading edge of the internal solitary wave accelerated, causing the front face to broaden.

Tail Waves
In addition to the change in the wave height, the waveform of the internal solitary wave changes when it propagates on a slope terrain at different moments at t = 2 s, 4 s, 5 s, 7 s, as shown in Figure 13. In the process of climbing the slope, the lower-layer water continues to become shallower, which decreases the particle velocity of the leading face and results in a steeper waveform. When it passes the top of the slope, the leading edge of the internal solitary wave accelerated, causing the front face to broaden. When internal solitary waves propagate on the slope, in addition to changes in waveform and amplitude, after passing the shoulder of the slope, the front face broadened, whereas the trailing face remained steep and the periodic trailing waves of the initial depression wave formed. This is because the lower water depth of the stratified fluid becomes shallower when passing through the slope terrain, resulting in the redistribu- When internal solitary waves propagate on the slope, in addition to changes in waveform and amplitude, after passing the shoulder of the slope, the front face broadened, whereas the trailing face remained steep and the periodic trailing waves of the initial depression wave formed. This is because the lower water depth of the stratified fluid becomes shallower when passing through the slope terrain, resulting in the redistribution of energy to produce periodic tail waves. Therefore, the original wave is redistributed to a new prominent wave following a trailing wave.
The morphological characteristics of the tail waves' time history profile at the observation point (horizontal location x = 1 m) under different slope height conditions are shown in Figure 14. Meanwhile, the moment when wave trough reaches the observation point is set to zero for the intuitional comparison. It was found that the higher the slope height, the more apparent the tail wave, which has a larger wave amplitude and a higher frequency. The maximum displacement of tail waves can reach 50% of ISWs amplitude for the case h s = 0.18 and a 0 /h 1 = 0.4. The characteristics of the tail wave under different initial wave amplitude conditions at the horizontal location (x = 1 m) for the slope height h s = 0.18 and slope angle β = 20 • are shown in Figure 15. It was found that the characteristics of the tail waves are also closely related to the value of the initial wave amplitude. The larger the initial wave amplitude, the more apparent the tail wave, which has a larger wave amplitude. However, the tail waves generated under different initial amplitude conditions had the same frequency. Meanwhile, it was also found that the larger the initial wave amplitude, the more apparent the asymmetry of the evolved waveform, and the smoother the front face.

Conclusions
In this study, we developed a numerical model to calculate the propagation and evolution of internal solitary waves on slope topography, and a series of calculations were performed for different slope geometrical parameters, initial amplitudes, and densities jump of two-layer fluid. The influence of different slope heights and slope angles on the evolution characteristics of internal solitary waves was discussed, and the following conclusions were obtained.

Conclusions
In this study, we developed a numerical model to calculate the propagation and evolution of internal solitary waves on slope topography, and a series of calculations were performed for different slope geometrical parameters, initial amplitudes, and densities jump of two-layer fluid. The influence of different slope heights and slope angles on the evolution characteristics of internal solitary waves was discussed, and the following conclusions were obtained.
In the process of internal solitary waves propagating over slope topography, the dimensionless wave amplitude first increased from the slope start point and reached its maximum value at the shoulder of the slope. Then, the amplitude decreased while propagating over the plateau. The maximum amplitude increased with the height and angle of the slope, and the attenuation of amplitude was more apparent after passing the steeper slope. Meanwhile, the larger initial amplitude and densities jump can induce larger wave amplitude above the shoulder of the slope topography. In the process of climbing the slope, the lower water continues to become shallower, which decreases the particle velocity of the leading face and results in a steeper waveform.
When it passes the top of the slope, the leading edge of the internal solitary wave accelerated, causing the front face to broaden. After passing the slope, the redistribution of energy promoted the generation of periodic tail waves, and the original wave is redistributed to a new prominent wave following a trailing wave. The greater the initial amplitude of the internal solitary wave, the more apparent the tail wave produced. Meanwhile, when the slope height is larger, the tail wave produced also had a larger amplitude and higher frequency. The period of the tail wave is unrelated to the size of the initial wave amplitude; however, it has a closer relationship with the slope height and angle.