# On a Coupled System of Shallow Water Equations Admitting Travelling Wave Solutions

Authors: D. Burini; S. De Lillo; D. Skouteris
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

We consider three inviscid, incompressible, irrotational fluids that are contained between the rigid wallsy=−h1 and y=h+H and that are separated by two free interfaces η1 and η2. A generalized nonlocal spectral (NSP) formulation is developed, from which asymptotic reductions of stratified fluids are obtained, including coupled nonlinear generalized Boussinesq equations and (1+1)-dimensional shallow water equations. A numerical investigation of the (1+1)-dimensional case shows the existence of solitary wave solutions which have been investigated for different values of the characteristic parameters.


---

## Body

## 1. Introduction

Since some early studies started in the 1950s (see [1, 2] and refs. therein), much interest has been devoted over the years to the flow of stratified fluids, both from a theoretical and from applied point of view. In particular, multifluid flows characterized by gravitationally stable density interfaces are a useful mathematical model in geophysics and engineering [3–5]. On the other hand, many analytical studies have been dedicated to the equations that describe the system of two ideal fluids, separated by a free interface [6]. Indeed, a great effort has been devoted to modeling the evolution of the internal waves in a two-fluid system and to deriving asymptotic reductions (see [7] for a recent review). Among them, of particular relevance are the Benjamin-Ono (BO) equation [8, 9] and the intermediate long wave (ILW) equation [10, 11]. Moreover, in [12] model equations governing the evolution of fully nonlinear long waves at the interface of two immiscible ideal fluids were derived, while the weakly nonlinear limit was treated in [13]; finally in [14] a Hamiltonian formulation of the two-fluid system was derived and through a perturbative theory a systematic analysis of the long waves scaling regimes was carried out. On the other hand, in recent years, a nonlocal formulation of water waves for both (1+1) and (2+1) dimensions was presented in [15], where the original equations with unknown boundary conditions are replaced by an integrodifferential equation and a nonlinear partial differential equation, both of which are formulated in a known domain. The nonlocal formulation obtained in [15] is derived from the general approach to studying boundary value problems for linear and nonlinear PDEs introduced in [16]. A crucial role in such approach is played by a nonlocal equation, called the global relation [16]. A generalization of the results obtained in [15] was presented in [17], where a nonlocal formulation was derived, governing two ideal fluids separated by a free interface and bounded above either by a rigid lid or by a free surface [17]. Due to the dependence on a free spectral parameter, the corresponding equations are usually called the nonlocal spectral (NSP) equations of the two-fluid system. The NSP equations were particularly useful for deriving asymptotic approximations; we wish to point out an asymptotically (2+1)-dimensional generalization of the intermediate long wave (ILW) equation reported in [17] which includes the KP equation and the Benjamin-Ono equation as limiting cases. Numerical investigations indicated the existence of lump type solutions, with a speed versus amplitude relationship shown to be linear in the shallow, intermediate, and deep water regime.However, to the best of our knowledge, the phenomenological models for more than two fluids mentioned at the beginning of this section have not been paralleled by any analytical study. This prompted us to develop a generalization of the NSP formulation to the case of three ideal fluids, separated by two free interfaces and limited above by a rigid lid. Namely, we consider three inviscid, incompressible, irrotational fluids that are confined between the rigid lidsy=-h1 and y=h+H and are separated by two free interfaces η1(x,t) and η2(x,t). We derive in the following the NSP equations governing the evolution of the three-fluid system. Specifically, define the functions q(x,t)=φ(x,η1(x,t),t), Q(x,t)=Φ(x,η1(x,t),t), P(x,t)=Φ(x,h+η2(x,t),t) and Q~(x,t)=Φ~(x,h+η2(x,t),t), where φ(x,y), Φ(x,y), and Φ~(x,y) are, respectively, the velocity potentials in the lower, the intermediate, and upper layer (see Figure 1).Figure 1
Three fluids with two free internal interfaces.The NSP formulation is given by the following equations.First layer (bottom)(1)∫R2‍eikxcosh⁡kh1+η1η1tdx=i∫R2‍eikxsinh⁡kh1+η1kk·∇qdx,intermediate layer(2)∫R2‍eikxsinh⁡⁡kη2η2tdx-∫R2‍eikxsinh⁡kη1-hη1tdx=-i∫R2‍eikxcosh⁡⁡kη1-hkk·∇Qdx+i∫R2‍eikxcosh⁡⁡kη2kk·∇Pdx,(3)∫R2‍eikxsinh⁡⁡kη2+hη2tdx-∫R2‍eikxsinh⁡kη1η1tdx=-i∫R2‍eikxcosh⁡⁡kη1kk·∇Qdx+i∫R2‍eikxcosh⁡(k(η2+h))k(k·∇)Pdx,third layer (top)(4)∫R2‍eikxcosh⁡⁡kη2-Hη2tdx=i∫R2‍eikxsinh⁡kη2-Hk(k·∇)Q~dx,Bernoulli’s equations(5)ρ1qt+12∇q2+gη1-η1t+∇q·∇η122(1+∇η12)-ρ2Qt+12∇Q2+gη1-η1t+∇Q·∇η122(1+∇η12)=σ1∇·∇η11+∇η12,ρ2Pt+12∇P2+gη2-η2t+∇P·∇η222(1+∇η22)-ρ3Q~t+12∇Q~2+gη2-η2t+∇Q~·∇η222(1+∇η22)=σ2∇·∇η21+∇η22.In (1)–(5), x=(x1,x2), k=(k1,k2), kx≡k·x, where k is a free spectral parameter and the constants g, σ1, and σ2 denote gravity, surface tension of η1(x,t), and surface tension of η2(x,t), respectively.In the next section, we derive the above equations, starting from the classic equations governing three ideal fluids separated by two free interfaces bounded above and below by rigid lids. Moreover, we derive conservation laws and integral identities for the three-fluid system for the NSP formulation.Section3 is devoted to the derivation of the weakly nonlinear equations: after a suitable nondimensionalization of the variables, we obtain the reduction to a system of shallow water equations in the weakly nonlinear limit.In Section4, under the assumption of maximal balance, we introduce travelling wave variables moving only to the right and we study only the (1+1)-dimensional case. In terms of the new variables, we finally obtain a system of coupled nonlinear shallow water equations which we study numerically in terms of the parameters entering the theory and that we show to admit solitary wave solutions.


## 2. A Nonlocal Spectral Formulation of the Classic Three-Fluid Equations

### 2.1. A Weak Formulation of the Classic Three-Fluid Equations

We recall the classic equations governing three ideal fluids separated by two free interfacesη1 and η2 and bounded above and below by rigid lids. It is assumed that the lower fluid is of density ρ1, the intermediate fluid is of density ρ2, and the upper fluid is of density ρ3 with ρ1>ρ2>ρ3. The equations are given in terms of the interface variables and the velocity potentials φ, Φ, and Φ~ associated with the lower, intermediate, and upper fluid domains, respectively, (6)Δφ=0,for-h1<y<η1,(7)φyy=-h1=0,(8)η1t=(φy-∇xφ·∇xη1)y=η1,(9)ΔΦ=0,forη1<y<h+η2,(10)η1t=(Φy-∇xΦ·∇xη1)y=η1,(11)η2t=(Φy-∇xΦ·∇xη2)y=h+η2,(12)ΔΦ~=0,forh+η2<y<H+h,(13)Φ~yy=H+h=0,(14)η2t=(Φ~y-∇xΦ~·∇xη2)y=h+η2,(15)ρ1φt+12∇φ2+gη1-ρ2Φt+12∇Φ2+gη1=σ1∇·∇η11+∇η12,on⁡y=η1,(16)ρ2Φt+12∇Φ2+gη2-ρ3Φ~t+12∇Φ~2+gη2=σ2∇·∇η21+∇η22,on⁡y=h+η2.We also require that |∇φ|→0, |∇Φ|→0 and |∇Φ~|→0 as |x|→∞. In the previous equations, the constants g, σ1, and σ2 denote gravity, the surface tension associated with the free interface η1, and the surface tension associated with the free surface η2. Equations (6), (9), and (12) express the fact that fluids are divergence and curl-free. Equations (15)-(16) express that the jump in pressure across an interface is balanced by surface tension.We now obtain a weak formulation of (6)–(16), expressed in terms of η1, η2, q, Q, P, and Q~. Let us define the domains (17)D(η1)≡{(x,y)∈R2∣x∈R2,-h1<y<η1(x,t)},Dη1,η2≡(x,y)∈R2∣x∈R2,η1(x,t)<y<h+η2(x,t)(x,y)∈R2∣x∈R2,,Dη2≡(x,y)∈R2∣x∈R2,h+η2(x,t)<y<h+H(x,y)∈R2∣x∈R2.For any bounded, harmonic functions ψ,Ψ,Ψ~, defined in D(η1), D(η1,η2), and D(η2), respectively, and satisfying (18)Δψ=0,for(x,y)∈D(η1),ψyy=-h1=0,ΔΨ~=0,for(x,y)∈D(η2),Ψ~yy=H+h=0,∇ψ,∇Ψ,∇Ψ~⟶0,asx⟶∞,the following identities hold at any given time t: (19)∫R2‍ψy=η1η1tdx-∫R2‍q(ψy-∇xψ·∇xη1)y=η1dx=0,(20)∫R2‍Ψy=h+η2η2tdx-∫R2‍Ψy=η1η1tdx+∫R2‍QΨy-∇xΨ·∇xη1y=η1dx-∫R2‍P(Ψy-∇xΨ·∇xη2)y=h+η2dx=0,(21)∫R2‍Ψ~y=h+η2η2tdx-∫R2‍Q~Ψ~y-∇xΨ~·∇xη2y=h+η2dx=0.Equations (19) and (21), together with the Bernoulli equations (15) and (16), constitute a nonlocal system of equations describing the three-fluid system.


### 2.2. Derivation of the Nonlocal Spectral Equations

We derive from (19) the NSP equation (1) given in the Introduction. To motivate the derivation, suppose that ψ(x,y) satisfies(22)Δψ(x,y)=0,ψy(x,-h1)=0.Define ψ^(k,y) by (23)ψ^(k,y)=12π2∫R2‍e-ikxψ(x,y)dx.Then, conditions (22) turn into the following: (24)ψ^yy(k,y)=k2ψ^(k,y),ψ^y(k,-y)=0.The solution to the above ODE is given by (25)ψ^(k,y)=ξ^kcosh⁡ky+h1,reverting back to physical space, (26)ψ(x,y)=∫R2‍eikxξ^kcosh⁡ky+h1dk.The previous equation shows that any harmonic function ψ whose y-derivative at y=-h1 vanishes can formally be written as a sum of functions eikxcosh⁡(|k|(y+h1)). Therefore, by linearity it suffices to require that (19) holds for the parametrized family of functions (27)ψk(x,y)=eikxcosh⁡ky+h1.Substituting ψk(x,y) into (19), we obtain (28)∫R2‍eikxcosh⁡kη1+h1η1tdx=∫R2‍qeikx-icosh⁡kη1+h1(k·∇)η1ksinh⁡⁡kη1+h1-icosh⁡kη1+h1(k·∇)η1eikx-icosh⁡kη1+h1(k·∇)η1ksinh⁡⁡kη1+h1dx.We simplify (28) by noting that(29)keikxsinh⁡kη1+h1-ieikxcosh⁡⁡kη1+h1k·∇η1=-i∇·kksinh⁡kη1+h1eikx.Using this identity in (28) and integrating by parts, we finally get the nonlocal spectral equation (1) as follows: (30)∫R2‍eikxcosh⁡⁡kη1+h1η1tdx=i∫R2‍eikxsinh⁡kη1+h1k(k·∇)qdx.Along the same lines used for the first layer, let us derive (2)-(3) from (20).We define the basic functions(31)Ψk(x,y)=eikxsinh⁡ky-h,and putting Ψk into (20) gives(32)∫R2‍eikxsinh⁡⁡kη2η2tdx-∫R2‍eikxsinh⁡kη1-hη1tdx=-∫R2‍eikxQkcosh⁡⁡kη1-h-isinh⁡kη1-hk·∇η1dx+∫R2‍eikxPkcosh⁡⁡kη2-isinh⁡kη2k·∇η2dx.We simplify (32) by noting that (33)eikxkcosh⁡⁡kη1-h-isinh⁡⁡kη1-hk·∇η1=-i∇·kkcosh⁡⁡kη1-heikx,eikxkcosh⁡⁡kη2-isinh⁡kη2k·∇η2=-i∇·kkcosh⁡(kη2)eikx.Using these identities in (32) and integrating by parts, we finally get the nonlocal spectral equation (2).Similarly, let us consider the basic functions(34)Ψk(x,y)=eikxsinh⁡(ky)and put Ψk into (20), obtaining(35)∫R2‍eikxsinh⁡⁡kη2+hη2tdx-∫R2‍eikxsinh⁡⁡kη1η1tdx=-∫R2‍eikxQkcosh⁡⁡kη1-isinh⁡(kη1)(k·∇)η1dx+∫R2‍eikxPkcosh⁡⁡kη2+h-isinh⁡kη2+h(k·∇)η2dx.We simplify (35) by noting that (36)eikxkcosh⁡⁡kη1-isinh⁡kη1k·∇η1=-i∇·kkcosh⁡⁡kη1eikx,eikxkcosh⁡⁡kη2+h-isinh⁡kη2+hk·∇η2=-i∇·kkcosh⁡kη2+heikx.Using these identities in (35) and integrating by parts, we finally get the nonlocal spectral equation (3).For the top layer, we require that (21) holds for the parametrized family of functions (37)Ψ~k(x,y)=eikxcosh⁡ky-H+h.Substituting Ψ~k(x,y) into (21), we obtain (38)∫R2‍eikxcosh⁡⁡kη2-Hη2tdx=∫R2‍Q~eikxksinh⁡⁡kη2-H-icosh⁡kη2-Hk·∇η2eikxksinh⁡⁡kη2-Hdx.We simplify (38) by noting that (39)keikxsinh⁡⁡kη2-H-ieikxcosh⁡⁡kη2-Hk·∇η2=-i∇·kksinh⁡kη2-Heikx.Using this identity in (28) and integrating by parts, we finally get the nonlocal spectral equation (4).


### 2.3. Conservation Laws and Integral Identities

In [17], conservation laws and integral identities for a two-fluid system were derived from the nonlocal spectral formulation. We now derive the analogous conservation laws for the three-fluid system for the NSP formulation. We start by expanding (2), (3) for small k. At first order, by setting to zero the coefficient of k and putting k2=0, we get (40)∫R2‍(η2η2t-η1η1t+Hη1t)dx=∫R2‍x1(Qx1-Px1)dx,(41)∫R2‍(η2η2t-η1η1t+Hη2t)dx=∫R2‍x1(Qx1-Px1)dx.Subtracting (41) from (40), we obtain the result(42)∂t∫R2‍(η1-η2)dx=0,which corresponds to mass conservation in the intermediate fluid domain.Similarly, by setting to zero the coefficient ofk2, with k2=0, we get(43)∂t∫R2‍x1η222H-η122H+η1dx=∫R2‍η222HPx1-η122H+H2-η1Qx1dx,(44)∂t∫R2‍x1η222H-η122H+η2dx=∫R2‍η222H+H2+η2Px1-η122HQx1dx.When we subtract (43) from (44), we obtain(45)∂t∫R2‍x1η2-η1dx=∫R2‍H2+η2Px1+H2-η1Qx1dx;similarly, if we set k1=0, at the same order we get the corresponding equation:(46)∂t∫R2‍x1η2-η1dx=∫R2‍H2+η2Px2+H2-η1Qx2dx.Equations (45) and (46) describe the evolution of the center of mass; the right-hand side is the momentum of the fluid.Next, we obtain a last identity at third order, setting to zero the coefficient ofk3:(47)∂t∫R2‍xj2(η2-η1)-H4(η22+η12)-16(η23+η13)dx=∫R2‍xjH2+η2Pxj+H2-η1Qxjdx,with j=1,2; the above is a virial type formula, analogous to the one obtained in [17].


## 2.1. A Weak Formulation of the Classic Three-Fluid Equations

We recall the classic equations governing three ideal fluids separated by two free interfacesη1 and η2 and bounded above and below by rigid lids. It is assumed that the lower fluid is of density ρ1, the intermediate fluid is of density ρ2, and the upper fluid is of density ρ3 with ρ1>ρ2>ρ3. The equations are given in terms of the interface variables and the velocity potentials φ, Φ, and Φ~ associated with the lower, intermediate, and upper fluid domains, respectively, (6)Δφ=0,for-h1<y<η1,(7)φyy=-h1=0,(8)η1t=(φy-∇xφ·∇xη1)y=η1,(9)ΔΦ=0,forη1<y<h+η2,(10)η1t=(Φy-∇xΦ·∇xη1)y=η1,(11)η2t=(Φy-∇xΦ·∇xη2)y=h+η2,(12)ΔΦ~=0,forh+η2<y<H+h,(13)Φ~yy=H+h=0,(14)η2t=(Φ~y-∇xΦ~·∇xη2)y=h+η2,(15)ρ1φt+12∇φ2+gη1-ρ2Φt+12∇Φ2+gη1=σ1∇·∇η11+∇η12,on⁡y=η1,(16)ρ2Φt+12∇Φ2+gη2-ρ3Φ~t+12∇Φ~2+gη2=σ2∇·∇η21+∇η22,on⁡y=h+η2.We also require that |∇φ|→0, |∇Φ|→0 and |∇Φ~|→0 as |x|→∞. In the previous equations, the constants g, σ1, and σ2 denote gravity, the surface tension associated with the free interface η1, and the surface tension associated with the free surface η2. Equations (6), (9), and (12) express the fact that fluids are divergence and curl-free. Equations (15)-(16) express that the jump in pressure across an interface is balanced by surface tension.We now obtain a weak formulation of (6)–(16), expressed in terms of η1, η2, q, Q, P, and Q~. Let us define the domains (17)D(η1)≡{(x,y)∈R2∣x∈R2,-h1<y<η1(x,t)},Dη1,η2≡(x,y)∈R2∣x∈R2,η1(x,t)<y<h+η2(x,t)(x,y)∈R2∣x∈R2,,Dη2≡(x,y)∈R2∣x∈R2,h+η2(x,t)<y<h+H(x,y)∈R2∣x∈R2.For any bounded, harmonic functions ψ,Ψ,Ψ~, defined in D(η1), D(η1,η2), and D(η2), respectively, and satisfying (18)Δψ=0,for(x,y)∈D(η1),ψyy=-h1=0,ΔΨ~=0,for(x,y)∈D(η2),Ψ~yy=H+h=0,∇ψ,∇Ψ,∇Ψ~⟶0,asx⟶∞,the following identities hold at any given time t: (19)∫R2‍ψy=η1η1tdx-∫R2‍q(ψy-∇xψ·∇xη1)y=η1dx=0,(20)∫R2‍Ψy=h+η2η2tdx-∫R2‍Ψy=η1η1tdx+∫R2‍QΨy-∇xΨ·∇xη1y=η1dx-∫R2‍P(Ψy-∇xΨ·∇xη2)y=h+η2dx=0,(21)∫R2‍Ψ~y=h+η2η2tdx-∫R2‍Q~Ψ~y-∇xΨ~·∇xη2y=h+η2dx=0.Equations (19) and (21), together with the Bernoulli equations (15) and (16), constitute a nonlocal system of equations describing the three-fluid system.


## 2.2. Derivation of the Nonlocal Spectral Equations

We derive from (19) the NSP equation (1) given in the Introduction. To motivate the derivation, suppose that ψ(x,y) satisfies(22)Δψ(x,y)=0,ψy(x,-h1)=0.Define ψ^(k,y) by (23)ψ^(k,y)=12π2∫R2‍e-ikxψ(x,y)dx.Then, conditions (22) turn into the following: (24)ψ^yy(k,y)=k2ψ^(k,y),ψ^y(k,-y)=0.The solution to the above ODE is given by (25)ψ^(k,y)=ξ^kcosh⁡ky+h1,reverting back to physical space, (26)ψ(x,y)=∫R2‍eikxξ^kcosh⁡ky+h1dk.The previous equation shows that any harmonic function ψ whose y-derivative at y=-h1 vanishes can formally be written as a sum of functions eikxcosh⁡(|k|(y+h1)). Therefore, by linearity it suffices to require that (19) holds for the parametrized family of functions (27)ψk(x,y)=eikxcosh⁡ky+h1.Substituting ψk(x,y) into (19), we obtain (28)∫R2‍eikxcosh⁡kη1+h1η1tdx=∫R2‍qeikx-icosh⁡kη1+h1(k·∇)η1ksinh⁡⁡kη1+h1-icosh⁡kη1+h1(k·∇)η1eikx-icosh⁡kη1+h1(k·∇)η1ksinh⁡⁡kη1+h1dx.We simplify (28) by noting that(29)keikxsinh⁡kη1+h1-ieikxcosh⁡⁡kη1+h1k·∇η1=-i∇·kksinh⁡kη1+h1eikx.Using this identity in (28) and integrating by parts, we finally get the nonlocal spectral equation (1) as follows: (30)∫R2‍eikxcosh⁡⁡kη1+h1η1tdx=i∫R2‍eikxsinh⁡kη1+h1k(k·∇)qdx.Along the same lines used for the first layer, let us derive (2)-(3) from (20).We define the basic functions(31)Ψk(x,y)=eikxsinh⁡ky-h,and putting Ψk into (20) gives(32)∫R2‍eikxsinh⁡⁡kη2η2tdx-∫R2‍eikxsinh⁡kη1-hη1tdx=-∫R2‍eikxQkcosh⁡⁡kη1-h-isinh⁡kη1-hk·∇η1dx+∫R2‍eikxPkcosh⁡⁡kη2-isinh⁡kη2k·∇η2dx.We simplify (32) by noting that (33)eikxkcosh⁡⁡kη1-h-isinh⁡⁡kη1-hk·∇η1=-i∇·kkcosh⁡⁡kη1-heikx,eikxkcosh⁡⁡kη2-isinh⁡kη2k·∇η2=-i∇·kkcosh⁡(kη2)eikx.Using these identities in (32) and integrating by parts, we finally get the nonlocal spectral equation (2).Similarly, let us consider the basic functions(34)Ψk(x,y)=eikxsinh⁡(ky)and put Ψk into (20), obtaining(35)∫R2‍eikxsinh⁡⁡kη2+hη2tdx-∫R2‍eikxsinh⁡⁡kη1η1tdx=-∫R2‍eikxQkcosh⁡⁡kη1-isinh⁡(kη1)(k·∇)η1dx+∫R2‍eikxPkcosh⁡⁡kη2+h-isinh⁡kη2+h(k·∇)η2dx.We simplify (35) by noting that (36)eikxkcosh⁡⁡kη1-isinh⁡kη1k·∇η1=-i∇·kkcosh⁡⁡kη1eikx,eikxkcosh⁡⁡kη2+h-isinh⁡kη2+hk·∇η2=-i∇·kkcosh⁡kη2+heikx.Using these identities in (35) and integrating by parts, we finally get the nonlocal spectral equation (3).For the top layer, we require that (21) holds for the parametrized family of functions (37)Ψ~k(x,y)=eikxcosh⁡ky-H+h.Substituting Ψ~k(x,y) into (21), we obtain (38)∫R2‍eikxcosh⁡⁡kη2-Hη2tdx=∫R2‍Q~eikxksinh⁡⁡kη2-H-icosh⁡kη2-Hk·∇η2eikxksinh⁡⁡kη2-Hdx.We simplify (38) by noting that (39)keikxsinh⁡⁡kη2-H-ieikxcosh⁡⁡kη2-Hk·∇η2=-i∇·kksinh⁡kη2-Heikx.Using this identity in (28) and integrating by parts, we finally get the nonlocal spectral equation (4).


## 2.3. Conservation Laws and Integral Identities

In [17], conservation laws and integral identities for a two-fluid system were derived from the nonlocal spectral formulation. We now derive the analogous conservation laws for the three-fluid system for the NSP formulation. We start by expanding (2), (3) for small k. At first order, by setting to zero the coefficient of k and putting k2=0, we get (40)∫R2‍(η2η2t-η1η1t+Hη1t)dx=∫R2‍x1(Qx1-Px1)dx,(41)∫R2‍(η2η2t-η1η1t+Hη2t)dx=∫R2‍x1(Qx1-Px1)dx.Subtracting (41) from (40), we obtain the result(42)∂t∫R2‍(η1-η2)dx=0,which corresponds to mass conservation in the intermediate fluid domain.Similarly, by setting to zero the coefficient ofk2, with k2=0, we get(43)∂t∫R2‍x1η222H-η122H+η1dx=∫R2‍η222HPx1-η122H+H2-η1Qx1dx,(44)∂t∫R2‍x1η222H-η122H+η2dx=∫R2‍η222H+H2+η2Px1-η122HQx1dx.When we subtract (43) from (44), we obtain(45)∂t∫R2‍x1η2-η1dx=∫R2‍H2+η2Px1+H2-η1Qx1dx;similarly, if we set k1=0, at the same order we get the corresponding equation:(46)∂t∫R2‍x1η2-η1dx=∫R2‍H2+η2Px2+H2-η1Qx2dx.Equations (45) and (46) describe the evolution of the center of mass; the right-hand side is the momentum of the fluid.Next, we obtain a last identity at third order, setting to zero the coefficient ofk3:(47)∂t∫R2‍xj2(η2-η1)-H4(η22+η12)-16(η23+η13)dx=∫R2‍xjH2+η2Pxj+H2-η1Qxjdx,with j=1,2; the above is a virial type formula, analogous to the one obtained in [17].


## 3. Weakly Nonlinear Equations

### 3.1. Nondimensionalization of the NSP Equations

In order to derive weakly nonlinear equations, we nondimensionalize all physical variables in (1)–(5) according to(48)x1′=x1L,x2′=γx2L,t′=c0Lt,η1,2=aη1,2′,β=aβ′,q=agLc0q′,Q=agHc0Q′,P=agHc0P′,Q~=agLc0Q~′,where L is a characteristic wavelength, γ is a nondimensional parameter, and a and c0=gH are a characteristic amplitude and velocity, respectively.Then, (1)–(5) become, after dropping primes and letting k→-k,first layer (bottom)(49)∫R2‍e-ikxcosh⁡kα1+εμη1η1tdx=iμ∫R2‍e-ikxsinh⁡kα1+εμη1k(k·∇)qdx,intermediate layer(50)∫R2‍e-ikxsinh⁡⁡kεμη2η2tdx-∫R2‍e-ikxsinh⁡⁡kεμη1-α2η1tdx=-i∫R2‍e-ikxcosh⁡⁡kεμη1-α2kk·∇Qdx+i∫R2‍e-ikxcosh⁡⁡kεμη2kk·∇Pdx,(51)∫R2‍e-ikxsinh⁡⁡kα2+εμη2η2tdx-∫R2‍e-ikxsinh⁡⁡kεμη1η1tdx=-i∫R2‍e-ikxcosh⁡⁡kεμη1kk·∇Qdx+i∫R2‍e-ikxcosh⁡kα2+εμη2k(k·∇)Pdx,third layer (top)(52)∫R2‍e-ikxcosh⁡(k(εμη2-μ))η2tdx=iμ∫R2‍e-ikxsinh⁡kεμη2-μk(k·∇)Q~dx,Bernoulli’s equations(53)qt+12ε∇q2+η1-εμ2η1t+ε∇q∇η122(1+ε2μ2∇η12)-ρ~η1t+εμ∇Q∇η1221+ε2μ2∇η12μQt+12εμ2∇Q2+η1-εμ2η1t+εμ∇Q∇η1221+ε2μ2∇η12=μ2σ~1∇∇η11+ε2μ2∇η12,(54)μPt+12εμ2∇P2+η2-εμ2η2t+εμ∇P∇η222(1+ε2μ2∇η22)-ρQ~t+12ε∇Q~2+η2-εμ2η2t+ε∇Q~∇η222(1+ε2μ2∇η22)=μ2σ~2∇∇η21+ε2μ2∇η22.In the above relations, ε=a/H is thenonlinearity ratio and μ=H/L is the so-calledaspect ratio [12]; moreover, it is (55)α1=h1L,α2=hL,σ~1=σ1H2gρ1,σ~2=σ2H2gρ2,ρ~=ρ2ρ1,ρ=ρ3ρ2.In the following, we assume αj≪1, j=1,2 and μ≪1, which corresponds to the case oflong waves; moreover, we take the nonlinearity ratio ε to be ε=O(μ).The above assumptions imply that we are interested to derive asymptotic reductions of the NSP equations in the case of weakly nonlinear long waves.


### 3.2. Derivation of the Equations

We first expand (49) and (53) in μ and γ2. From (49), taking the Inverse Fourier Transform (IFT), we get(56)qx1=iμcoth⁡(α1D1)η1t+o(μ2),while (53) gives(57)η1=-qt1-ρ~-ε21-ρ~qx1+γ2qx22+μρ~(1-ρ~)Qt+μ2σ~1η1x1x1+o(εμ2).On the other hand, from (49), it is also useful to derive the following reduction: (58)η1t=-α1μqx1x1+γ2qx2x2-ε∂x1η1qx1+oεμ2.We now use (58) at leading order (L.O.) in (56) and, after some manipulations, we obtain from (57) the following equation:(59)η1t=iα1coth⁡⁡α1D11-ρ~qx1tt-ε21-ρ~∂tqx1+γ2qx22+μρ~(1-ρ~)Qtt+o(μ2).By (58) and (59), thefirst generalized Boussinesq equation obtains (60)α1qx1x1+γ2qx2x2+iμα1coth⁡⁡α1D11-ρ~qx1tt-εμ1-ρ~∂x1qtqx1-εμ21-ρ~∂tqx12+μ2ρ~1-ρ~Qtt+o(εμ2,μ3)=0.We now turn our attention to the intermediate layer.By equating the expansions of (50) and (51), we get the following expression for Qx1: (61)Qx1=-iη2tsinh⁡(α2D1)+icoth⁡(α2D1)η1t+o(εμ)which in turn gives (62)Qx1t=-iη2ttsinh⁡⁡α2D1+icoth⁡α2D1η1tt+oεμ,where, from (58) and (69) at L.O., we have (63)η1tt=-α1μqx1x1t,η2tt=Q~x1x1t.When (63) are used in (62), after integrating with respect to x1, we finally get (64)Qt=-iQ~x1tsinh⁡(α2D1)-iα1μcoth⁡(α2D1)qx1t+o(εμ).We now take the t-derivative of (57) and use (64), obtaining (65)η1t=-qtt1-ρ~+iα1ρ~1-ρ~coth⁡⁡α2D1qx1tt-ε21-ρ~∂tqx1+γ2qx22-iμρ~1-ρ~Q~x1ttsinh⁡⁡α2D1+o(εμ2);by equating (65) and (58), we get thesecond generalized Boussinesq equation: (66)α1qx1x1+γ2qx2x2-μ1-ρ~qtt-εμ1-ρ~∂x1(qtqx1)-εμ21-ρ~∂tqx12-iα1μρ~1-ρ~coth⁡⁡α2D1qx1tt-iμ2ρ~(1-ρ~)Q~x1ttsinh⁡(α2D1)+o(εμ2)=0.We now expand (52) and take the IFT, getting(67)Q~x1=-iμcoth⁡(μD1)η2t+o(μ2),while from the expansion of (54) we get (68)η2=ρ1-ρQ~t-μ1-ρPt+ε2ρ1-ρQ~x1+γ2Q~x22+oεμ2.On the other hand, from (52), it is also useful to derive the following expansion: (69)η2t=Q~x1x1+γ2Q~x2x2-ε∂x1η2Q~x1+oεμ.We use (68) at L.O. in (69) and get (70)η2t=Q~x1x1+γ2Q~x2x2-ερ(1-ρ)∂x1(Q~tQ~x1)+o(εμ).In order to get an expression for Pt in (68), we now expand (50), getting (71)Px1=isinh⁡α2D1η1t+cosh⁡α2D1Qx1+oεμ.By taking the t-derivative of (71) and integrating with respect to x1, we get (72)Pt=-iα1μsinh⁡(α2D1)qx1t+cosh⁡(α2D1)Qt+o(εμ),where the L.O. of (58) has also been used.When we substitute (72) back into (68) and take the t-derivative, we obtain (73)η2t=iα11-ρsinh⁡⁡α2D1qx1tt+ρ1-ρQ~tt-μ1-ρcosh⁡⁡α2D1Qtt+ε2ρ(1-ρ)∂tQ~x12+o(εμ2)=0.By equating (70) and (73), we finally obtain thethird generalized Boussinesq equation:(74)Q~tt-1-ρρQ~x1x1+γ2Q~x2x2+iα1ρsinh⁡⁡α2D1qx1tt+ε2∂tQ~x12-μρcosh⁡⁡(α2D1)Qtt+ε∂x1(Q~tQ~x1)+o(εμ)=0.The system of 3 generalized Boussinesq equations (60), (66), and (74) will now be reduced to a system of 2 independent equations. We first take the shallow water limit, given by α1≪1 in (60) and α2≪1 in (66), obtaining the following reductions: (75)qx1x1+γ2qx2x2-μα111-ρ~qtt-εμα11(1-ρ~)∂x1(qtqx1)-εμα1121-ρ~∂tqx12+μ2α1ρ~(1-ρ~)Qtt+o(εμ2)=0,(76)qx1x1+γ2qx2x2-μα111-ρ~qtt-εμα11(1-ρ~)∂x1(qtqx1)-εμα1121-ρ~∂tqx12-μα2ρ~1-ρ~qtt-μ2α2α1ρ~(1-ρ~)Q~tt+o(εμ2)=0.Next, the combination of (75) and (76) gives (77)Qtt=-α1α21μqtt-1α2Q~tt,which we substitute back into (74) and we obtain (78)Q~tt-1-ρρQ~x1x1+γ2Q~x2x2=α1α2ρqx1x1tt-ε2∂tQ~x12-α1α21ρqtt-μα21ρQ~tt-ε∂x1(Q~tQ~x1).Rearranging (66), we get (79)qtt-α1μ1-ρ~qx1x1+γ2qx2x2=α1α2ρ~qtt-ε∂x1(qtqx1)-ε2∂tqx12+μα2ρ~Q~tt;then the system of weakly nonlinear, shallow water equations is given by (78) and (79).


## 3.1. Nondimensionalization of the NSP Equations

In order to derive weakly nonlinear equations, we nondimensionalize all physical variables in (1)–(5) according to(48)x1′=x1L,x2′=γx2L,t′=c0Lt,η1,2=aη1,2′,β=aβ′,q=agLc0q′,Q=agHc0Q′,P=agHc0P′,Q~=agLc0Q~′,where L is a characteristic wavelength, γ is a nondimensional parameter, and a and c0=gH are a characteristic amplitude and velocity, respectively.Then, (1)–(5) become, after dropping primes and letting k→-k,first layer (bottom)(49)∫R2‍e-ikxcosh⁡kα1+εμη1η1tdx=iμ∫R2‍e-ikxsinh⁡kα1+εμη1k(k·∇)qdx,intermediate layer(50)∫R2‍e-ikxsinh⁡⁡kεμη2η2tdx-∫R2‍e-ikxsinh⁡⁡kεμη1-α2η1tdx=-i∫R2‍e-ikxcosh⁡⁡kεμη1-α2kk·∇Qdx+i∫R2‍e-ikxcosh⁡⁡kεμη2kk·∇Pdx,(51)∫R2‍e-ikxsinh⁡⁡kα2+εμη2η2tdx-∫R2‍e-ikxsinh⁡⁡kεμη1η1tdx=-i∫R2‍e-ikxcosh⁡⁡kεμη1kk·∇Qdx+i∫R2‍e-ikxcosh⁡kα2+εμη2k(k·∇)Pdx,third layer (top)(52)∫R2‍e-ikxcosh⁡(k(εμη2-μ))η2tdx=iμ∫R2‍e-ikxsinh⁡kεμη2-μk(k·∇)Q~dx,Bernoulli’s equations(53)qt+12ε∇q2+η1-εμ2η1t+ε∇q∇η122(1+ε2μ2∇η12)-ρ~η1t+εμ∇Q∇η1221+ε2μ2∇η12μQt+12εμ2∇Q2+η1-εμ2η1t+εμ∇Q∇η1221+ε2μ2∇η12=μ2σ~1∇∇η11+ε2μ2∇η12,(54)μPt+12εμ2∇P2+η2-εμ2η2t+εμ∇P∇η222(1+ε2μ2∇η22)-ρQ~t+12ε∇Q~2+η2-εμ2η2t+ε∇Q~∇η222(1+ε2μ2∇η22)=μ2σ~2∇∇η21+ε2μ2∇η22.In the above relations, ε=a/H is thenonlinearity ratio and μ=H/L is the so-calledaspect ratio [12]; moreover, it is (55)α1=h1L,α2=hL,σ~1=σ1H2gρ1,σ~2=σ2H2gρ2,ρ~=ρ2ρ1,ρ=ρ3ρ2.In the following, we assume αj≪1, j=1,2 and μ≪1, which corresponds to the case oflong waves; moreover, we take the nonlinearity ratio ε to be ε=O(μ).The above assumptions imply that we are interested to derive asymptotic reductions of the NSP equations in the case of weakly nonlinear long waves.


## 3.2. Derivation of the Equations

We first expand (49) and (53) in μ and γ2. From (49), taking the Inverse Fourier Transform (IFT), we get(56)qx1=iμcoth⁡(α1D1)η1t+o(μ2),while (53) gives(57)η1=-qt1-ρ~-ε21-ρ~qx1+γ2qx22+μρ~(1-ρ~)Qt+μ2σ~1η1x1x1+o(εμ2).On the other hand, from (49), it is also useful to derive the following reduction: (58)η1t=-α1μqx1x1+γ2qx2x2-ε∂x1η1qx1+oεμ2.We now use (58) at leading order (L.O.) in (56) and, after some manipulations, we obtain from (57) the following equation:(59)η1t=iα1coth⁡⁡α1D11-ρ~qx1tt-ε21-ρ~∂tqx1+γ2qx22+μρ~(1-ρ~)Qtt+o(μ2).By (58) and (59), thefirst generalized Boussinesq equation obtains (60)α1qx1x1+γ2qx2x2+iμα1coth⁡⁡α1D11-ρ~qx1tt-εμ1-ρ~∂x1qtqx1-εμ21-ρ~∂tqx12+μ2ρ~1-ρ~Qtt+o(εμ2,μ3)=0.We now turn our attention to the intermediate layer.By equating the expansions of (50) and (51), we get the following expression for Qx1: (61)Qx1=-iη2tsinh⁡(α2D1)+icoth⁡(α2D1)η1t+o(εμ)which in turn gives (62)Qx1t=-iη2ttsinh⁡⁡α2D1+icoth⁡α2D1η1tt+oεμ,where, from (58) and (69) at L.O., we have (63)η1tt=-α1μqx1x1t,η2tt=Q~x1x1t.When (63) are used in (62), after integrating with respect to x1, we finally get (64)Qt=-iQ~x1tsinh⁡(α2D1)-iα1μcoth⁡(α2D1)qx1t+o(εμ).We now take the t-derivative of (57) and use (64), obtaining (65)η1t=-qtt1-ρ~+iα1ρ~1-ρ~coth⁡⁡α2D1qx1tt-ε21-ρ~∂tqx1+γ2qx22-iμρ~1-ρ~Q~x1ttsinh⁡⁡α2D1+o(εμ2);by equating (65) and (58), we get thesecond generalized Boussinesq equation: (66)α1qx1x1+γ2qx2x2-μ1-ρ~qtt-εμ1-ρ~∂x1(qtqx1)-εμ21-ρ~∂tqx12-iα1μρ~1-ρ~coth⁡⁡α2D1qx1tt-iμ2ρ~(1-ρ~)Q~x1ttsinh⁡(α2D1)+o(εμ2)=0.We now expand (52) and take the IFT, getting(67)Q~x1=-iμcoth⁡(μD1)η2t+o(μ2),while from the expansion of (54) we get (68)η2=ρ1-ρQ~t-μ1-ρPt+ε2ρ1-ρQ~x1+γ2Q~x22+oεμ2.On the other hand, from (52), it is also useful to derive the following expansion: (69)η2t=Q~x1x1+γ2Q~x2x2-ε∂x1η2Q~x1+oεμ.We use (68) at L.O. in (69) and get (70)η2t=Q~x1x1+γ2Q~x2x2-ερ(1-ρ)∂x1(Q~tQ~x1)+o(εμ).In order to get an expression for Pt in (68), we now expand (50), getting (71)Px1=isinh⁡α2D1η1t+cosh⁡α2D1Qx1+oεμ.By taking the t-derivative of (71) and integrating with respect to x1, we get (72)Pt=-iα1μsinh⁡(α2D1)qx1t+cosh⁡(α2D1)Qt+o(εμ),where the L.O. of (58) has also been used.When we substitute (72) back into (68) and take the t-derivative, we obtain (73)η2t=iα11-ρsinh⁡⁡α2D1qx1tt+ρ1-ρQ~tt-μ1-ρcosh⁡⁡α2D1Qtt+ε2ρ(1-ρ)∂tQ~x12+o(εμ2)=0.By equating (70) and (73), we finally obtain thethird generalized Boussinesq equation:(74)Q~tt-1-ρρQ~x1x1+γ2Q~x2x2+iα1ρsinh⁡⁡α2D1qx1tt+ε2∂tQ~x12-μρcosh⁡⁡(α2D1)Qtt+ε∂x1(Q~tQ~x1)+o(εμ)=0.The system of 3 generalized Boussinesq equations (60), (66), and (74) will now be reduced to a system of 2 independent equations. We first take the shallow water limit, given by α1≪1 in (60) and α2≪1 in (66), obtaining the following reductions: (75)qx1x1+γ2qx2x2-μα111-ρ~qtt-εμα11(1-ρ~)∂x1(qtqx1)-εμα1121-ρ~∂tqx12+μ2α1ρ~(1-ρ~)Qtt+o(εμ2)=0,(76)qx1x1+γ2qx2x2-μα111-ρ~qtt-εμα11(1-ρ~)∂x1(qtqx1)-εμα1121-ρ~∂tqx12-μα2ρ~1-ρ~qtt-μ2α2α1ρ~(1-ρ~)Q~tt+o(εμ2)=0.Next, the combination of (75) and (76) gives (77)Qtt=-α1α21μqtt-1α2Q~tt,which we substitute back into (74) and we obtain (78)Q~tt-1-ρρQ~x1x1+γ2Q~x2x2=α1α2ρqx1x1tt-ε2∂tQ~x12-α1α21ρqtt-μα21ρQ~tt-ε∂x1(Q~tQ~x1).Rearranging (66), we get (79)qtt-α1μ1-ρ~qx1x1+γ2qx2x2=α1α2ρ~qtt-ε∂x1(qtqx1)-ε2∂tqx12+μα2ρ~Q~tt;then the system of weakly nonlinear, shallow water equations is given by (78) and (79).


## 4. Multiple Scale Derivation: Special Case

In order to obtain some interesting limiting equations, we now make the assumption ofmaximal balance (the small terms are of the same order) ε=μ=γ2. This reflects a balance of weak nonlinearity and weak dispersion.In this section, we study aSpecial Case, that is, when the waves velocities in the shallow water equations (78) and (79) are the same: (80)α1μ(1-ρ~)=1-ρρ=:co2.Here below, we assume the following asymptotic expansions forq, Q, and Q~: (81)q=q0+εq1+⋯,(82)Q=Q0+εQ1+⋯,(83)Q~=Q~0+εQ~1+⋯.We also introduce new variables: (84)ξ=x1-c0t,ζ=x1+c0t,T=εt,y=x2,which aretravelling waves variables and describe the direction of waves propagation along the positive x-axis: ξ represents right moving waves and ζ represents left moving waves.In terms of the new variables, one obtains(85)∂t=-c0∂ξ+c0∂ζ+ε∂T;∂x1=∂ξ+∂ζ.Considering the assumption of maximal balance together with (80), (79) becomes (86)qtt-co2qx1x1-ρ~α2εQ~tt+qtt+ε2∂t(qx12)-εc02qx2x2+ε∂x1(qtqx1)=0.Substituting the expansions (81), (83) into (86) and equating L.O. terms yield the wave equations: q0tt-q0xx=0 and Q~0tt-Q~0xx=0, whose solutions, respectively, are q0=F(ξ,T)+G(ζ,T) and Q~0=F~(ξ,T)+G~(ζ,T).We also assume unidirectional waves andonly work with the right moving waves, so that (87)q0=F(ξ,T),Q0=F^(ξ,T),Q~0=F~(ξ,T).Then, keeping the o(ε) terms, (86) becomes (88)-4c02q1ξζ=2c0FξT+ρ~c02α2F~ξξ+Fξξ+3c0FξFξξ+c02Fyy.In order to remove secular terms, the right-hand side of (88) needs to be set to zero (see [18]): (89)2c0FξT+ρ~c02α2F~ξξ+Fξξ+3c0FξFξξ+c02Fyy=0.Following the same steps as above for (78), together with (89), and studying only the (1+1)-dimensional case, we obtain a set of two equations:(90)2c0FξT+ρ~c02α2F~ξξ+Fξξ+3c0FξFξξ=0,2c0F~ξT-1ρρ~2c0FξT+c0FξFξξ+α2c02ρFξξξξ+3c0F~ξF~ξξ=0.In the following, we will study only the(1+1)-dimensional case of the derivation.Finally, define the functionsU=Fξ, U~=F~ξ; then, differentiating (90) with respect to ξ, we can write the resulting equations in terms of U and U~ as (91)2UξT+c0ρ~α2U~ξξ+Uξξ+32∂ξ2(U2)=0,(92)2U~ξT-1ρρ~2UξT+12∂ξ2(U2)+α2c0ρUξξξξ+32∂ξ2(U~2)=0.Substitute (91) into (92) and obtain (93)2ρρ~U~ξT+c0ρ~α2U~ξξ+c0ρ~α2Uξξ+∂ξ2(U2)+α2ρ~c0Uξξξξ+32ρρ~∂ξ2(U~2)=0.Our multiple scale derivation leads therefore to having a system of two equations: (91) and (93) constitute a system of coupled nonlinearshallow water equations.The two coupled shallow water equations will be studied numerically as a function of the parameters entering the theory, in order to prove the existence of solitary waves and analyze their behaviour.


### 4.1. Numerical Investigation of Traveling Wave Solutions of the Coupled Shallow Water Equations

We now investigate whether (91) and (93) possess traveling wave solutions. For convenience, take c0=1. Passing to a traveling coordinate system moving with velocity c in the ξ direction, (91) and (93) become, respectively, (94)-2cUξξ+ρ~α2U~ξξ+Uξξ+32∂ξ2(U2)=0,-2cρρ~U~ξξ+ρ~α2U~ξξ+Uξξ+α2ρ~Uξξξξ+∂ξ2(U2)+32ρρ~∂ξ2(U~2)=0.To solve (94) numerically, we use the spectral renormalization (SPRZ) method, developed in [19]. We take the Fourier transform of the previous equations and get, upon rearranging, (95)U^=2cρρ~U~^+1/2U2^-3/2ρρ~U~2^2c+ρ~k12α2,(96)U~^=-2cα2U^-ρ~U^+3/2α2U2^ρ~.In general, we cannot find a solution to (95)-(96) by naive iteration. Instead, we assume that U=λv and U~=λ~v~, where λ, λ~ are unknown parameters and v, v~ are unknown functions (this step is the renormalization part). Then, (95) can be written in terms of λ and v as (97)v^=1λ2cρρ~(λ~v~^)+1/2λ2v2^-3/2ρρ~(λ~2v~2^)2c+ρ~k12α2.Note that by multiplying (97) by v^¯ (where v^¯ denotes the conjugate of v^), rearranging and integrating the result we get (98)λ=∫R‍[2cρρ~(λ~v~^)+1/2λ2v2^-3/2ρρ~(λ~2v~2^)]v^¯dk1∫R‍[2c+ρ~k12α2]v^v^¯dk1.Along the same line, (96) can be written in terms of λ~ and v~ as (99)v~^=1λ~-2cα2(λv^)-ρ~(λv^)+3/2α2(λ2v2^)ρ~.By multiplying (99) by v~^¯ (where v~^¯ denotes the conjugate of v~^), rearranging and integrating the result we get(100)λ~=∫R‍[-2cα2(λv^)-ρ~(λv^)+3/2α2(λ2v2^)]v~^¯dk1∫R‍ρ~v~^v~^¯dk1.We use the above SPRZ scheme to solve for the modesU and U~ when α2=0.2, c=-1, ρ=0.5, and ρ~=0.7.Figure2 displays the resulting speed versus amplitude relationship, which is nearly linear for α2=0.2 (shallow water regime). Note that the horizontal axis in Figure 2 is -c.Figure 2
Amplitude versus negative speed forα2=0.2, μ=0.1.Figure3 shows the x cross sections of U and U~; on the left there are the solutions starting from a Gaussian function and on the right the ones starting from a step function, both for a typical speed c=-1 in the shallow water regime.Figure 3
x cross sections of U and U~ for α2=0.2, μ=0.05, c=-1, ρ=0.5, and ρ~=0.7.After obtaining these results, it would now be of particular interest from the applied point of view to address the issue of a local breakdown of the nonlinear internal waves propagating in the density stratified fluid [20]. Phenomena of this kind can induce localized turbulence in the stratified fluid and are relevant for their occurrence in the atmosphere and the oceans. We plan to address this issue in the future together with a generalization of our model to the (2+1)-dimensional case.


## 4.1. Numerical Investigation of Traveling Wave Solutions of the Coupled Shallow Water Equations

We now investigate whether (91) and (93) possess traveling wave solutions. For convenience, take c0=1. Passing to a traveling coordinate system moving with velocity c in the ξ direction, (91) and (93) become, respectively, (94)-2cUξξ+ρ~α2U~ξξ+Uξξ+32∂ξ2(U2)=0,-2cρρ~U~ξξ+ρ~α2U~ξξ+Uξξ+α2ρ~Uξξξξ+∂ξ2(U2)+32ρρ~∂ξ2(U~2)=0.To solve (94) numerically, we use the spectral renormalization (SPRZ) method, developed in [19]. We take the Fourier transform of the previous equations and get, upon rearranging, (95)U^=2cρρ~U~^+1/2U2^-3/2ρρ~U~2^2c+ρ~k12α2,(96)U~^=-2cα2U^-ρ~U^+3/2α2U2^ρ~.In general, we cannot find a solution to (95)-(96) by naive iteration. Instead, we assume that U=λv and U~=λ~v~, where λ, λ~ are unknown parameters and v, v~ are unknown functions (this step is the renormalization part). Then, (95) can be written in terms of λ and v as (97)v^=1λ2cρρ~(λ~v~^)+1/2λ2v2^-3/2ρρ~(λ~2v~2^)2c+ρ~k12α2.Note that by multiplying (97) by v^¯ (where v^¯ denotes the conjugate of v^), rearranging and integrating the result we get (98)λ=∫R‍[2cρρ~(λ~v~^)+1/2λ2v2^-3/2ρρ~(λ~2v~2^)]v^¯dk1∫R‍[2c+ρ~k12α2]v^v^¯dk1.Along the same line, (96) can be written in terms of λ~ and v~ as (99)v~^=1λ~-2cα2(λv^)-ρ~(λv^)+3/2α2(λ2v2^)ρ~.By multiplying (99) by v~^¯ (where v~^¯ denotes the conjugate of v~^), rearranging and integrating the result we get(100)λ~=∫R‍[-2cα2(λv^)-ρ~(λv^)+3/2α2(λ2v2^)]v~^¯dk1∫R‍ρ~v~^v~^¯dk1.We use the above SPRZ scheme to solve for the modesU and U~ when α2=0.2, c=-1, ρ=0.5, and ρ~=0.7.Figure2 displays the resulting speed versus amplitude relationship, which is nearly linear for α2=0.2 (shallow water regime). Note that the horizontal axis in Figure 2 is -c.Figure 2
Amplitude versus negative speed forα2=0.2, μ=0.1.Figure3 shows the x cross sections of U and U~; on the left there are the solutions starting from a Gaussian function and on the right the ones starting from a step function, both for a typical speed c=-1 in the shallow water regime.Figure 3
x cross sections of U and U~ for α2=0.2, μ=0.05, c=-1, ρ=0.5, and ρ~=0.7.After obtaining these results, it would now be of particular interest from the applied point of view to address the issue of a local breakdown of the nonlinear internal waves propagating in the density stratified fluid [20]. Phenomena of this kind can induce localized turbulence in the stratified fluid and are relevant for their occurrence in the atmosphere and the oceans. We plan to address this issue in the future together with a generalization of our model to the (2+1)-dimensional case.



---

*Source: 197978-2015-09-17.xml*