Hamiltonian approach to modelling interfacial internal waves over variable bottom

We study the effects of an uneven bottom on the internal wave propagation in the presence of stratification and underlying non-uniform currents. Thus, the presented models incorporate vorticity (wave-current interactions), geophysical effects (Coriolis force) and a variable bathymetry. An example of the physical situation described above is well illustrated by the equatorial internal waves in the presence of the Equatorial Undercurrent (EUC). We find that the interface (physically coinciding with the thermocline and the pycnocline) satisfies in the long wave approximation a KdV-mKdV type equation with variable coefficients. The soliton propagation over variable depth leads to effects such as soliton fission, which is analysed and studied numerically as well.


Introduction
It is well known that ocean wave dynamics displays very complex features, partly due to the effects of the bottom topography which quite often deviates from the convenient scenario of a flat bed. Waves over variable bottom are an active area of research, and various scales, geometries and approximations have been examined. For a survey of results we refer the reader, for example, to the review by Kirby [39], the book by Dingemans [27] and the references therein. The best known nonlinear wave models have been extended for fluids with uneven bottom as well: we recall here the KdV equation for long waves, for example, which has been generalised and studied thoroughly by Johnson [38], and the NLS equation for modulated waves, treated by Djordjević and Redekopp [28].
While most of the studies of internal waves involve irrotational flow, shear background currents have been included as well [13, 15, 17-19, 25, 36, 47]. However, the combined effects of a variable bottom topography, sheared currents and stratification on the arising internal waves are rather less investigated. We attempt to bring our contribution toward filling this gap by a derivation of a model equation (with variable coefficients) of KdV type which describes the interface in a flow with a variable depth and a flat surface in the presence of currents, density stratification and geophysical effects. A significant part in our endeavour is played by a variational approach based on the Zakharov's Hamiltonian formulation [52], the subsequent developments such as [1][2][3] as well as other irrotational scenarios, like surface and internal waves [24] or variable bottoms [16,23]. We advance here another feature of ocean dynamics: the presence of (non-uniform) underlying currents modeled by a specific choice of vorticity function.
The Hamiltonian formulation in our approach makes an extensive use of the Dirichlet-Neumann operators (DNO) [21,23,24] and provides a convenient setting for incorporating a series of features like interacting fluid layers, topography effects, stratification and underlying currents. In the presented study we illustrate the derivation and application of the Hamiltonian framework based on DNO for the case of uneven bottom. The illustrative example concerns internal wave propagation in the equatorial region in the presence of the Equatorial Undercurrent. The equatorial internal waves are special in a sense due to the effect of the Coriolis force. This effect keeps the waves propagating along the Equator like in a wave guide. Some further details could be found for example in [15,18]. We would like to note also the recent study by Guyenne [35] who proposed a numerical model for nonlinear surface waves in the presence of a vertically sheared current utilizing a Hamiltonian formulation combined with a series expansion of the DNO.
The rest of the paper is organised as follows. In Section 2 we formulate the model from the governing equations of the fluid mechanics. In Section 3 we describe the Hamiltonian representation of the evolution equations and then in Section 4 we present a long wave approximation which reduces to a KdV equation with variable coefficients which represents a generalisation of the flat bottom scenario. The effects on the solitons such as fission due to the variable depth are studied in Section 5. A special case with small or vanishing coefficient of the quadratic nonlinearity, taking into account the cubic nonlinearity, is analysed in Section 6. The analogues of the basic conserved quantities are given in Section 7. The Dirichlet-Neumann operators for the lower layer which is bounded by an uneven bottom is derived in the Appendix A. The numerical scheme for the finite-difference implementation of KdV-type equation with variable coefficients is outlined in Appendix B.

Equations of motion for an internal wave
We consider a two-dimensional water flow, moving under the influence of gravity, such that the x-axis is oriented along the horizontal direction and the y-axis is pointing vertically upwards. The time variable will be denoted with t. The flow is composed of two domains, Ω and Ω 1 , consisting of water with different constant densities and separated by a common interface, denoted y = η(x, t), which represents y x y = η(x, t) Figure 1. The system under study involving the two domains of different density (2.2), (2.1) and the four layers of different vorticity (2.5). The flat surface is y = h 1 , the elevation of the internal wave corresponds to the curve y = η(x, t) and the variable bottom corresponds to y = −h + β(x). The current profile (2.5) (for the special case of undisturbed interface η = 0) is given as well, the magnitude U is oriented along the horizontal axis, and the vertical dashed line corresponds to U = 0.
an interfacial internal wave. More specifically, we assume that, adjacent to the bed, the lower water domain is given as being situated below the near-surface region where β(x) is some given function, indicative of the unevenness of the bottom, while h and h 1 are positive constants, such that y = h 1 is the flat surface, and y = −h is the average depth of the seafloor, cf. Figure 1. Moreover, we assume that the internal wave has the properties the last one related to the fact that the average depth of the interface is at y = 0. 1 Taking into account the Equator's peculiar feature of behaving like a wave-guide we will look at two-dimensional inviscid and incompressible fluid motion confined near the Equator by the action of the Coriolis forces. We would like to point out that most of the physical variables that we employ here (like the density, the generalised velocity potentials or the components of the velocity field) display discontinuities across the interface y = η(x, t) that separates the two fluid regions. To make the reader observant of this aspect, we use the index 1 as a label for the upper layer. Whenever we refer to the overall physical variable without specification of the layer, we shall use bold face symbol.
Therefore, denoting with (u(x, y, t), v(x, y, t)) the velocity field, the equations of motion are Euler's equations where P = P (x, y, t) denotes the pressure, ω is the rotational speed of Earth, g is the gravitational acceleration and ρ denotes the density of the fluid which is assumed to be piece-wise constant being distributed as with the understanding that we look at a stable stratification, that is ρ > ρ 1 . Throughout the paper we will consider rotational water flows, and moreover, we presuppose that the vorticity is given as where U (y) represents the background shear current and has the following profile: the current is piece-wise linear with respect to y with the exception of the layer near the bottom where it decays to zero. More specifically, we introduce five layers as follows : Here γ and γ 1 represent constant vorticities in the corresponding sub-domains, κ is a constant component of the current. The choice of this current is made to model the Equatorial undercurrent which is formed by the winds blowing to the west [13,18,19] so that the current on the surface is negative. We note that the current is not continuous at y = η, and has a jump when η = 0 (continuous is only the normal component of the velocity field). Between the other layers U (y) is a continuous function such that U (−m) = 0, U (−l) = −lγ + κ. On Fig. 1 the current is illustratively shown for the situation of an undisturbed fluid when η ≡ 0. The layer between depths y = −m and y = −l is with prescribed current shape U (y) which could be chosen, for example, as necessary to match the field data. It could be shown that the energy per unit length of this layer is constant and thus the layer does not contribute to the equations on the interface (the pycnocline) [13]. However, the reconstruction of the velocity field in the body of the fluid would, of course, depend on U (y). The current gradually weakens towards the bottom layer, where it is zero, and there the motion is irrotational. With respect to the described stratification we have the following notations for the velocity field (2.6) u(x, y, t) := u(x, y, t), in Ω, In addition, we have the equation of mass conservation for incompressible fluid Complementing the equations of motion are the boundary conditions, of which the dynamic boundary condition (2.9) P = P atm on y = h 1 , (with P atm being the constant atmospheric pressure) decouples the motion of the water from that of the air. In addition, the impermeability of the surface boundary, of the interface and of the bottom, leads to the kinematic boundary conditions. On the surface and the interface respectively they are (2.10) v 1 = 0 on y = h 1 , while on the bottom we have the condition stating that the normal component of the velocity field vanishes.
To ease the notation, in the further developments of the paper, we introduce sub-indices b and s for the evaluations of physical quantities at the bottom y = B(x) and at the common interface y = η(x, t), respectively. In order to describe propagation of solitary waves we make the assumption that all considered functions η(x, t), ϕ(x, y, t), ϕ 1 (x, y, t), β(x) are in the Schwartz class with respect to the x variable, that is declining fast enough when x → ±∞ (for all values of the other variables), which we denote by η(x, ·) ∈ S(R), etc.
The equation of mass conservation (2.8) ensures the existence of a stream function determined up to an additive term that depends only on time, by in Ω, u 1 = ψ 1,y , v 1 = −ψ 1,x , in Ω 1 We will impose the condition that the stream function is continuous on the interface y = η(x, t), condition which can be expressed as (2.14) ψ(t, x, η(x, t)) = ψ 1 (t, x, η(x, t)).
The latter condition implies that the normal velocity components equal across the interface y = η(x, t).
We introduce now the (generalized) velocity potential by means of The kinematic boundary conditions (2.10) and (2.11) can now be written as and respectively, as The following notation will be used later on in the paper. Namely, we set Euler's equations can be expressed by means of the stream function and of the generalized velocity potential as (2.19) ∇ ϕ t + 1 2 |∇ψ| 2 + P ρ − (γ + 2ω)ψ + gy = 0.
We have therefore From the condition P = P atm on the top surface and availing also of the Schwartz property of the stream function ψ 1 and of the velocity potential ϕ 1 we infer that P atm /ρ 1 + gh 1 = f 1 (t) for all t. Moreover, utilizing the continuity of the pressure at the interface y = η(x, t) and making the choice f (t) := ρ 1 f 1 (t)/ρ we obtain the Bernoulli type equation (2.20) ρ (ϕ t ) s + |∇ψ| 2 where χ := ψ(t, x, η(x)) = ψ 1 (t, x, η(x)) is the stream function evaluated at the interface.

The Hamiltonian functional and the Hamiltonian formulation
We start this section by indicating the most physical choice for the Hamiltonian functional which is the total energy of the flow. This functional is then written in terms of the canonical variables: the interface defining function η and a suitable combination of the velocity potentials ϕ and ϕ 1 evaluated on the interface. We then introduce an almost-Hamiltonian formulation of the considered wave-current system that becomes Hamiltonian in the absence of sheared underlying currents (that is, for γ = γ 1 = 0). Lastly, it will be proven that an appropriate change of variables renders the almost-Hamiltonian formulation into a bona fide Hamiltonian one.
3.1. The Hamiltonian. The outset of this section's endeavour is the evaluation of the kinetic and of the potential energy for each domain. For the lower domain the potential energy is Observing that (3.1) u 2 + v 2 = (ϕ x + U (y)) 2 + ϕ 2 y = div(ϕ∇ϕ) + 2ϕ x U (y) + (U (y)) 2 we use Green's Theorem (Divergence Theorem) and find that the kinetic energy of the lower layer can be calculated as where n s = (−η x , 1)/ 1 + η 2 x is the outward-pointing unit normal vector (with respect to Ω) to the wave surface, while n b = (B x , −1)/ 1 + B 2 x is the outwardpointing unit normal vector on the bottom. Since Let us introduce the Dirichlet-Neumann operators G ij (β, η) given by This operator depends on the bottom variations through β(x). The details of the computation of these Dirichlet-Neumann operators is given in the Appendix.
For the upper layer, similarly, The minus sign is because the outward normal for the domain Ω 1 is −n s . Recall that Some of the integrals above are not convergent due to the constant densities at infinity. However, the Hamiltonian is the energy difference from the unperturbed state, whose energy is the energy of the current (which is infinite due to the infinite domain): therefore the Hamiltonian will be evaluated from (3.11) The function above is not yet the Hamiltonian. As the Bernoulli equation (2.20) suggests, the momentum-type variable, conjugate to the coordinate-type variable η is [1-3] Using (3.10) and (3.9) as well as (2.17) we have From (3.12) and (3.14) we express Φ and Φ 1 in terms of ξ and η : Utilizing (3.14), (3.17) and (3.11) we express the Hamiltonian in the form This expression formally coincides with the expression for the flat bottom [14], apart from the fact that G(η, β) depends now on the bottom topography. Considering the case where γ 1 = γ and µ = 0 and the Hamiltonian reduces further to 3.2. Hamiltonian structure. The Hamiltonian structure for two-layer domains allowing for currents and an internal wave is derived in [11,12] for the case of a flat bottom. Moreover, the Hamiltonian formulation for the irrotational scenario for surface waves over a rough bottom was derived in [23], developing the perturbative technique for the Dirichlet-Neuman operators for non-even bottom. The currents in the layer −m ≤ y ≤ −l do not contribute to the Hamiltonian, as it could be seen from (3.7) and (3.18). Indeed, as noted in [13,36] the motion of the interface is affected only by the motion in the layers adjacent to the interface, so that the equations of motion of the internal wave (2.17) and (2.20) can be represented in the (non-canonical) Hamiltonian form is a constant and is the stream function, evaluated at y = η(x, t), (see [11,20] for details).
Introducing the variable u = ξ x one can write down (3.20) in the equivalent form Remark 3.1. There is a formal transformation of the equations (3.20) into a canonical form by the following change of one of the variables (cf. [11,12,46]) Therefore, the system (3.20), which could be also represented as is Hamiltonian too [11,12].
and hence ζ(x, t) ∈ S(R). We note that In what follows we will employ these equations with an approximation for the Hamiltonian functional H.

3.3.
Series expansion of the Dirichlet-Neumann operator. The Hamiltonian (3.11) depends on the Dirichlet-Neumann operator This is a self-conjugate operator. Some details about G ij at different orders can be found in Appendix A.1. The operator can be expanded over the powers of η and β. Let us now introduce appropriate scales. The interfacial waves are assumed of small amplitude, relative to h, i.e. |η max |/h = ε 1. The bottom variations are also considered small, but |β max |/h of order˜ ≤ ε 1/3 . The rationale of this choice will become evident later (see the expansion (3.28)). We point out that the magnitude of the bottom variations is not fixed by ε, varying within the described limits, with all derivations being valid for the flat bottom as well. In order to keep terms up to the order of ε we keep in the expansion η 0 , η 1 and β 0 , . . . , β 3 . Therefore, in all expansions we are keeping the contributions from the following entries: ij , whose orders are explicitly is of order˜ ε and is neglected. The next assumption is the assumption of the slow variations of the bottom profile. Mathematically, we assume that β = β(εx). Then the commutator of β and the differentiation operator D := −i∂ x is proportional to εβ (εx) which itself is of order ε. Thus, if we keep only terms of order ε, we can write ε a Dβ ≈ ε a βD (where 0 < a ≤ 1) since the difference ε a Dβ − ε a βD ∼ ε a+1 ε and could be neglected. In other words, with the exception of the leading order term, we can interchange βD and Dβ. The expansion involves also the long-wave parameter δ = h/λ 1 where λ is the typical wavelength. Since k = 2π/λ is the wave number, sometimes we write also symbolically δhk instead of hk to remember the fact that the quantity is of order δ. Moreover, we write δhD, instead of hD since k is the eigenvalue of D = −i∂ x when acting on functions representing plane waves exp(ikx). With these assumptions the truncated expansion is is the local depth and X = εx indicates that the bottom depth varies slowly with x. This of course coincides with the result from [16] which was obtained following a slightly different approach based on the framework from [23]. Assuming that O(h) = O(h 1 ), for the operator G 1 we have as usual [24] (3.29) or with the hyperbolic tangent functions expanded,

The long wave approximation
We are concerned in this section with the KdV-like long-wave regime which arises when the relation between the scales is ε ∼ δ 2 , ξ = O(δ), and then both u and η are of order δ 2 . Then for the operator B we have (4.1) which entails that the approximate Hamiltonian whose expansion, including terms of order δ 6 , is: where, using the notation X = δ 2 x, The Hamiltonian equations (3.23) for the Hamiltonian (4.2) in terms of η and u = ξ x are (note that the equations written with scales will bring a factor of δ 2 for each Hamiltonian variable, which will compensate the overall factor δ 4 of the Hamiltonian) The x-derivatives of α k (X) produce quantities of smaller order. Substituting the explicit form of α 5 we have (4.5) Since ω = 7.3 × 10 −5 rad/s, κ ∼ 1 m/s, then g 2ωκ and the 2ωκ term will be neglected.
In the leading order η t + κη x + (α 1 u) x = 0, The monochromatic solutions for η and u can be obtained in the form where c(X) is the wave speed, which depends on the "slowly varying" variable X. From (4.6) -(4.7) it is straightforward to obtain the following quadratic equation for the wave speed c: The solutions are and c − ≈ −3.39 m/s (left running waves). It is evident also that c(X) is of the same order as κ. Another observation is that the presence of vorticity does not change considerably the wave speed, which for the irrotational case (γ = 0, γ 1 = 0) for example, is c ± ≈ 1 ± 4.37 m/s. In contrast, the wave speed of the surface waves (whose effect is neglected here) depends significantly on the vorticity near the surface [17].
As in the previous studies, following [16,37,38], in addition to the variable X, we introduce the characteristic variable in the form where R(X) is a function such that R (X) = 1/c(X). The (x, t) coordinate partial derivatives change according to (4.11) The equations then can be transformed from (x, t) variables to the slow variables (θ, X). This way, of course, two sets of equations arise (for the left and for the right running waves). The equations written in terms of the new variables are From the second equation and (4.9) Therefore, in the leading order we have Next we substitute (4.14) in (4.12) and then we substitute (4.15) in the terms of order δ 2 to obtain the single equation for η which is a KdV-type equation [41] with variable coefficients that depend on functions, slowly varying with X. From (4.8) one can establish a connection between c X and α 1,X : One possible limit of equation (4.16) is the irrotational case where ω = γ = γ 1 = 0. Hence, it follows that Γ = α 4 = α 6 = 0, c 2 = g(ρ − ρ 1 )α 1 and thus, equation (4.16) acquires the form which showed up in non-dimensional form in [29] lacking, however, details of its derivation.
Then the equation (4.16) becomes This situation corresponds to surface waves over one layer of fluid and has been analysed in [17]. Further reduction of (4.20) could be obtained for the irrotational case with κ = 0, b = c 2 /g. The equation acquires the form (4.21) (2cη X + c X η) + c 2 3g 2 η θθθ + 3g c 2 ηη θ = 0, which is the equation derived by Johnson [37,38] for one layer of irrotational fluid over variable bottom; we refer the reader to [16] for its derivation by means of a Hamiltonian approach. Therefore equations (4.16), (4.20) provide generalisations of Johnson's equation. These equations resemble the KdV equation, the important distinction being that they exhibit variable coefficients. While the integrability of KdV type systems is well established in numerous investigations for a long time [49], it seems that these particular models are, for a general choice of b(X), not integrable [50] .

Fission of solitons moving over a step-like bottom
We study an example where the step-like bottom is modeled by The constants are taken as the actual values in the SI system of units as follows: h = 2000, h 1 = 200, g = 9.81 ρ = 1037, ρ 1 = 1026. The constants α,β are given in each case. First, we study the irrotational case (4.18). The bottom threshold is located at X = 0, the initial condition is an exact KdV soliton coming from When α > 0 the soliton moves from deep to shallow regions, when α < 0 -from shallow to deep. The depth profile is given on Fig. 2a.
When κ = 0, considering right-moving waves, we have The initial condition is the soliton that satisfies the unperturbed KdV equation with b = b 0 = h(1 + α) : where all quantities with sub-index 0 are evaluated with b = b 0 = h(1+α), including c 0 = [c(X)] b=b 0 : Writing the previous KdV equation (5.3) in the form [c(X)] 2 [2c(X) + α 1 Γ] are constants, we note that it is well known that it has an one-soliton solution where K, θ 0 are arbitrary constants. Therefore, we can take an initial condition at X = X 0 < 0 (ideally modelling a soliton coming from X → −∞ ) which is an exact solution of the equation for X → −∞: Actually for the choice of X 0 it is sufficient that tanh(βX 0 /h) ≈ −1, e.g. |βX 0 /h| > 2.5 ) noting that X plays the role of a "time" variable. The initial soliton profile is given on Fig. 2b. The explanation of the soliton fission when the initial soliton (5.7) reaches the threshold at X = 0 follows from the Quantum Mechanical theory of the Pöschl-Teller potential sech 2 (Kθ), cf. [31]. The Lax operator for the KdV-equation has the form of a Schrödinger equation type spectral problem for the eigenfunction ψ(θ) with potential η(θ, 0) = EK 2 sech 2 (Kθ), see the details for example in [49]. Then, since the spectral problem is iso-spectral, that is X-independent, we note the following. The eigenfunction ψ corresponds to a potential (KdV-solution) with N discrete eigenvalues if E = −N (N + 1)/2. This corresponds also to an N -soliton KdV solution. Let us now take an initial condition at X 0 such that η(θ, X 0 ) is the one-soliton (N = 1) solution for X 0 → −∞, and N > 1 soliton solution for X → ∞. Then taking into account the constant coefficients of the corresponding KdV equations in the two situations, (when X → ±∞) we have an equality of the coefficient in front of the sech 2 -potential, which is X-independent, Some more details are available for example in [16]. Introducing the notation b * = h(1 − α) in the irrotational case (4.18) we have This formula gives N = 2.07 while in reality we observe at least 3 solitons. The numerical solution is presented at Fig. 3. Actually, formula (5.10) could be improved by noticing that with an integrating factor, introducing F (θ, X) = c(X)η(θ, X) the equation      Fig. 2a and 2b, −1000 < X < 6000, ω = γ = γ 1 = 0. The horizontal axis is for the variable θ. The waterfall plot corresponds to increasing values of the variable X.
which also appears in [29] where it is obtained from the arguments of Johnson [38], see also [16]. In particular, for our data it gives N = 2.14. The discrepancy could be due to the fact that the derivation of (5.8) is based on the equality of the KdVamplitudes at the moment of the hitting of a step-like threshold. The threshold in our case however is not sharp, it is modelled by the smooth tanh function. In addition, over the region of the obstacle the equations are not exactly KdV equations, since their coefficients depend in general on the bottom variations in the region of the smooth threshold. Moreover, at the obstacle there is always a reflected wave, which is not taken into account. So the formula for N should be considered only as an estimate.
Next, we study numerically the situation with nonzero vorticities. We take γ = 0.1, γ 1 = −0.1 and ω = 0, being very small in comparison to the other vorticities. The depth b(X)-profile is the same as on Fig. 2a however when K = 0.002, the amplitude of depression of the initial condition rises to 100 meters (which is   not unusual for internal waves) as shown on Fig. 5a. The soliton fission is less pronounced giving two solitons, as it could be seen from Fig. 4 and Fig. 5b.
The approximate ratio (5.9) gives N ≈ 1.98, which is in an agreement with the results. One could hope for an improved formula like in the irrotational case, however the integration factor is apparently not in a simple form and we are going to limit ourselves with the approximation (5.9). We provide the dependence of the ratio (5.9) on the step magnitude α on Fig. 6a and 6b. It is evident that in both cases the maximum of the ratio is reached for α ≈ 0, 7 which is already quite an extreme value. In the rotational case the ratio N (N + 2)/2 barely reaches the value of 3 which corresponds to just two solitons, which are actually observed in the numerical experiments. The findings show that the solitons of the internal waves are quite robust (in comparison to those of the surface waves, [16,17]) and remain stable for relatively mild bottom variations.
In the case α < 0, as predicted from the theory, no soliton fission is observed, the initial soliton very slowly decays loosing its energy through waves of radiation in the region X > 0.

Special case when the coefficient of the nonlinear term is of smaller order or vanishing
There is a particular case when the coefficient of the ηη θ term in (4.16), is of order of ε or smaller. In fact, there are situations with values of X and the parameters such that this coefficient could be zero or very closed to zero, so that the main nonlinearity term is η 2 η θ . This could be seen immediately in the irrotational case (4.18) when ρ 1 b 2 = ρh 2 1 . In this section we explore the situation with small or vanishing coefficient (6.1). Since the dispersive term η θθθ matches the order of the nonlinear term η 2 η θ , this indicates that the scaling is η and u of order δ, that is ε and δ of the same order. Since we have taken X = δ 2 x, where the scale of the wave variations is over δx, we have now X = ε 2 x. The variable changes to the characteristic variables like in (4.10), (4.11) become The extended Hamiltonian with terms of leading order ε 2 up to order ε 4 could be obtained from (3.18) and the expansions (3.28) and (3.30) with δ = ε where the new terms of order ε 4 are with coefficients 3) The two Hamiltonian equations lead to the analogues of (4.12) and (4.14), however with extra terms, arising from the contributions of the terms with β 1 , β 2 and β 3 in the Hamiltonian: From (6.5) we have the following relation in the leading order and from (6.5) and (6.6) we obtain We substitute (6.5) in (6.4) and then in the so obtained equation we eliminate u with the help of (6.7), thus obtaining the following equation for η : Finally, if the order of the coefficient (6.1) itself is of order ε or smaller, we observe that all terms are of the same order, hence giving a mKdV-type equation The equation (6.9) generalises (4.16). Indeed, it works in the scaling used in the derivation of (4.16) as well, because in this scaling the term with η 3 could be neglected. Other authors also suggest the mKdV-type equation as a suitable generalisation avoiding the problem with the vanishing term in front of the η 2coefficient, [33,42]. Finally we point out that the mKdV equation with constant coefficients is integrable, [49] so that, one can embark on developing soliton perturbation theory for (6.9).

Conserved quantities
With an integrating factor (7.1) the equation (6.9) acquires the following form for the quantity We note that for Γ = κ = 0, the above relationship is just E = c(X)η. The conserved quantity (mass conservation) from (7.2) is The interpretation of this type of results is not straightforward -see the comments in the variable bottom section of the Johnson's book [37]. The mass conservation makes sense for the full system and not just for the solution describing the pycnocline. Multiplying (7.2) by E leads further to and therefore we obtain an analogue of the "energy" conservation, which reads as (7.7) 1 2 R E 2 (X, θ)dθ = const.

Conclusions
The main achievement of this paper is the derivation of a KdV type equation (4.16) which describes the propagation of interfacial internal waves in two-layer domains bounded below by a variable bottom and above by a flat surface. Our analysis includes the shear currents in the two domains and hinges on a consistent derivation of the Dirichlet-Neumann (DN) operators in the Boussinesq approximation. While the final result for the bottom-dependent DN operator (3.28) recovers the one from Craig et al. [23], the setup of its derivation opens up new possibilities towards an application to a multi-layer system of fluids, which will be explored in forthcoming publications. The bottom-dependent DN operator allows the application of the "nearly"-Hamiltonian formulation, developed for the configuration of two fluid layers in [12,14] following the DN approach of Craig et al. [24].
An example for a possible realistic situation are the equatorial waves and currents in the equatorial Pacific Ocean, where the so-called Equatorial Undercurrent resides and where the abyssal hills are the most abundant seabed structures near the equator. These Pacific Ocean hills are typically 50-300 m in height, with a width of 2-5 km and a length of 10-20 km [26]. Other seabed structures are the seamounts which are higher, but with horizontal dimensions of the same order, that is, bottom structures with horizontal diameters of 2-20 km are typical. Since the bottom length scale over wavelength ratio is of order ∼ 1/δ, which could be a factor of 2-10, the modelling setup is a realistic scenario for waves of 0.5-5 km wavelength.
The general equation (4.16) in various limits leads to several known simplified cases with variable bottom, like the irrotational case (4.18), the single layer case with background current (4.20) and without current (4.21), which goes back to the well-known work of Johnson [38]. In addition, it has been noted that for some values of the parameters the coefficient of the nonlinear term of the model equation might be close to zero, rendering the next order term η 2 η θ of increased significance. The scaling for this special case is identified and a model equation of mKdV type is derived (6.9), utilising the Hamiltonian method.
Wave-breaking of solitary waves is another very significant and interesting topic [51]. However, its analytical studies will require modelling beyond the KdV and Boussinesq-type models.
where n s = (−η x , 1)/ 1 + η 2 x and n b = (β x , −1)/ 1 + β 2 x are the outward unit normal vectors corresponding to the interface and to the bottom, respectively. Therefore, the last equality is written as In what follows we work with a fixed wave number k and a special associated harmonic function ϕ k (x, y) = (a(k)e ky + b(k)e −ky )e ikx . The corresponding values on the surface and at the bottom are 2 Φ k (x) = (ae kη(x) + be −kη(x) )e ikx , Φ b,k (x) = (ae −k(h−β(x)) + be k(h−β(x)) )e ikx Expanding we have (summation over j ≥ 0) Since ϕ k (x, y) = (ae ky + be −ky )e ikx it follows that We assume that each entry G mn could be written as an infinite series where G (p,q) mn (η, β) is a homogeneous expression of its arguments in a sense that for any two constants C 1 , C 2 The coefficients a(k) and b(k) in the expressions above are arbitrary, therefore we can equate the corresponding coefficients in (A.7) and we have We proceed to compute now the operator G (1,0) , that is the one which exhibits powers of type η 1 and β 0 . First, we have From the above system we derive We would like to note now that the operator corresponding to iββ Hence, we have (A.32) Availing of the identity we obtain that Although it is not obvious, the operator 1 2 β 2 D 2 β − 1 6 β 3 D 2 is self-conjugate, see the identity (A.46), thus G = ikβ x (βk) 2 ae −kh + be kh 2 − k(βk) 3 ae −kh + be kh 6 . We assume an initial condition of the form of one-soliton solution (5.7). We consider an uniform mesh in the interval [−L 1 , L 2 ] θ i = (i − 1)∆θ, spatial step ∆θ = (L 1 +L 2 )/(N −1), and X n = n∆X, where N is the total number of grid points in the interval and ∆X is the time increment. Respectively, η n i and η n+1 i denote the value of η at the i -th spatial point and "time" stages X n and X n+1 correspondingly. We construct the following nonlinear difference scheme which is convergent: Such a scheme is stable when time-stepping with respect to the physical time, provided the iterative procedure to resolve the nonlinear terms is convergent. Since the scheme given by (B.3) cannot be implemented directly, because it is nonlinear we follow the idea of [9] to introduce internal iterations, namely This way, for the current iteration of the unknown function (superscript n+1; k+1) we have an implicit system with five-diagonal band matrix. We begin from an initial condition η n+1,0 = η n and conduct the internal iterations (repeating the calculations for the same time step (n + 1) with increasing value of the superscript k) until convergence.
Note that the initial condition is a very good guess which is within O(∆X) of the sought solution for η. This makes the convergence of the internal iterations very fast. We have performed the numerical experiments to verify this fact. Even for very large values of the time increment ∆X we have not encountered any instability of the internal iterations. In each case under consideration we selected ∆X such that no more than six internal iterations were required to reach the precision of 10 −12 . After the internal iterations converge, one gets the solution of the nonlinear scheme by setting η n+1 ≡ η n+1,k+1 . In such a way we have fully implicit, nonlinear and conservative scheme. For the inversion of the five-diagonal N × N matrix we use a generalized algorithm based on Gaussian elimination with pivoting (for details see [8,10,48]).
The scheme was thoroughly validated through the standard numerical tests involving halving the spacing and time increment. The global truncation error of time approximation was verified as by Runge principle and confirmed the second order of accuracy in time. In a similar fashion we found that the global spatial truncation error is also second-order, O((∆θ) 2 ).