Numerical Solution of Ostrovsky Equation over Variable Topography Passes through Critical Point Using Pseudospectral Method

Internal solitary waves have been documented in several parts of the world. This paper intends to look at the effects of the variable topography and rotation on the evolution of the internal waves of depression. Here, the wave is considered to be propagating in a two-layer ﬂuid system with the background topography is assumed to be rapidly and slowly varying. Therefore, the appropriate mathematical model to describe this situation is the variable-coefﬁcient Ostrovsky equation. In particular, the study is interested in the transition of the internal solitary wave of depression when there is a polarity change under the inﬂuence of background rotation. The numerical results using the Pseudospectral method show that, over time, the internal solitary wave of elevation transforms into the internal solitary wave of depression as it propagates down a decreasing slope and changes its polarity. However, if the background rotation is considered, the internal solitary waves decompose and form a wave packet and its envelope amplitude decreases slowly due to the decreasing bottom surface. The numerical solutions show that the combination effect of variable topography and rotation when passing through the critical point affected the features and speed of the travelling solitary waves.


Introduction
The classical model for weakly nonlinear waves propagating over a uniform topography, h, yields to the formation of the Korteweg-de Vries (KdV) equation. The effect of the variable topography, however, must be taken into consideration when retrieving the mathematical model. In consequence, the KdV is succeeded by the variable-coefficient Korteweg-de Vries (vKdV) equation. Johnson [1] was the first person who derive the vKdV equation for surface waves, in which Q = c, and were discussed recently by Grimshaw et al. [2] for internal waves. In the form of the vKdV equation, Grimshaw [3,4] carried out a detailed analysis and an effective model of solitary wave propagation over the variable topography. On the assumption that the flow is two-dimensional, with x and z indicating horizontal and vertical coordinates, respectively, the end result is where A(x, t) is the amplitude of the modal function φ(z), defined by {ρ 0 (c − u 0 ) 2 φ z } z + ρ 0 N 2 φ = 0, for − h < z < 0, φ = 0, at z = −h, (c − u 0 ) 2 φ z = gφ, at z = 0, that can also be utilized to find the linear phase speed, c. Here, g is gravitational acceleration, ρ 0 (z) is the density of the background that can be defined by the buoyancy frequency N (z) in which, ρ 0 N 2 = −gρ 0z and u 0 (z) represent the background current, while t is time variables. The coefficient µ in equation (1) represents the nonlinear term, while λ is the coefficient of 826 Numerical Solution of Ostrovsky Equation over Variable Topography Passes through Critical Point Using Pseudospectral Method the dispersive term. These coefficients µ and λ are defined by the functionality of the different physical structure where and I is written as which is obtained by the derivation from the basic equations (see [5] for details). Some theory regarding the solitary waves over the variable topography on the water surface as well as internal solitary waves are discussed [2,6,7]. According to [8,9], the KdV equation is actually augmented with the additional effects known as a background rotation and lead to the formation of the Ostrovsky equation. Several numerical experiments with these equations have shown that the rotational effect inhibits the evolution of internal solitary waves and forms an envelope nonlinear wave packet [8,10]. An extension of the KdV equation which is the Ostrovsky equation (refer Ostrovsky [11] and Grimshaw et al. [12]), is provided by where the terms inside the bracket are correspond to the KdV equation, while the coefficient γ on the right-hand side represent the rotational effect determined by where f is the Coriolis parameter and The linear dispersion relation of (4) of wavenumber k and frequency ω for sinusoidal waves, sin(kX −ωt), is ω = γ k −λk 3 . Thus, the phase velocity, c and the group velocity, c g are written as According to the dispersion relation, the KdV case (γ = 0) can support the solitary waves solution for all c > 0. Meanwhile, the solitary waves solution could not be supported for Ostrovsky equation (4) due to the rotational effect [13]. The situation with λγ > 0, and background current, u 0 = 0, that is typical for surface and internal waves, results in unsteady envelope solitary wave solutions [13]. Meanwhile, the circumstances with λγ < 0 which is also described as anomalous case in sufficiently strong shear has been lead to a formation of steady envelope solitary wave [14,15].
Nevertheless, the combination of the variable topography and rotation resulted in the development of the variablecoefficient Ostrovsky equation [16,17,18,19]. These effects is believed can led significant impacts on the propagation of the solitary wave. However, the internal waves usually have a negative polarity which has been proved by the observation of internal wave packets in the deep ocean [20]. This is due to changes in the sign of the coefficient in the quadratic nonlinear term, ν, which is well explained by the KdV theory, [2,4]. A cubic nonlinear term is sometimes applied to the KdV equation [2] and to Ostrovsky equation [21,22] in context of the internal wave with a large amplitude, particularly, when the nonlinearity term, µ = 0 is zero. Then the Ostrovsky equation is extended to Gardner Ostrovsky equation The coefficient µ 1 can be found in [23]. In this paper, such cubic nonlinearity is neglected. To conclude, in order to comprehend the Gardner-Ostrovsky equation, one could start by examining the case where the cubic nonlinear term, µ 1 is zero and the case when µ 1 = 0 will be a topic for future study.
The main interest of this study is to investigate the evolution of the internal solitary wave over variable topography when the nonlinearity term slowly and rapidly passes the critical point in the presence of the background rotation, which is an extension from the studies by [17,24]. Similarly, the interest here is still the circumstances where λγ > 0, which is supposed to form an unsteady wave packet. Under the existence of variable topography as the fluid depth, h, is slowly changing, the Ostrovsky equation (4) is retrieved by

Problem Formulation
Equation (9) is acknowledged as the variable-coefficient Ostrovsky equation with the additional components corresponding to the slow variability in the basic state. µ, λ, γ, c, and Q are now in function of x as the modal equation depends parametrically on x, where c(x) is the linear long wave speed while Q(x) is the linear magnification factor, so that QA 2 becomes the wave action flux for linear long waves given by Through the transformation equation (9) is now written as in which Both coefficients ν(τ ) and δ(τ ) are now functions that are defined by the set of the basic state of the fluid, h = h(τ ). Normally, we indicate A(x, t) = A(X, τ ) with h(x) = h(τ ) in which the depth varies gradually along the propagation direction of x. Like the KdV equation, (12) (when σ = 0 and the coefficient is constant) has the initial solitary wave solutions [25,26] that can be written as where the speed, c is proportional to the wave amplitude, a and to the square of the wavenumber, k 2 . In defining the internal wave, two-layer fluids are a typical model that can be used. Consider the density in the upper layer, ρ 1 and lower layer, ρ 2 are fixed for interfacial waves in a two-layer fluid, say ρ 1 = 0.001 and ρ 2 = 0.002. The top layer height is H 1 , and H 2 = H − H 1 is the lower layer depth, where H is the total depth of fluid (see Figure 1). For typical internal waves, the free boundary is replaced with a rigid boundary, and hence, φ(0) ≈ 0 will thus be the upper boundary condition for φ(z).
Assuming u 0 (z) = 0 (no background flow) and setting ρ 0 = 1, then the expressions are attained Substitute (15) into (2), (3) and (5), the coefficients µ, λ, and γ are now where g = g(ρ 2 − ρ 1 )/ρ 1 is reduced gravity. Following [19], the numerical solutions of the equation (12) are obtained on the assumption that the initial parameters for H 1 and H 2 are 1.5 and 1.0, respectively, with f = 1. Noted that, the rotational effect reduced the initial wave amplitude in every case. Therefore, to make it more visible, the initial amplitude, a 0 = 3 is used for all cases in this study. Hence, a typical plot of linear dispersion relations for phase velocity (6) and group velocity (7) of the proposed model (12), where the coefficients described by (16) are constants is shown in Figure 2. According to the dispersion relation in the Figure 2, similarly, the gap in c > 0 for σ = 0 shows that solitary waves can exist with the travelled velocity being positive for all wavenumbers, k. On the other hand, no solitary waves exist for all c when σ = 0. The solitary waves are disintegrated here, and the wave packets emerges with group velocity, c g is negative for every k. Here, the local maximum group velocity of the problem is c gm = −1.5811 at k m = 0.7953.
In most physical problems, the coefficient λ is remains positive for any X. However, an interesting condition occurs when the nonlinearity term, µ, passes through a critical point and changes sign at a particular location where H 2 > H 1 in the deeper water, resulting in an internal wave of depression. It is believed that, in a non-rotating variable medium, σ = 0, the formation of a solitary wave train is categorized into two situations, according to whether µ change its sign or otherwise. These situations have been discussed in detail by Grimshaw [27]. Similarly, the depth of the upper layer here is considered constant for all X and the depth of the lower layer varies monotonically from h 0 to h 1 (refer to Figure 1) in the interval X 0 ≤ X ≤ X 1 .

The Pseudospectral Method
The numerical approach to solve equation (12) is Pseudospectral (PS) method. The numerical solution of the equations is captivating and delicate task, and because of this, only a few numerical studies for Ostrovsky equations have been conducted. This method has been selected among many researches in solving the nonlinear equation especially the integral part when the rotation is added. PS method can considerably speed up the calculation when using Fast Fourier Transform (FFT) which is proven to be a very effective algorithm for computing the Discrete Fourier Transform (DFT). The PS approach uses the Fourier transform to transform spatial derivatives of the PDEs and substitutes the temporal derivative by finite-difference approximation, resulting in a 3-level scheme that needs to be numerically solved [7]. Many numerical codes were generated to solve the constant Ostrovsky equation, see 828 Numerical Solution of Ostrovsky Equation over Variable Topography Passes through Critical Point Using Pseudospectral Method [9,13,15]. The equation (12) is integrated in space τ by the leapfrog finite-difference framework in the temporal time X. Then, the infinite interval is supplemented by −L < X < L with a substantial L to preserve the assumptions of periodicity. Initially, by implementing ξ = sX + π, the solution interval It is now convenient to use W (ξ, τ ) = 1 2 sU 2 notation for the nonlinear terms. Therefore, nonlinear term in (17) is simplified to To obtain the numerical solution of (18), [0, 2π] interval is then discretised by N + 1 equidistant points. Let ξ 0 = 0, ξ 1 , ξ 2 , · · · , ξ N = 2π, so that ∆ξ = 2π N . In the following, we assume N is power of two. By letting m = N 2 , the DFT of where p = −m, −m + 1, −m + 2, · · · , m − 1 while i = √ −1, the typical imaginary number. The inverse Fourier transform of where j = 0, 1, 2, · · · , N − 1, while F ( U ) and F −1 ( U ) are DFT and inverse Fourier transform respectively. The given transformations (19) and its inverse, (20) can be accomplished efficiently using the FFT algorithm, in which greatly minimizes the total number of calculations needed. The derivatives of U with respect to X can be determined using Then, the DFT with respect to X of (18) provides The hat represents the Fourier transform. The approximations listed below can be used Then, substitute (23) into (22), the forward scheme for (12) is attained in the following format .
(24) Equation (24) is a three-level scheme, where to get the third level, U k+1 , one needs to identify the first level, which is the initial condition that shall subsequently refer it as U k−1 and subsequent second level, U k . The step is redo until the required U k+1 is achieved. In obtaining the second level, U k , the interval between U k−1 and U k is divided into ten sub intervals. Hence, ∆τ in (24) is substituted by ∆τ 10 to obtain the equation for U k as .
(25) Since the interval between U k−1 and U k is divided by ten sub intervals, (25) is evaluated ten times to get U k . The values of L and N , as well as the required time step ∆τ is determined accordingly. The numerical resolution for this problem is τ = 0.001 and N = 13684 on the large time scale to avoid any possible issues. Note that the current simulations are just for the research framework. Table 1 shows the comparison of the solitary wave solution obtained by using (14) and the PS method where the effect of the variable topography is omitted.

Results and Discussions
The effect of nonlinearity term on the formation of wave in the presence and absence of background rotation is shown in the next Subsection 4.1 and 4.2. The initial KdV solitary wave solution in (14) is used as initial condition: Generally, a common initial condition for the KdV equation will form a train of solitary waves, but under the background rotation, this wave is diminished and supplanted by an envelope wave packet. According to [9,13], the initial solitary wave decays due to rotation by forming inertia-gravity waves and is absolutely emitted on a time scale of Based on (27), the extinction time for our problem here is τ e = 0.4821, hence, by τ ≈ 1, the solitary waves is said to be completely extinguished and formed a wave packet. Under the influence of variable topography, it should be noted that the coefficient ν in (16) will change its sign if H 2 > H 1 . The numerical results are categorized into two, i.e., one that the nonlinearity term changed its sign in Subsection 4.2 and one that the nonlinearity term did not change its sign in Subsection 4.1. The case when there are polarity changes is separated into two cases where the depth varies slowly and rapidly.

Transformation of Solitary Wave Without Passes Through Critical Point
The transformation of the solitary wave when the nonlinearity term, ν did not change it sign where, H 2 < H 1 after the slope is first examined. The cross-section of long time evolution of formation of waves propagating over gradually increasing depth, in the absence and presence of rotational effect is shown in Figure 3. The depth profile of lower layer is given by As can bee seen clearly from the Figure 3, the rotation effect, σ removed the spectral gap where the solitary wave solution can be occurred and formed the unsteady envelope wave packet. The method used is parallel to the result obtained by [13] which has a good agreement with the theoretical solution and laboratory experiment. The propagation of a smallamplitude trailing shelf behind the solitary wave as it propagates over the deeper region can be observed in Figure 3 (Left), which is in excellent agreement with the result found by [6,19]. Here, the number of waves within the envelope is inversely proportional to the wave amplitude where the amplitude of the waves is determined by the greatest amplitude in the envelope [13]. On a much larger temporal scale, the trailing shelf decays into another trailing radiating wave and another wave packets is emerges when the background rotation exists, for example, see Figure 3 (Right) at τ = 3400.
However, if the depth continually increases, the amplitude of the waves also decreases, and the waves transform into waves of depression when ν passes the critical point (where H 1 = H 2 ). The wave solutions actually depend on how rapidly ν changes its sign [3]. It is believed that the solitary wave is transformed into a radiating wavetrain if ν passes through zero rapidly. If ν changes slowly, the amplitude of the opposite polarity of trailing develops indefinitely as the solitary wave amplitude reduces. Consider the fact that the number of soliton N is obtained by the ratio Z = ν + δ − /ν − δ + , where if Z > 0, the polarity does not change, and N = 1 + [( √ 8Z + 1 − 1)/2]. However, when Z < 0, there is polarity change and no soliton is formed, indicating that the entire solitary wave disintegrates into radiation, according to [6].

Transformation of Solitary Wave Passes Through Critical Point
This section examined the propagation of waves passing the critical points where a dramatic change in the formation of waves is expected to occur. Normally, the solitary wave decays adiabatically and begins to disperse once it propagates over the decreasing slope. However, according to (16), observed that, as τ climbs, ν in (12) shifts from positive to negative values when H 2 > H 1 . The depth profile of the lower layer is taken by The adiabatic behaviour breaks when ν approaches the critical point. When the polarity of the internal solitary wave changes its sign, the leading solitary wave deforms, its amplitude gradually decreases, and an elevation rarefaction wave is produced [24,28]. At this time, its amplitude decreases while the trailing shelf formed grows in amplitude, see Figure 4 (Left). Based on (29), it is not until τ > 261 when ν start to change its sign and the initial internal solitary wave of the negative trailing shelf (in H 2 < H 1 region) becomes an initial disturbance to generate an undular bore, which can also be regarded as a secondary solitary wave [29,30]. This formation, also known as a dispersive shock wave, is formed as a result of the steepening caused by the nonlinear effect, which triggered the dispersive effects. This pattern has been observed on the free-surface in Northern Australia and the Andaman Sea, for example, with the leading wave followed by a train of welldeveloped undulations where the leading solitary wave travels at a different velocity than the trailing edge. It happens once the fronts of the nonlinear long waves become very steep. The existence of this phenomenon is due to the reaction between travelling solitary waves and bottom topography, hence being easily disturbed by the change of polarity caused by topography conditions. This undular bore transforms into a linear wave on a large time-scale, due to the diminishing pedestal. Hence, the amplitude of the waves decreases over time. 831 a packet of waves followed by trailing waves. The effect of rotation on internal solitary waves has been observed in the South China Sea [16]. The bigger value of Coriolis frequency or rotational effect chosen, however, has totally extinguished the initial solitary wave and formed a wave packet, a localized disruption caused by the accumulation of numerous different wave forms. The variation in the topography then affects the features and speed of the travelling wave packet. Clearly, the wave packet disperses while propagating due to variable topography. The packet structure with the largest wave amplitude has nonlinear features and satisfied the theoretical solution. Within all the emitted radiation, the wave packet formed has the slowest absolute speed [8,9]. In this case, the effect of rotation tends to be dominant, where the nonlinearity effect is generally suppressed, as the variance of the nonlinear coefficient in the previous case and in this case does not seem to have a significant influence on the subsequent evolution except for the amplitude and speed. It follows that the smaller amplitude of the leading wave and trailing edge of the wave trains lead to the formation of the smaller amplitude of the wave packets and its secondary trailing wave packet, and hence, travelled faster. The numerically obtained speed of nonlinear wave packets (in absolute value) propagating over gradually increasing depth when the nonlinearity term change, (29) and did not change its sign, (28) is given in Table 2 and the comparison is shown in Figure 5. Next, the wave generation over a rapidly changing slope is examined. The interest is still when there is polarity change, which is when ν changed its sign rapidly from positive to negative. Both depth profile conditions for H 2 satisfied h 1 > H 1 .
As previously discussed, when the nonlinearity term, ν, rapidly passes through the critical stage, the solitary waves are absorbed and completely transformed into a radiating wave train [3]. Here, the whole structure reduced its amplitude into a radiating wave train. The results obtained shown in Figure 6 (Left) are in agreement with previous work. Although others have been looked at this problem intensively, here we were able to examine further the transformation of a solitary wave at the turning point whenever the background rotation is taken into consideration. Similarly, the amplitudes of the wave trains, as well as the amplitudes of the wave packets and their trailing radiating waves, decrease over time on a large time scale, as shown in Figure 6 (Right).
These two cases from (29) and (30) involve polarity change in (12) as the waves propagates over the rapidly and slowly increasing depth. When H 2 > H 1 , a train of solitary waves of negative polarity is formed. Here, as σ = 0, the leading solitary wave has identical amplitude with the trailing shelf close to the critical point (ν = 0), and they are usually smaller than their counterparts without a change in polarity. Likewise, in the case of σ = 0, the combined influence of the nonlinearity and rotation is applied over the solitary wave, and one of the predicted formation outcomes is that if ν is positive (τ ≤ 100), the result is distinguished by a greater amplitude with a few number of waves in the envelope, as opposed to that in the case of ν is negative (at τ = 1000, 3400). When the polarity changes, the leading solitary wave train cannot maintain its shape while passing through the critical point and is slowly defeated by the evolving rarefaction waves but the result is distorted due to the combined rotational effect; see the right Figures 4 and 6.
Typically, in this case, the internal waves decay into several wave packets followed by a few residual waves. Any nearly localized packets of waves eventually evolve, each containing a long-wave envelope that is shorter and travelled faster. After a long generation, a different number of wave packets emerge and the formed wave packet becomes more complicated with the presence of more wave trains. The formation of a nonlinear wave packet when σ = 0 depends on the wavetrain formed when σ = 0, as shown in Figures 4 and 6. If the amplitude of the wavetrain is small, then the amplitude of the nonlinear wave packet is also reduced [13]. As a result, the nonlinear wave packet formed will be narrower and the number of the waves inside the envelope wave packet will increase. Furthermore, noticed that the extinction time, (27) for the proposed model is τ e = 0.4821, which is too short compared to the total run time in our numerical solution, τ ≈ 3400 to take the variation of the nonlinear term into account, resulting in the dominant rotational effect.

Conclusions
In this study, the formation of an internal KdV-type solitary wave of negative polarity over a rapidly and slowly varying topography in a two-layer fluid system in the presence of background rotation are simulated. The numerical result indicates that the polarity of the internal wave may vary following the depth of the region. If the depth after the slope is greater than the depth of the upper layer, where there is a polarity change, the solitary wave is converted into a train of solitary waves of negative polarity and is accompanied by a radiation wave. At a larger timescale, the amplitude of the solitary waves decreases and diminishes due to the pedestal [24].
However, if there is a background rotation, the wave is gradually replaced by a wave packet accompanied by another secondary trailing wave packet. The existence of secondary trailing wave packet become more obvious after long time evolution. In this situation, the rotational effect seems dominant, with the nonlinear effect indicating some small influence on the evolution of the wave packet which has been supported by [17]. As described, the presented findings could be due to much longer model run-time, as the model runs much longer than the extinction time in (27). Unlike in the KdV case, where nonlinearity has a tremendous impact on the formation of solitary waves, this structure leads to some disruption when the background rotation is applied. As a conclusion, for the case when σ = 0, in both cases (when ν > 0 and ν < 0), the re-sulting speed of the wave packet shows that the wave packets travel faster if the depth is deeper (see Figure 5).
However, the amplitude of the wave packet is still reduced from the initial amplitude, a 0 , due to the rotational effect that makes the energy of the soliton slowly decrease in both cases. In contradiction, the KdV soliton keeps maintaining its amplitude once it propagates into another constant depth when ν > 0, while the leading waves continuously reduce its amplitude when ν < 0. Besides, the wave destroys itself when it attempts to pass through the critical point which has been studied numerically by [28]. Nonetheless, the results show that they are consistent with observational studies, that the rotation and nonlinearity terms have a significant impact on the propagation of internal waves [31]. However, these conclusions, of course, must be revised whenever the cubic nonlinear term, µ 1 is considered.
As stated earlier, a cubic nonlinear term is frequently applied to the equations, particularly when there is a polarity change. The implications of a solitary wave have been studied thoroughly, which depend on the sign of the cubic term [2,4]. But then, in such an extended KdV equation, a similar theory for the internal waves still has to be developed and therefore would tend to be even more complex. Lastly, it should note that the study has decided to consider the variable-coefficient Ostrovsky equation rather than incorporating other effects to compensate for the loss of the quadratic nonlinear term at the turning point. In fact, the existence of a cubic nonlinear term is likely to have a major impact on the process of transformation. The addition of the cubic nonlinear term would almost certainly result in phenomena that are not comparable to those reported here.