Solitary-Wave Solutions of Benjamin–Ono and Other Systems for Internal Waves: II. Dynamics

Considered here are two systems of equations modeling the two-way propagation of long-crested, long-wavelength internal waves along the interface of a two-layer system of fluids in the Benjamin–Ono and the Intermediate Long-Wave regime, respectively. These systems were previously shown to have solitary-wave solutions, decaying to zero algebraically for the Benjamin–Ono system, and exponentially in the Intermediate Long-Wave regime. Several methods to approximate solitary-wave profiles were introduced and analyzed by the authors in Part I of this project. A natural continuation of this previous work, pursued here, is to study the dynamics of the solitary-wave solutions of these systems. This will be done by computational means using a discretization of the periodic initial-value problem. The numerical method used here is a Fourier spectral method for the spatial approximation coupled with a fourth-order, explicit Runge–Kutta time stepping. The resulting, fully discrete scheme is used to study computationally the stability of the solitary waves under small and large perturbations, the collisions of solitary waves, the resolution of initial data into trains of solitary waves, and the formation of dispersive shock waves. Comparisons with related unidirectional models are also undertaken.


Introduction
In the precursor [1] to this work, the authors considered Benjamin-Ono and Intermediate Long Wave systems of the form for x ∈ R and t ≥ 0. These were derived in [2] and are here written in unscaled, dimensionless variables.The system (1) is a one-spatial-dimensional model for the two-way propagation of long-crested internal waves along the interface of a two-layer system of homogeneous, inviscid fluids of densities ρ j , j = 1, 2, with the bottom layer density ρ 2 > ρ 1 for static stability, and depths d j , j = 1, 2. The upper layer is taken to be bounded above by a rigid surface (the so-called rigid lid approximation) and the lower layer is bounded below by a horizontal, flat, impermeable bottom.As mentioned, these models allow for bi-directional propagation, unlike the original Benjamin-Ono and Intermediate Long-Wave equations (abbreviated to BO and ILW henceforth).However, they do assume the waves are long-crested, so there is no appreciable variation in the horizontal direction orthogonal to the principle direction of propagation.The independent variable x is proportional to position in the direction of propagation while t is proportional to elapsed time, ζ = ζ(x, t) is the deviation of the interface from its rest position at the point x at time t.The dependent variable u = u(x, t) is a horizontal velocity-like variable while γ = ρ 1 /ρ 2 < 1 is the density ratio and α is a modeling parameter.
In terms of the non-dimensional parameters (where a and λ denote, respectively, a typical amplitude and wavelength of the interfacial wave and δ = d 1 /d 2 is the depth ratio), the physical regimes under which the Euler equations for internal waves are consistent with the two-dimensional version of ( 1) are (see [2]): • Intermediate Long-Wave regime: µ ∼ ǫ 2 ∼ ǫ 2 ≪ 1, µ 2 ∼ 1. (This means that the upper layer is shallow, and the interfacial deformations are small with respect to the depths of both layers.) • Benjamin-Ono regime: µ ∼ ǫ 2 ≪ 1, ǫ 2 ≪ 1, µ 2 = +∞ (corresponding to a very deep lower layer).
The two regimes are mathematically distinguished in (1) by the definition of the Fourier multiplier operator H.This takes the form H = ∂ x T √ µ 2 , with in the ILW case (where P.V. stands for the Cauchy principal value of the integral), and H = ∂ x H, where is the Hilbert transform, in the BO case.In terms of the corresponding Fourier symbols, if f denotes the Fourier transform of f , then Hf (k) = g(k) f (k), k ∈ R, where Briefly described now are some mathematical properties of the initialvalue problem for (1) (see [1] for details).The system (1) is linearly well posed if and only if α ≥ 1 [2] and the nonlinear initial-value problem for the ILW case has been shown in [3] to be well posed, locally in time, in suitable Sobolev spaces if α > 1.In the same paper, again when α > 1, similar well-posedness results for the BO system are suggested to hold because of the convergence of solutions of the ILW system to those of the BO system as µ 2 → ∞.
On the other hand, to the best of our knowledge, the system (1) for α = 0 has only the simple conservation laws whereas the case α = 0, which is ill posed [4] admits two additional invariant quantities, and a Hamiltonian structure emerges.
Another important property of (1) is the existence of solitary-wave solutions of the form ζ(x, t) = ζ(x − c s t), u(x, t) = u(x − c s t), moving with constant speed c s = 0, having profiles ζ = ζ(X), u = u(X) which, along with their derivatives, decay to zero as X = x − c s t → ±∞.They must be solutions of the system Existence of smooth, small-amplitude solutions of (6) was recently proved by Angulo-Pava and Saut in [5].In the same paper, the decay as |X| → ∞ was shown to be exponential in the ILW case and algebraic (as 1/|X| 2 ) in the BO case, just as for the solitary-wave solutions of their undirectional counterparts.
The numerical generation of solitary-wave solutions of (1) was discussed by the authors in [1], which is Part I of the present study (and see also [6]).Three iterative methods, one of Newton type, the classical Petviashvili iteration, and a modification of it, were proposed to solve iteratively a discretization of (6) based on Fourier collocation approximation of the corresponding periodic problem.The three schemes were compared in accuracy and some properties of the resulting computed solitary waves that emerged from the experiments were pointed out.These concerned the speed-amplitude relationship, the convergence of the solitary waves of the ILW system to those of the BO system, and comparisons with solitary-wave solutions of both the classical and regularized versions of the unidirectional BO and ILW.The regularized versions of these unidirectional equations, which result from using the lowest-order relation, ∂ x = −∂ t to modify the dispersion relation, will be denoted henceforth by rBO and rILW, respectively.
With a reasonably good grasp of the solitary-wave solutions in hand, it is natural to turn to their dynamics as solutions of the time-dependent problem.To this end, the corresponding periodic initial-value problem (IVP) is discretized in space with a spectral method and in time with the explicit, fourth-order Runge-Kutta (RK) scheme.Error estimates of the spectral semi-discretization were proved in [6].These estimates naturally depend upon the regularity of the solutions of (1) and, in particular, lead to spectral convergence in the smooth case.In addition, the fully discrete method was also used in Part I of this project to check the accuracy of the computed solitary-wave profiles.Its very satisfactory performance gave the confidence needed to make use of it for the developments in the present essay.
From this computational perspective, several properties of (1) are studied.The first group of experiments analyzes the stability of the solitarywave solutions.Small perturbations of the traveling-wave profiles, which are determined as in [1], are taken as initial conditions of the time-dependent numerical scheme.The evolution of the corresponding numerical approximation suggests that the solitary waves are asymptotically stable.Indeed, what is observed is that a perturbed initial solitary wave develops into a principal part which appears to converge rapidly to a new solitary wave, together with smaller waves of a purely dispersive nature.There are two groups of these dispersive waves, traveling in opposite directions.
One would expect that internal solitary waves would be unstable or lead to series of solitary waves when perturbed by a large amount.Surprisingly, the experiments suggest that this is not the case for the BO and ILW systems.Large perturbations of their solitary waves appear to develop into a single, large solitary wave.The amplitude of the emergent solitary wave is related to the size of the perturbation in a way to be discussed later.The solitary wave is followed by comparatively small amplitude, dispersive tails quite similar to those appearing when small perturbations are considered.
Another issue studied in the present paper, which appears related to their stability, is the interaction of solitary waves.Since (1) admits twoway propagation, both experiments of head-on and of overtaking collisions are carried out.As the interactions are expected to be inelastic, additional information provided by the experiments will concern the emerging waves and the nature of the dispersion resulting from the collisions.We complete the computational study of properties related to the stability of solitary waves by analyzing numerically the resolution of general initial data into trains of solitary waves along with dispersive tails.In particular, the role of the energy and mass of the initial data, represented in terms of its amplitude and wavelength, is examined in the context of this resolution property.
The computational study of (1) will be concluded with a discussion of two additional issues.The first one concerns the relationship between the system and the corresponding unidirectional equations.More specifically, we are interested in the dynamics of (1) when presented with data corresponding to unidirectional propagation.To this end, we analyze computationally the evolution of the interface according to the system using initial data derived from solitary-wave solutions of the unidirectional model and vice-versa.We find that the systems compare qualitatively, and to some extent quantita-tively, with their unidirectional counterparts in the unidirectional regime.Indeed, in both cases, the numerical experiments reveal that the evolution that emerges is similar to that which obtains from small perturbations of the exact solitary-wave solution of the relevant system.
The second point concerns the formation of dispersive shock waves (DSW henceforth).When they exist, one expects such waves to involve two scales, one fast and oscillatory and a second one, slow and modulational, cf.[7,8].Thus, such waves would consist of a wavetrain with a local, rapidly oscillating structure while the envelope wave parameters themselves are changing on a much slower time scale.In the present paper, the formation of DSWs associated to (1) is investigated by integrating numerically two problems typically related to DSW formation: the Riemann problem and the dam break problem.
The paper is structured as follows: Section 2 is devoted to a description and an accuracy study of the numerical method used for the computational work to follow.In Section 3, the stability of the solitary waves is analyzed via a group of experiments concerning the evolution of small and large perturbations of the solitary-wave initial data.The investigation of head-on and overtaking collision of solitary waves, as well as the resolution property is presented in Section 4. Discussion of some relationships between unidirectional models and the associated bi-directional systems is the subject of Section 5, while Section 6 contains the computational study of DSWs.Some concluding remarks may be found in the final Section 7.
Overall, the dynamics for the ILW and the BO systems appear qualitatively quite similar.To keep the length of the script under control, the present paper will mostly report experiments on the dynamics of the BO system, partly because BO is computationally more challenging owing to the slower decay of its solitary waves.

Description of the numerical method
In this section the fully discrete scheme for the numerical approximation of (1) is described.The particular method was already used in [1] to check the accuracy of the computed solitary-wave profiles.This involved considering the computed solitary-wave profiles as initial conditions and monitoring by how much the resulting solutions evolve away from an appropriate translation of the profile.
While the problems under consideration are set on the spatial domain R, the computations take place in a spatially periodic setting.The commonplace practice of approximating IVP's on R having localized initial data via periodic problems goes back at least to the 1960's.(see e.g.Zabusky and Kruskal [9] and Tappert [10]).Preliminary arguments justifying such approximation for unidirectional models can be found in [11], but see the more recent work of H. Chen [12], where explicit error estimates, valid over long time scales depending on the length of the period, are obtained.The same type of analysis can be applied to the bidirectional case.
Taking this point as settled, attention is focused on the periodic IVP for (1) on a long enough interval [−l, l], with initial conditions given by two smooth periodic functions ζ 0 , u 0 , and where the dispersion operator H takes the form of the corresponding periodic operator, which is to say, H is a multiplier on the Fourier coefficients by the function whose discrete symbol is given as in (4) (cf.[1]).
For an even integer N ≥ 1, let S N be the space of trigonometric polynomials of degree at most N, viz.
Let {x j = −l + jh, j = 0, . . ., 2N − 1} be a uniform mesh of nodes on [−l, l] with h = l/N.Let P = P N denotes the L 2 -projection onto S N while If = I N f denotes the interpolating trigonometric polynomial to f based on the nodes x j .
For T > 0, the semidiscrete Fourier-Galerkin approximation is defined as the mapping for t ≥ 0, with initial data The spectral discretization ( 7) was introduced in [6] and L 2 -error estimates were established there.For the experiments below, and to take advantage of the numerical generation of solitary waves performed in [1], the method was implemented in pseudospectral form.Algebraically, this is equivalent to the collocation formulation [13].With a small abuse of notation, (ζ N , u N ) : [0, T ] → S N × S N will also denote the semidiscrete Fourier collocation approximation, which is to say, the mapping The semidiscrete system ( 8)-( 9) is conveniently formulated using a nodal where •'s connote Hadamard products, D N is the Fourier pseudospectral differentiation matrix of order 2N and H N stands for the corresponding discrete version of the operator H.In terms of the k th discrete Fourier coefficients ζ N (k) and u N (k) of ζ N and u N , the discrete Fourier transform allows us to write (10) in the form where k = πk l , g k = g(| k|), −N ≤ k ≤ N and g is as given in ( 4).The initial-value problem for the semidiscrete system (11) is integrated in time using the explicit, 4th-order Runge-Kutta method.
Remark 2.1.It is worth mentioning some additional properties of the semidiscrete schemes.Note first that for the continuous, periodic initial-value problem for the system (1) on (−l, l), the analogs of ( 5) are also independent of time.A natural discretization of ( 12) is given by hI 1,h , hI 2,h where and 1 N denotes the vector with all components equal to one.The angle brackets •, • here denote the Euclidean inner product in C N .Note that if F N stands for the N × N matrix associated to the discrete Fourier transform and and where H N (respectively D N ) denotes the diagonal matrix with diagonal entries Therefore, taking the inner product with 1 N in the first and second equations of (10) reveals that the quantities ( 13) are preserved up to round-off error by the semidiscrete approximation.

Stability of solitary waves under small and large perturbations
In our study of the stability of solitary waves of both ILW and BO systems, several types of numerical experiments are performed.The general procedure is always the same: first generate an approximate solitary-wave profile by using one of the methods described in [1].Then, this computed profile is perturbed in one or both of its components ζ and u.Finally, the perturbed profile is used as an initial condition for the fully-discrete scheme described in Section 2 and the evolution of the resulting numerical approximation is monitored.
In all cases, the solitary waves appear to be stable under small perturbations.The behavior is illustrated by the following example: Consider the BO system with α = 1.2, γ = 0.8 and ε = √ µ = 0.1 and the solitary wave with c s = 0.57 and amplitude approximately a = 3.9129.This solitary wave is perturbed by multiplying the function ζ 0 (x) by the factor r = 1.1.Thus the perturbed initial data is (rζ 0 , u 0 ) where (ζ 0 , u 0 ) is the unperturbed solitarywave profile.This represents a perturbation of about 10%. Figure 1 shows the evolution of the perturbed solitary wave to a new, stable solitary wave followed by a rather interesting dispersive tail.The dispersive tail consists of a right and a left-traveling part as indicated in Figure 1.In this figure we only show the ζ-component of the solution since u behaves similarly.Analogous behaviour has been observed in the case of Boussinesq systems for surface waves in [14].The accuracy of the results is guaranteed in this case by using N = 32, 768 nodes in the interval [−512, 512] (so h = 3.125 × 10 −2 ) and ∆t = 0.01.The results reported here are representative of a whole series of tests with r = 0, but |r − 1| fairly small.The results obtained for the ILW case are quite similar to those presented in Figure 1.The main difference is that, because the solutions decays ex-ponentially to 0, the simulations are more easily carried out since a smaller spatial interval is required.This stable behaviour seems to persist for larger perturbations.For example, Figure 2 presents the evolution of a solitary wave of the ILW system with an amplitude perturbation r = 1.5 in the ζ component only.Keeping the parameters of the ILW system the same as in the BO system, the unperturbed speed of the solitary wave was taken to be c s = 0.51.With these specifications, the amplitude was approximately a = 1.7.The evolution of the perturbed initial condition leads to a new solitary of amplitude approximately a = 2.35 and two counter-propagating dispersive tails.A series of simulations with larger values of r showed similar behavior for both BO and ILW systems.The emerging solitary wave had roughly the amplitude rA, where A is the amplitude of the unperturbed solitary wave.These results speak to the extreme stability of the solitary-waves solutions of both systems.

Dispersive tails
Observe in Figure 1(d) that the tail generated by the perturbation of a solitary wave consists of two dispersive groups, a left-and a right-propagating part.These seem to be well separated and there appears to be no interaction between them or with the emerging solitary wave.This phenomenon has also been observed in Boussinesq systems for surface waves.It is worth remarking that in the Boussinesq system context, there were cases where an interaction between the two parts of the dispersive tails was observed; this seemed to occur when blowup phenomena was imminent (see [14]).
A little analysis based on observation of the outcome of numerical simulation casts some light on the remarks above.Define two quantities (Note that R < 1.)As seen in [1], the speed of the solitary waves seems to be bounded from below by c γ in case of the BO system and by Rc γ in case of the ILW systems.The solution of these systems that start from a perturbed solitary wave evolves into a solitary wave of speed c s , say, and a small remainder, which we are calling the dispersive tail.We suppose c s > c γ in the BO case and that c s > Rc γ in the ILW situation.Small solutions of  the system (1) satisfy approximately the linearized system where X = x − c s t is a frame of reference moving with a solitary wave of speed c s > c γ in the BO case and c s > Rc γ in the ILW case.
Since the operators H and (15) and the second equation is used to solve for (∂ t − ∂ X )u in terms of ζ, there obtains Solutions of ( 16) are superpositions of plane waves ζ(X, t) = e i(kX−ω(k)t) where .
As before, g is given in (4).Note that g is positive and increasing for x ≥ 0. Define g(0) = 1/ √ µ 2 in the ILW case.If one assumes α ≥ 1, which is required for well posedness, then φ has the following properties: Therefore, for all wavenumbers k > 0, the phase speed of the dispersive tail, relative to the speed of the solitary wave, satisfies The inequalities (19) show in which direction the plane-wave components of the dispersive tails propagate, relative to the solitary wave.Since φ(0) < 1, the absolute phase speed v + (k) + c s of the right-traveling waves e i(kX−ω + (k)t) is less than or equal to c γ , and similarly the left-traveling waves e i(kX−ω − (k)t) have |v − (k) + c s | ≤ c γ .Moreover, components corresponding to longer wavelength (smaller k) are faster than those of shorter wavelength.
Examine now the associated group velocities.Observe that, 2 , x > 0, Note first that (18) implies that ψ(x) ≤ φ(x) for x ≥ 0. On the other hand, according to (4), xg ′ (x) = g(x) + x 2 C(x) where Consequently C(x) ≤ 0 for x > 0 and so xg ′ (x) ≥ g(x) for x > 0. This and the hypothesis α ≥ 1 imply that In consequence, it transpires that 0 ≤ ψ(x) ≤ φ(x), for x ≥ 0, and Therefore, for all wavenumbers k > 0, By using ( 14), (17), and the hypotheses on c s , it is deduced that c s > φ(0)c γ in both the BO and ILW cases.Then (20) means that there are two dispersive groups, one traveling to the left and one to the right (with the solitary wave), but with a group velocity smaller than c s .

Solitary wave interaction and the resolution property
The present section is concerned with the interactions of solitary waves.We remind the reader that for the unidirectional BO and ILW models, overtaking interactions are elastic owing to the complete integrability of the equations.For the regularized versions of these models, apparently the collisions are not elastic as witnessed by the simulations in [15,16].Our purpose now is to study computationally the interactions of solitary waves of the BO and ILW systems (1).Since the models support propagation in both directions, two types of collisions can be considered.

Head-on collisions
The experiments begin with a symmetric head-on collision for the BO system.This is illustrated by the following example.Two computed solitary waves with the same speed c s = 3.2 (with α = 1.2, γ = 0.1, ε = 0.5, √ µ = 0.1) on [−4096, 4096], propagating in opposite directions are computed.(To reverse the direction of propagation, simply remark that if (ζ, u) is a solitary wave propagating to the right with speed c s , then (ζ, −u) is a solitary wave propagating the the left with speed −c s .).The location of the peak amplitudes of the two solitary pulses was initially placed at x = −800 and x = 800, respectively.(In this and all the rest of the experiments of this section, h = 0.125 and ∆t = 0.01.)The head-on collision is illustrated in Figure 3 at several time instances of the evolution of the numerical approximation to the deviation of the interface.The interaction appears to be inelastic and similar to the head-on collision of Boussinesq systems for surface water waves, (see eg. [17]).During the interaction, the solitary waves change slightly in shape and their speed decreases.After the inelastic head-on collision the emerging solitary waves have altered amplitudes and speeds and they shed a symmetric dispersive tail propagating in both directions.(A magnification of the dispersive tail is shown in Figure 3   The peak amplitude of the numerical approximation during the interaction is presented in Figure 4(a), while Figure 4(b) shows the locations of the peak amplitudes.The amplitudes and their locations have been computed using Newton's method as described in [1].The amplitude before the collision was A = −0.362898while after the (symmetric and inelastic) interaction it has decreased to A = −0.362912which is a 0.39% change.A similarly small phase change can be observed.
In addition, we performed asymmetric head-on collisions, with similar results and the same conclusions.This inelastic character is also observed in the case of the ILW system.Due to the exponential decay rate of the solitary wave profiles in this case, the numerical simulations do not require such a large spatial domain for accurate approximation.

Overtaking collisions
We also consider a second type of interaction between solitary waves, namely overtaking collisions.These occur when a wave of one speed is placed behind another with a smaller speed.Recall again the elastic behavior of this type of collision for the unidirectional ILW and BO equations.On the other hand, the regularized versions of these unidirectional models featured inelastic collisions (according to the experiments reported in [15,16] for the BO case).
For the BO system, the typical behavior for this interaction is illustrated by the following experiment.Taking α = 1.2, γ = 0.1 and ε = √ µ = 0.1 on the interval [−4096, 4096], two solitary waves with speeds c s = 3.5 and c s = 3.1 (traveling to the right) with amplitudes A s = −4.7157and A s = −0.8978respectively, have been considered as initial data for the numerical integration of (1).The peaks of these two soliary waves were sited initially at x = −1, 000 and x = 1, 000.The numerical computations used N = 65, 536 modes (corresponding to h = 0.125) and ∆t = 0.01.For comparison, the same experiment for the corresponding solitary-wave solutions (with the same speeds) of the rBO model is also performed.The solitary-wave solutions for rBO have an exact formula, and it is not necessary to generate them numerically.Having the same speeds as those for the BO system, the amplitudes are not the same due to the different speed-amplitude relation, see [1]. Figure 5 presents the ζ profile during the interaction for different propagation times.The experiment suggests an inelastic overtaking collision of the solitary waves for both the BO system and rBO equation.We also note (see especially 5(d)) that the interaction of the solitary waves in   (f): numerical approximation at times t = 0, 4800, 5000, 5500, 6000, 8000.After the interaction, there is a small difference in the speeds of the larger emergent solitary waves.
the case of the rBO equation looks more inelastic (more nonlinear) and lasts longer.Perhaps as a result of this longer interaction time, the solitary waves of the rBO equation after the interaction appear to possess a larger phase change.This more inelastic character in the rBO case is noted, for example, by the fact that after the interaction, the larger solitary wave of the rBO equation has a change in amplitude of about 0.68%, in comparison with the BO system, where the change is about 0.51%.
As expected, after the inelastic interaction, dispersive tails were generated.For the experiments under discussion, these emerging tails are shown in more detail in Figure 6 for the BO system and in Figure 7 for the rBO equation.The magnification given in Figure 6(b) shows that after the interaction of the solitary waves of the BO system, an N-shaped wavelet traveling to the left is generated.This does not appear in the case of the rBO equation, see Figure 7(b).The formation of these types of waves has also been observed in unidirectional and bidirectional surface water-wave models, cf.[18,19].The rest of the dispersive tails are comparable in both size and wavelength for the two models, see Figures 6(d    The dispersive tails generated in the case of the ILW system consist of the same two parts as in the BO system and include the N-shaped wavelet.The collision is thus clearly inelastic, as is the case for the rILW equation.The result of these interactions is illustrated in Figure 8.

Resolution into solitary waves
Another important property, which is related to the stability of the solitary waves, is the process whereby solitary waves emerge from the evolution of more general initial conditions.This is often referred to as resolution into solitary waves, and explains why these special solutions attract so much attention.To illustrate this resolution property in the current context, a series of runs is reported.The initial conditions in this section are Gaussian perturbations of the interface with zero initial velocity, i.e.
where A and λ are constants.Such initial conditions are reminiscent of the early stages of tsunami generation [20].For these initial conditions, the   resolution seems to be determined by the energy and mass of the pulse, which in this case is related to the amplitude A and the wavelength λ, as seen in what follows.
Consider the BO system with α = 1.2, γ = 0.1 and ε = √ µ = 0.1 with initial conditions (21), where A = 2 and λ = 10.These initial conditions are unable to trigger the generation of solitary waves.Instead, two counterpropagating dispersive wavetrains appear as shown in Fig. 9. Indeed, we monitored the L ∞ -norm of the solution as a function of t and observed that it decreased steadily as t increased from 0 to 2, 000.This decay of the maximum norm is one of the hallmarks of a solution containing no solitarywave component.The result is not surprising, as the branch of solitary-wave solutions of the BO system appears to have energies that are bounded away from 0 Increasing either the amplitude or wavelength of the initial condition has the effect of activating the generation of solitary waves.For example, the evolution of ( 21) with A = 5 and λ = 10 is shown in Figure 10.Note that two counter-propagating solitary waves are generated followed by dispersive tails.Similarly, increasing the wavelength of the initial condition by taking λ = 20 and holding A = 5 produces the four solitary waves shown in Figure 11, while Figure 12 shows the resolution of the initial condition with A = 1 and λ = 100 into at least four, and very probably six, solitary waves.

Comparisons with solitary waves of the BO equation
As discussed in Part I of this work [1], the rBO and rILW equations describe the propagation of internal waves in mainly one direction and have solitary-wave solutions found in closed form.The BO and ILW systems should be able to describe waves propagating in one direction similarly to the unidirectional models provided that the initial data are 'unidirectional'.In the present Section, the ability of the two-way systems to model the waves of their unidirectional counterparts, and especially solitary waves of the rBO and rILW equations is investigated.The results for the BO system are shown below; the experiments for the ILW system are easier to perform and the overall conclusions are the same.
The tricky point here is that unidirectional models require initial data only for the interfacial variable ζ(x, 0) = ζ 0 (x), whereas the two-way propagation systems also require an initial velocity u(x, 0) = u 0 (x).Guided by the rigorous theory of surface-wave Boussinesq systems [21], we use the higherorder approximation to the unknown velocity that applies in the unidirectional case (see Equation (11) of [1]).Thus if we consider a solitary wave ζ 0 (x) of the rBO equations, given exactly in Equations ( 13)-( 16) of [1], then the unknown initial velocity u 0 (x) for the BO system is consistently chosen by the higher-order approximation (22).
For the first reported experiment, we took α = 1.2, γ = 0.8, ε = √ µ = 0.1.Figure 13 presents the evolution of the BO system with initial data for ζ 0 taken to be a solitary-wave solution of the rBO equation with c s = 0.53 and amplitude A = −1.6,together with the u 0 derived from (22).The numerical approximation evolves into a principal profile of solitary-wave type, along with a trailing tail.Most of the resulting wave travels to the right with the rBO-solitary wave, while a very small portion travels in the opposite direction.After a simulation which is long enough for the tail to separate and disperse, the numerical approximation appears to be a new solitary wave of amplitude A ≈ −1.78, differing from that of the initial profile by about 11.25% in the maximum norm.What is particularly notable here is that this 'unidirectional' data does indeed produce a signal that, in the main, propagates in only one direction.22) for the BO system.The maximum negative excursion of the solution as a function of time is also displayed.Observe that the negative excursion levels off as t grows, consistent with a traveling-wave structure.
Similarly, a simulation was run of the rBO equation with initial condition given by a solitary-wave solution of the BO system.The solitary-wave profile was generated following the numerical techniques introduced in [1] and with speed c s = 0.53.The computed speed-amplitude relation gives a numerical profile of the amplitude A ≈ −1.63.The IVP for the rBO equation is approximated by a fully discrete scheme based on the same procedures described in Section 2 for the BO system.The evolution of the corresponding numerical approximation is shown in Figure 14.As in the previous experiment, this shows a solution of the rBO equation consisting of a solitary-wave profile followed by a tail.As time evolves and the tail disperses, the approximation tends to an emerging solitary wave with amplitude A ≈ −1.47, which is a difference in uniform norm of about 10% when compared with the amplitude of the initial profile.

Dispersive shock waves
This section is devoted to studying numerically the formation of dispersive shock waves in the regularized ILW and BO equations, as well as the ILW and BO systems.Since the results are broadly similar for ILW and BO models, only those corresponding to the BO-type models will be shown.
The formation of DSW's is associated to the propagation of a shock wave through a medium where dispersion effects dominate dissipation (see e.g.Hoefer and Ablowitz [7]).Under these conditions, the changes in the medium represented by the shock wave typically occur in the form of a structure involving two scales: a rapidly oscillatory wave scale and a slow modulational scale.Examples of DSWs appear in many disciplines such as Water Waves, Plasma Physics, Optics, etc. (see [22,23,8] and the references therein).
From the mathematical point of view, media featuring DSW's are often described by a system of conservation laws together with dispersive and dissipative terms.In the present study, consideration is given to the situation where dispersion dominates and so dissipation can be safely ignored.For completely integrable equations, the slow modulation of a rapidly oscillating wave can be described using Whitham's modulation theory, [22,23,24,8].In such problems, exact solutions for the corresponding Whitham modulation equations can be constructed in terms of the Riemann invariants of the conservation law.These in turn can be used to obtain, from suitable initial and boundary data, asymptotic representations of dispersive shock waves.This was applied, in [25] for the dispersive Riemann problem associated with the KdV equation; see also [23,26].Other important integrable systems have been similarly analysed; see for example [27] for the mKdV equation, [28] for the BO equation and [29] and [30] for the NLS equation.
In non-integrable systems, it is usually not possible to write the modulation equations in Riemann invariant form and several alternate strategies to study the formation of DSWs have emerged, e.g. the method of El [31] or the approximate method of Marchant and Smyth discussed in [32].
In many cases, DSW's are generated from discontinuities in the initial data.Our computational study will consider two initial-value problems for the rBO equation and the BO system with discontinuous initial data, namely the Riemann problem and the dam break problem.
Note first that if one neglects the dispersive terms in the BO system (1), the hyperbolic conservation laws emerge.This system can be written in a matrix-vector form, viz.
where v = (ζ, u) T and The distinct eigenvalues of A are and the associated eigenvectors, also known as Riemann invariants of the system (24), are They both satisfy the equations Thus, along the characteristic curves of the system, the quantities R ± are constant in time.

The Riemann problem
Consider first the BO system with Riemann-type initial data, which is to say the initial conditions are of the form with ζ ± , u ± satisfying the compatibility condition These initial data usually generate what is referred to as a simple DSW.
Taking the slightly easier case ζ + = u + = 0, it transpires that To simulate Riemann-type problems with the numerical method described in Section 2, the initial data must be slightly smoothed.Otherwise our Fourier-spectral method cannot handle the problem.Instead of (25), initial data of the form with is posited.Both ζ 0 and u 0 , while smooth, feature large gradients reminiscent of the pure Riemann data (25) (see Figure 15).Since ζ(−x, 0) = ζ(x, 0), x ∈ R and ζ(x, 0) decays to zero exponentially as |x| → ∞, the corresponding periodic IVP with initial conditions ( 28)-( 29) for the BO system and (28) for the rBO equation are integrated for x in a long spatial interval (−l, l).
For l = 1024 and step sizes h = 0.125 and ∆t = 0.001, the numerical simulation with these initial conditions is shown in Figures 15 and 16 respectively.In the case of the BO system, we observe the formation of a right-propagating simple DSW plus a dispersive rarefaction wave.Recall that in the absence of the dispersion, one expects a classical shock wave followed by a classical rarefaction wave.Since the compatibility condition (26) is not exact for the BO system, the initial conditions generate also leftward propagating dispersive tails.This does not seem to be the case of the rBO equation, as shown in Figure 16.Note that the formation of the simple DSW and the rarefaction wave does not look to involve additional structures.

The dam break problem
Next is introduced the so-called dam break problem.This is a special case of the Riemann problem in which u ± = 0 in (25).For the experiment below we will consider the periodic IVP for the BO system with initial conditions ζ(x, 0) on [−1024, 1024] given by (28) and u(x, 0) = 0.The resulting numerical simulation, see Figure 17, appears as a waveform consisting of two wave packets, each of which has a dispersive shock plus a rarefaction wave.The packets propagate symmetrically in opposite directions.

Concluding remarks
Studied here has been two asymptotic models for the propagation of internal waves along the interface of a two-layer fluid system, with the upper   layer bounded above by a rigid lid and the lower layer bounded below by a rigid, featureless horizontal bottom.The systems were derived in [2] in the Intermediate Long Wave and the Benjamin-Ono regimes, respectively.The Benjamin-Ono system is the limiting case of the Intermediate Long-Wave system when a very deep (theoretically infinite) lower layer is assumed.The corresponding one-dimensional pseudo-differential systems (1) have the same quadratic-type nonlinearities, while the linear parts for both of them are nonlocal.
In Part I of this study [1], the authors presented several numerical techniques to generate approximate solitary wave solutions of these systems whose existence, in the small-amplitude case, was proved in [5].A comparative study with solitary wave-solutions of related, unidirectional models of ILW and BO type was also developed.
In the present essay, we further study the solitary-wave solutions of (1), analyzing by computational means some aspects of their dynamics.To this end, the periodic IVP is solved numerically on a long enough spatial interval that it approximates well the IVP on R. The spatial approximation of our scheme is via a spectral Fourier-Galerkin method (implemented as pseudospectral; the resulting method is algebraically equivalent to a Fourier collocation scheme).This is coupled with an explicit, 4th-order Runge-Kutta time-stepping.The L 2 -convergence of the spectral semidiscretization was proved in [6] and the resulting fully discrete scheme was already used in [1] to check the accuracy of the computed solitary-wave profiles.These facts, together with the convergence studies rerported in [1], provide confidence in the accuracy of the simulations reported herein.
Specifically, our computational study is concerned with the following aspects: the dynamics of the solitary waves under small and large perturbations, the behaviour of overtaking and head-on collisions of solitary waves, and the resolution of initial data into solitary waves.We also compare computationally the bi-directional and uni-directional models corresponding to the same physical regime.The study is completed by analyzing numerically the formation of dispersive shock waves.Some of the conclusions of the study are the following: • Under small perturbations, the solitary waves seem to be nonlinearly stable, in the sense that an initially perturbed solitary-wave profile evolves into a modified solitary wave with slightly different amplitude and speed, along with a dispersive tail with two components, one trav-eling to the left and one to the right.The existence of these two dispersive groups was analysed via small-amplitude, plane-wave solutions of (1), linearized about the rest state.
• Testing the stability of solitary waves of the BO and ILW systems under large perturbations we found that they are still extremely stable.They evolve into a new solitary wave of an enhanced magnitude plus a dispersive tail of relatively small amplitude.
• The resolution property is illustrated with experiments involving the evolution of initial data of Gaussian type.The initial condition develops into a train of waves whose number depends on the energy of the initial Gaussian pulse (represented through its amplitude and wavelength parameters) and with a dispersive tail behind.
• The interactions are observed to be inelastic, something that is expected on account of the lack of Hamiltonian structures and the fact that there are only linear conserved quantities.Since the models are bidirectional, two types of collisions are considered.In the case of overtaking collisions, the behaviour of the waves after the interaction was compared to that of similar experiments for the uni-directional regularized BO equation.The main difference observed in the experiments seems to be the formation, in the case of the BO system, of N-shaped wavelets in the tail, traveling in the direction opposite to that of the emerging solitary waves.In the case of head-on collisions, the interactions develop dispersive tails propagating in both directions, as expected from the analysis of the small solutions of the linearized system, mentioned above.
• Comparisons between the unidirectional and bidirectional models when the initial data for the bidirectional model is well prepared for unidirectionality show remarkable similarities.This suggests that a theorem of the sort derived in [21] for surface waves likely holds for internal waves as well.
• The numerical experiments concerning the formation of DSWs were focused on the classical Riemann and dam break problems.In the case of the Riemann problem, suitable initial data generated a simple DSW along with a dispersive rarefaction wave, with an additional small dispersive structure traveling in the opposite direction.We checked that from the same initial condition, the BO and regularized BO equations developed similar structures but without the small additional dispersion.In the case of the dam break problem, the experiments suggest that from initial data with zero flow, the approximate solution of the BO system develops two DSWs plus rarefaction structures traveling in opposite directions.

Figure 1 :
Figure 1: Nonlinear stability of a BO-solitary wave with c s = 0.57 with amplitude approximately a = 3.9129, perturbed initially by about 10% in only the ζ-component of the system.

Figure 2 :
Figure 2: Nonlinear stability of a solitary wave with c s = 0.51 perturbed initially by about 50% in only the ζ-profile of the ILW-system.

Figure 3 :
Figure 3: Symmetric head-on collision of solitary-wave solutions of the BO system with c s = 3.2: (a)-(e): numerical approximation at times t = 100, 250, 313, 400, 1500; (f) is a magnification of the marked part of (e).

Figure 4 :
Figure4: Location of peak amplitudes of the numerical approximation during the symmetric head-on collision of the BO system, see Figure3.Notice the small, retarded phase shift after the interaction.

Figure 5 :
Figure 5: Overtaking collision of two solitary waves of the BO system (solid lines) and the rBO equation (dashed lines) with speeds c s = 3.5 and c s = 3.1.(a)-(f): numerical approximation at times t = 0, 4800, 5000, 5500, 6000, 8000.After the interaction, there is a small difference in the speeds of the larger emergent solitary waves.

Figure 6 :
Figure 6: The dispersive tails after the interaction of the solitary waves of the BO system appearing in Figure 5.(b) and (d) are magnifications of the marked part of (a) and (c), respectively.Notice the N-shaped, left-propagating wave.The next experiment provides an indication of the similarities between the BO-and the ILW-systems.An overtaking collision of two solitary waves

Figure 7 :
Figure7: The dispersive tails after the interaction of the solitary waves of the rBO equation from Figure5.Graphs (b) and (d) are magnifications of the marked part of (a) and (c), respectively.There is no N-shaped residual here.

Figure 8 :
Figure 8: The dispersive tail and the N-shaped wavelet produced after the overtaking collision of two solitary waves of the ILW system compared to the solution of the rILW equation with the analogous solitary waves, cf.Figure 5.
Figure 8: The dispersive tail and the N-shaped wavelet produced after the overtaking collision of two solitary waves of the ILW system compared to the solution of the rILW equation with the analogous solitary waves, cf.Figure 5.

Figure 12 :
Figure 12: Resolution of a Gaussian initial condition (21) with A = 1 and λ = 100 into solitary-wave solutions of the BO system.

Figure 13 :
Figure13: Evolution of a solitary wave of the rBO equation used as initial condition for ζ 0 with u 0 determined by(22) for the BO system.The maximum negative excursion of the solution as a function of time is also displayed.Observe that the negative excursion levels off as t grows, consistent with a traveling-wave structure.

Figure 14 :
Figure 14: Evolution of a solitary wave of the BO system used as an initial condition for the rBO equation.The maximum negative excursion of the solution as a function of time is also shown.