( PDF ) A reduced model for internal waves interacting with submarine structures at intermediate depth

Download PDF

ads:

INSTITUTO NACIONAL DE MATEM

ATICA PURA E APLICADA

A REDUCED MODEL FOR

INTERNAL WAVES INTERACTING

WITH SUBMARINE STRUCTURES

AT INTERMEDIATE DEPTH

Author: Ail

ın Ruiz de Z

arate F

abregas

Adviser: Prof. Dr. Andr

e Nachbin

2007

ads:

Livros Grátis

http://www.livrosgratis.com.br

Milhares de livros grátis para download.

To my family

ads:

Acknowledgements

First, I would like to thank Professor Nachbin for his support, encouragement and

guidance along these years and for introducing me to the fascinating world of

water waves. It has been a privilege to work with him.

It has been also a pleasure to study at IMPA, for the institution provides all the

conditions for it. In particular, the Fluid Dynamics Group at IMPA has developed

an excellent atmosphere for work.

I appreciate the eﬀort of my teachers in Brazil and in Cuba. I owe them my

best results, now and ever.

I would like to thank Professor Wooyoung Choi for his useful comments re-

garding this work and future extensions and Professor Uri Ascher for his observa-

tions about the stability of numerical schemes.

I would like to express my deepest gratitude to my husband, Rodrigo Morante,

for his support, help and love along these years. Without him I would not have

been able to complete this project. Inparticular, I thank him for helping me typing,

generating part of the ﬁgures and correcting the work presented here.

I thank my Mom for have been always ready to help me, specially in the most

diﬃcult moments.

I thank my friends, they are few but very trustworthy.

I thank Professors Daniel Alfaro, Stefanella Boatto, Roberto Kraenkel, Dan

Marchesin and Jorge Zubelli for kindly accepting to participate in my thesis com-

mittee.

I thank the ANP/PRH-32 scholarship program and CAPES for ﬁnancial sup-

port.

Contents

Abstract 1

Resumo 2

1 Introduction 3

2 Derivation of the reduced model 8

2.1 Reducing the upper layer dynamics to the interface . . . . . . . . 11

2.2 Connecting the upper and lower layers . . . . . . . . . . . . . . . 18

2.3 Dispersion relation for the linearized model . . . . . . . . . . . . 27

2.4 Unidirectional wave regime . . . . . . . . . . . . . . . . . . . . . 29

2.5 Solitary wave solutions . . . . . . . . . . . . . . . . . . . . . . . 34

3 A higher-order reduced model 36

3.1 Higher-order upper layer equations . . . . . . . . . . . . . . . . . 36

3.2 Improved approximation for pressure at the interface . . . . . . . 41

3.3 Dispersion relation for the higher-order model. Comparison with

the previous model . . . . . . . . . . . . . . . . . . . . . . . . . 47

4 Numerical results 52

4.1 Hierarchy of one-dimensional models . . . . . . . . . . . . . . . 52

4.2 Method of lines . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4.3 Flat bottom experiments . . . . . . . . . . . . . . . . . . . . . . 66

4.4 Periodic topography experiments . . . . . . . . . . . . . . . . . . 75

4.5 Computing solitary waves solutions . . . . . . . . . . . . . . . . 81

Conclusions and future work 86

A Approximation for the horizontal derivatives at the unperturbed in-

terface 88

B The Dirichlet-to-Neumann operator 91

C The periodic counterpart of the operator T 93

Bibliography 98

Abstract

A reduced one-dimensional strongly nonlinear model for the evolution of inter-

nal waves over an arbitrary bottom topography is derived. The reduced model is

aimed at obtaining an eﬃcient numerical method for the two-dimensional prob-

lem. Two layers containing inviscid, immiscible, irrotational ﬂuids of diﬀerent

densities are deﬁned. The upper layer is shallow compared with the characteristic

wavelength at the interface of the two-ﬂuid system, while the depth of the bottom

region is comparable to the characteristic wavelength. The nonlinear evolution

equations obtained describe the behaviour of the internal wave elevation and mean

upper-velocity for this water conﬁguration. The system is a generalization of the

one proposed by Choi and Camassa for the ﬂat bottom case in the same physical

settings. Due to the presence of topography a variable coeﬃcient accompanies

each space derivative. These Boussinesq-type equations contain the Intermediate

Long Wave (ILW) equation and the Benjamin-Ono (BO) equation when restricted

to the unidirectional wave regime. We intend to use this model to study the inter-

action of waves with the bottom proﬁle. The dynamics include wave scattering,

dispersion and attenuation among other phenomena. The research is relevant in

oil recovery in deep ocean waters, where salt concentration and diﬀerences in

temperature generate stratiﬁcation in such a way that internal waves can aﬀect

oﬀshore operations and submerged structures.

Resumo

E obtido um modelo reduzido unidirecional fortemente n

ao linear para a evoluc¸

de ondas internas sobre topograﬁas de fundo arbitr

ario. Com o modelo reduzido

busca-se obter m

etodos num

ericos eﬁcientes para resolver o problema bidimen-

sional. S

ao consideradas duas camadas contendo dois ﬂuidos inv

ıscidos, imis-

ıveis e irrotacionais de densidades diferentes. A camada superior

e delgada

se comparada

a longitude de onda caracter

ıstica. As equac¸

oes de evoluc¸

ao n

lineares obtidas descrevem o comportamento da elevac¸

ao da onda interna e a

velocidade superior m

edia para esta conﬁgurac¸

ao da

agua. O sistema

e uma

generalizac¸

ao daquele proposto por Choi e Camassa para o caso de fundo plano

nas mesmas condic¸

oes f

ısicas. Devido

a presenc¸a da topograﬁa, cada derivada es-

pacial est

a acompanhada por um coeﬁciente vari

avel. Estas equac¸

oes de Boussi-

nesq cont

em a equac¸

ao da Onda Longa Intermedi

aria (Intermediate Long Wave,

ILW) e a equac¸

ao de Benjamin-Ono (BO) se restritas ao regime unidirecional de

propagac¸

ao de ondas. Pretendemos utilizar este modelo para estudar a interac¸

das ondascomo perﬁl do fundo. A din

amica inclui reﬂex

ao, dispers

ao e atenuac¸

das ondas entre outros fen

omenos. A pesquisa

e de import

ancia na recuperac¸

ao de

petr

oleo em

aguas profundas oce

anicas onde a concentrac¸

ao de sal e as diferen-

c¸as de temperatura geram estratiﬁcac¸

ao de tal forma que as ondas internas podem

afetar as operac¸

oes oﬀshore e as estruturas submersas.

Chapter 1

Introduction

Modelling waves is of great interest in the study of ocean dynamics. Internal

ocean waves, for example, appear when salt concentration and diﬀerences in tem-

perature generate stratiﬁcation. They can interact with the bottom topography and

submerged structures as well as with surface waves. In particular, in oil recovery

in deep ocean waters, internal waves can aﬀect oﬀshore operations and submerged

structures. Accurate reduced models are a ﬁrst step in producing eﬃcient compu-

tational methods for engineering problems in oceanography. This was the goal in

[24, 1].

To describe this nonlinear wave phenomenon in deep waters there are several

bidirectional models containing the Intermediate Long Wave (ILW) equation and

the Benjamin-Ono (BO) equation, starting from works such as [3, 9, 25, 14, 17]

to more recent papers such as [20, 6, 7, 8, 13]. In these models two fundamen-

tal mechanisms, nonlinearity and dispersion, are responsible for the main features

of the propagating wave. One of the most interesting behaviours observed is the

existence of solitary wave solutions with permanent shape. They are observed

when the steepening of a given wave front due to the nonlinearity and the attenua-

tion and ﬂattenning promoted by the dispersion are balanced on a particular scale.

Usually the contribution of nonlinearity is quantiﬁed by the non-dimensional non-

linearity parameter α, which is the ratio between the wave amplitude and the ﬂuid

layer thickness. It appears as a small non-zero parameter in the so-called weakly

nonlinear regime, and accompanies the nonlinear terms. On the other hand, the

dispersion parameter β is the squared ratio between the ﬂuid layer thickness and

the typical wavelength. It appears in the dispersion relation, making the phase

velocity a function of the wavenumber k. The balance that creates a solitary wave

is commonly obtained through a scaling relation between α and β, in the form of

a power law, for asymptotic values α ≪ 1 and β ≪ 1. In the water conﬁguration

considered here, it is the scaling α = O(

√

β) that leads to the ILW [14, 17]. In

the limit when one layer thickness tends to inﬁnity, the ILW equation becomes the

BO equation [3, 9, 25].

For all these models, the dependence on the vertical coordinate has been elim-

inated by focusing on speciﬁc regimes and using systematic asymptotic expansion

methods in small parameters. This results in a considerable simpliﬁcation of the

original Euler equations that leads to more eﬃcient computational methods than

the integration of the Euler system in the presence of a free interface. However,

the approximation needs to be accurate even for large values of the parameters

α and β. In other words, the model needs to be robust enough to cover several

regimes in which the viscosity eﬀects are negligible, justifying the use of the Eu-

ler equations. In [8], the authors compared weakly nonlinear models with experi-

mental data obtained by Koop and Butler in [16]. They found a divergence. This

motivated them to propose a strongly nonlinear model for ﬂat bottom that shares

the simplicity of the weakly nonlinear ones and extends its domain of validity.

The numerical results agree very well with the experimental data. This model is

generalized in the present work to consider an arbitrary sea bottom. We also im-

prove the asymptotic expansion to the next order of approximation in the pressure

term by taking a nonhydrostatic correction term. The resulting strongly nonlinear

model of higher order is more complicated than the previous one mentioned here,

but it has a weakly nonlinear version very similar to the strongly nonlinear model

of lower order. This fact implies that the weakly nonlinear higher-order model

should serve as a good model for moderate amplitude internal waves in a deep

water conﬁguration. We remark that the new models support bidirectional wave

propagation, so they are able to capture the reﬂected wave from the propagation

over a nonuniform sea bottom.

The models found in the literature consider ﬂat or slowly varying bottom to-

pography. Here, the model of Choi and Camassa is generalized to the case of

an arbitrary bottom topography by using the conformal mapping technique de-

scribed in [24]. We obtained a strongly nonlinear long-wave model like Choi and

Camassa’s, which is able to describe large amplitude internal solitary waves. A

system of two layers constrained to a region limited by a horizontal rigid lid at the

top and an arbitrary bottom topography is considered, as described in Fig. 2.1. The

upper layer is shallow compared with the characteristic wavelength at the interface

of the two-ﬂuid system, while the lower region is deeper. The nonlinear evolution

equations describe the behaviour of the internal wave elevation and mean upper-

velocity for this water conﬁguration. These Boussinesq-type equations contain

the ILW equation and the BO equation in the unidirectional wave regime. We

intend to use this model to study the interaction of waves with the bottom proﬁle,

in particular that of solitary waves. This is part of our future goals. The dynamics

described include wave scattering, dispersion and attenuation among other phe-

nomena.

The work is organizedas follows. In Chapter 2 the physical setting is presented

and as the main result, a reduced strongly nonlinear one-dimensional model is pro-

posed. Section 2.1 is devoted to obtaining a set of upper layer averaged equations

that will be completed with information provided by the lower layer in order to

derive the reduced model. The continuity of pressure at the interface establishes

a connection between both layers, as shown in Section 2.2. Through this condi-

tion we add the topography information to the averaged upper layer system. The

case when the depth of the bottom layer approaches inﬁnity is also considered. In

Section 2.3 the dispersion relations for the linearized models are computed. An

ILW equation with variable coeﬃcient and the BO equation are obtained from

the reduced models as unidirectional wave propagation models in Section 2.4. In

Section 2.5 theoretical solitary wave solutions are presented for the ILW equation

and for the Regularized ILW equation. The purpose of Chapter 3 is to exhibit a

model that improves the order in the asymptotic approximation in the pressure

term of the reduced model obtained in Chapter 2. To that end, in Section 3.1

one more term of the asymptotic expansion of the mean horizontal derivative of

pressure is added to the upper layer averaged equations. Then, in Section 3.2, the

approximation of the pressure at the interface is improved and a reduced strongly

nonlinear one-dimensional model of higher-order is obtained. The dispersion re-

lations for the higher-order model and for the previous model are compared with

the full dispersion relation originating from the Euler equations in Section 3.3.

Chapter 4 is devoted to the numerical resolution of the reduced model obtained

in Chapter 2. A hierarchy of one-dimensional models can be derived from this

strongly nonlinear model as shown in Section 4.1 by considering the diﬀerent

regimes (linear, weakly nonlinear or strongly nonlinear) as well as the ﬂat or cor-

rugated bottom cases. Numerical schemes based on the method of lines for all

models are described in Section 4.2 together with the study of their stability prop-

erties. The results from the Matlab implementations are shown in Sections 4.3, 4.4

and 4.5, including periodic topography experiments and solitary wave solutions.

Technical justiﬁcations for the manipulations done in Section 2.4 are provided in

Appendix A. The relation between the Hilbert transform on the strip (involved

in the models considered) and the Dirichlet-to-Neumann operator is presented in

Appendix B. Due to the nonlocal deﬁnition of the Hilbert transform on the strip,

it must be redeﬁned on the periodic domain used for numerical implementations,

as done in Appendix C.

Chapter 2

Derivation of the reduced model

In this chapter we generalize the work by Choi and Camassa [8]. Their asymptotic

technique for reducing a pair of two-dimensional (2D) systems of nonlinear partial

diﬀerential equations (PDEs) to a single one-dimensional (1D) system of PDEs

at an interface, is generalized to include very general submarine structures and

topographies at the bottom of the lower ﬂuid layer.

We start with a two-ﬂuid conﬁguration. Deﬁne the density of each inviscid,

immiscible, irrotational ﬂuid as ρ

for the upper layer and ρ

for the lower layer.

For a stable stratiﬁcation, ρ

> ρ

. Similarly, (u

, w

) denotes the velocity com-

ponents and p

the pressure, where i = 1, 2. The upper layer is assumed to have

an undisturbed thickness h

, much smaller than the characteristic wavelength of

the perturbed interface L > 0, hence the upper layer will be in the shallow water

regime. At the lower layer the irregular bottom is described by z = h

(h(x/l) −1).

The function h needs not to be continuous neither univalued, see for example

Fig. 2.1 where a polygonal shaped topography is sketched. We can assume that

h has compact support so the roughness is conﬁned to a ﬁnite interval. More-

η(x, t)

Figure 2.1: Two-ﬂuid system conﬁguration.

over h

is the undisturbed thickness of the lower layer outside the irregular bottom

region and it is comparable with the characteristic wavelength L, that character-

izes an intermediate depth regime. In the slowly varying bottom case we deﬁne

ε = L/l ≪ 1; when a more rapidly varying bottom is of concern, the horizontal

length scale for bottom irregularities l is such that h

< l ≪ L. The coordinate

system is positioned at the undisturbed interface between layers. The displace-

ment of the interface is denoted by η(x, t) and we may assume that initially it has

compact support.

The corresponding Euler equations are

+ w

= 0,

+ u

+ w

= −

+ u

+ w

= −

− g,

for i = 1, 2. Subscripts x, z and t stand for partial derivatives with respect to

spatial coordinates and time. The continuity condition at the interface z = η(x, t)

demands that

+ u

= w

, p

= p

namely, a kinematic condition for the material curve and no pressure jumps al-

lowed.

At the top we impose a rigid lid condition,

(x, h

, t) = 0,

commonly used in ocean and atmospheric models, while at the irregular imper-

meable bottom

−

′





+ w

= 0.

Introducing the dimensionless dispersion parameter β =





, it follows from

the shallowness of the upper layer that







= O





≪ 1.

From the continuity equation for i = 1 we have,

= O





= O







Let U



be the characteristic shallow layer speed. According to these

scalings, physical variables involved in the upper layer equations are non-dimen-

sionalized as follows:

x = L˜x, z = h

˜z, t =

t, η = h

˜η,

= U

˜u

, w



βU

˜w

, p

= (ρ

) ˜p

In a weakly nonlinear theory, η is usually scaled by a small parameter. Note that

here we have an O(1) scaling. This will lead to a strongly nonlinear model.

2.1 Reducing the upper layer dynamics to the inter-

face

The dimensionless equations for the upper layer (the tilde has been removed) are:

+ w

= 0,

+ u

+ w

= −p



+ u

+ w



= −p

− 1. (2.1)

The boundary conditions are

+ u

= w

and p

= p

at z = η(x, t), (2.2)

(x, 1, t) = 0.

Focusing on the upper shallow region, consider the following deﬁnition: for any

function f(x, z, t), let its associated mean-layer quantity f be

f(x, t) =

1 − η



f(x, z, t) dz.

By averaging we will reduce the 2D Euler equations to a 1D system.

Let η

= 1 − η. From the horizontal momentum equation we have

+ η

= −η

. (2.3)

We need to express each of these mean-layer quantities in terms of

and η.

The diﬃculty at this stage is breaking up the mean of square, and other general

quadratic terms, into individually averaged terms. To begin with, note that

(

)



dz − η

(x, η, t),

= η

− η

where u

is evaluated at the interface (x, z, t) = (x, η(x, t), t). So,

(

)

+ η

. (2.4)

Similarly

= u

+ (η

)

, (2.5)

and

2η

= η





. (2.6)

Therefore at z = η(x, t),

(

) =

(

)

+ u





From the kinematic condition Eq. (2.2)

= u

−

and by substitution,

(

) =

(

)

+ u

−





. (2.7)

On the other hand, integration by parts and incompressibility give

= −w

−



dz = −w



dz.

From Eq. (2.6),

= −w





. (2.8)

Substituting Eqs. (2.7) and (2.8) in Eq. (2.3), the following mean-layer equation

is derived

(η

)





= −η

. (2.9)

The incompressibility condition gives w

= η

at z = η(x, t). This, together

with Eq. (2.5) shows that

= u

+ (η

)

Substitution of w

into Eq. (2.2) leads to η

+ u

= u

+ (η

)

and

−η

= (η

)

. (2.10)

As pointed in [8], the system of Eqs. (2.9)–(2.10) was already considered in

[31, 5] for surface waves. In reducing (averaging) the 2D Euler equations to this

1D system no approximations have been made up to this point. Nevertheless, the

quantities u

· u

and p

prevent the closure of the system of Eqs. (2.9)–(2.10).

These quantities will be expressed in terms of η and

up to a certain order in

the dispersion parameter β. Note that until now, we still have not used the vertical

momentum equation and the continuity of pressure boundary condition. We start

by approximating

and then proceed to do the same for u

· u

The vertical momentum equation over a shallow layer suggests the following

asymptotic expansion in powers of β

f(x, z, t) = f

(0)

+ β f

(1)

+ O(β

)

for any of the functions u

, w

, p

. In fact, from Eq. (2.1), p

= −1 + O(β).

Integrating from η to z we have that

(x, z, t) − p

(x, η, t) = −(z −η) + O(β),

and the pressure continuity across the interface gives

(x, z, t) = p

(x, η, t) −(z − η) + O(β).

The pressure p

(x, η, t) should be non-dimensionalized in the same fashion as p

that is,

= ρ

˜p

Deﬁne P(x, t) = p



x, η(x, t), t



. Then

= P(x, t) − (z − η) + O(β),

which immediately yields

= P

(x, t) + η

+ O(β).

By averaging we get



(x, t) dz + η

+ O(β),

= P

(x, t) + η

+ O(β),





x, η(x, t), t





+ η

+ O(β). (2.11)

An approximation for P

will be obtained later from the Euler equations for

the lower ﬂuid layer. We now approximate the mean squared horizontal velocity

in terms of

and η.

In order to express

· u

as a function of u

and η, it should be pointed out

that the irrotational condition in non-dimensional variables is



i. e. u

= βw

. Hence



(0)



= 0 and as expected for shallow water ﬂows u

(0)

independent from z:

(0)

= u

(0)

(x, t). (2.12)

We now correct this ﬁrst order approximation. By using

= u

(0)

+ βu

(1)

+ O(β

) (2.13)

and Eq. (2.12) it is straightforward that

= u

(0)

+ 2βu

(1)

(0)

+ O(β



dz =



(0)

dz + 2β



(1)

(0)

dz + O(β

and

· u

= u

(0)

(1 − η) + 2η

βu

(0)

(1)

+ O(β

so that

· u

= u

(0)

· u

(0)

+ 2βu

(0)

(1)

+ O(β

). (2.14)

Also from Eq. (2.13),



dz = u

(0)

+ β

(1)

+ O(β

= u

(0)

+ β

(1)

+ O(β

· u

= u

(0)

· u

(0)

+ 2βu

(0)

(1)

+ O(β

). (2.15)

Thus Eqs. (2.14) and (2.15) lead to our desired approximation, namely that

· u

= η

· u

+ O(β

). (2.16)

Using Eq. (2.16), the nonlinear interfacial system (2.9) becomes

+ η



· u

+ O(β

)



= −η

+ η

(η

)

+ η

· u

= −η

+ O(β

From Eq. (2.10) one obtains that

+ η

· u

= −η

+ O(β

). (2.17)

After substitution of Eq. (2.11), the following set of approximate equations for the

upper layer was derived from Eqs. (2.9), (2.10):

+ (η

)

= 0,

· u

= −η

−





x, η(x, t), t





+ O(β),

or equivalently,











−η



(1 − η)



= 0,

· u

= −η

−





x, η(x, t), t





+ O(β).

(2.18)

We have almost closed our system of PDEs. Now we need to get an expression

for p

in order to close the system and also to establish a connection with the lower

ﬂuid layer.

2.2 Connecting the upper and lower layers

The coupling of the upper and lower layers is done through the pressure term. To

get an approximation for P

(x, t) =





x, η(x, t), t





from the Euler equations for

the lower ﬂuid layer, notice that out of the shallow water approximation

= O(1),

so that the following scaling relation

= O





= O(1)

follows from the continuity equation. At the interface, from the kinematic condi-

tions and the relations above, we have that

= O







= O







at z = η. Following these scalings introduce the dimensionless variables for the

lower region (with a tilde)

x = L˜x, z = L˜z, t =

t, η = h

˜η,

= (ρ

)˜p

, u



βU

˜u

, w



βU

˜w

This naturally suggests that we introduce the velocity potencial φ =

√

βU

φ.

Note that the deﬁnition for ˜z is diﬀerent from the one for the upper region, since

it involves the characteristic wavelength L instead of the vertical scale h

. This is

consistent since both scales are of the same order.

In these dimensionless variables, the Bernoulli law for the interface reads



β φ



+ φ



+ η +

P = C(t),

where the tilde has been ignored. C(t) is an arbitrary function of time. Then, up

to order β, the pressure derivative P

= −







(x,



βη, t)



+ O(β), (2.19)

where φ satisﬁes the Neumann problem with a free upper surface boundary con-

dition, given as











+ φ

= 0, on −

h(Lx/l)

≤ z ≤



β η(x, t),

= η

√

β η

, at z =

√

β η(x, t),

−

′

(Lx/l)φ

+ φ

= 0, at z = −

h(Lx/l)

(2.20)

Furthermore,





βφ





βη, t





βφ





βη, t



+ βφ





βη, t





βφ





βη, t



+ O(β),



βφ

(x, 0, t) + O(β),

where a Taylor expansion about z = 0 was performed.

Therefore,

= −





β φ

(x, 0, t)



+ O(β). (2.21)

As in the ﬂat bottom case [8], from Eq. (2.21) it is clear that it is suﬃcient to

ﬁnd the horizontal velocity φ

at z = 0 in order to obtain P

at the interface. Due

to the presence of the small parameter

√

β in problem (2.20), φ

(x, 0, t) can be

approximated by the horizontal velocity at z = 0 that comes from the linearized

problem around z = 0,











+ φ

= 0, on −

h(Lx/l)

≤ z ≤ 0,

= η

, at z = 0,

−

′

(Lx/l)φ

+ φ

= 0, at z = −

h(Lx/l)

(2.22)

In this systematic reduction we use Taylor expansion to ensure that

(x, 0, t) = φ





β η, t



+ O







= η



β η

+ O







z = 0

z = −

h(Lx/l)

ζ = 0

ζ = −

Figure 2.2: Conformal mapping, (x, z) =



x(ξ, ζ), z(ξ, ζ)



and therefore,

(x, 0, t) = η

+ O







To ﬁnd the horizontal velocity φ

(x, 0, t) in problem (2.22), a conformal map-

ping between the ﬂat strip ζ ∈



−

, 0



and the lower layer at rest is performed.

See Fig. 2.2.

The problem in conformal coordinates is











ξξ

+ φ

ζζ

= 0, on −

≤ ζ ≤ 0,

(ξ, 0, t) = M(ξ)η



x(ξ, 0), t



, at ζ = 0,

= 0, at ζ = −

(2.23)

where the previous Neumann condition at the top is now modiﬁed by M(ξ) =

(ξ, 0) which is the nonzero element of the Jacobian of the conformal mapping at

the unperturbed interface. As shown in [24], its exact expression is:

M(ξ

) = 1 −

∞



−∞



Lx(ξ, −h

/L)/l



cosh



πL

(ξ − ξ

)



dξ.

Moreover, the Jacobian along the unperturbed interface is an analytic function.

Hence a highly complex boundary proﬁle has been converted into a smooth

variable coeﬃcient in the equations.

To obtain the Neumann condition at the unperturbed interface in problem

(2.23), consider

= φ

+ φ

evaluated at z = 0 (equivalently ζ = 0):

(ξ, 0, t) = φ

(x, 0, t)x

(ξ, 0) + φ

(x, 0, t)z

(ξ, 0).

The Cauchy-Riemann relations and the fact that z(ξ, 0) = 0 and z

(ξ, 0) = 0 imply

that x

(ξ, 0) = 0, which leads to the Neumann condition in problem (2.23).

Since a conformal mapping was used in the coordinate transformation and

(ξ, 0) = 0, it is guaranteed that z

(ξ, 0) = x

(ξ, 0) is diﬀerent from zero. From

(ξ, 0, t) = φ

(x, 0, t)x

(ξ, 0), the velocity φ

(x, 0, t) is recovered as

(x, 0, t) =

(ξ, 0, t)

M(ξ)

The bottom Neumann condition is trivial in these new coordinates.

Notice that the terrain-following velocity component φ

(ξ, 0, t) is a tangential

derivative on the boundary for problem (2.23). Hence it can be obtained as the

Hilbert transform on the strip (see [15]) applied to the Neumann data. Namely

(ξ, 0, t) = T



(

ξ, 0, t)



(ξ),

and substituting the Neumann data from problem (2.23),

(ξ, 0, t) = T



ξ)η



ξ, 0), t





(ξ),

where

T[f](ξ) =



ξ)coth



(

ξ − ξ)



ξ (2.24)

is the Hilbert transform on the strip of height h. In this case, h = h

/L. The

singular integral must be interpreted as a Cauchy principal value. The eﬀect of

the two-dimensional undisturbed layer below the interface is being collapsed onto

a one-dimensional singular integral without any approximation. The results above

are used in (2.21) by noting that φ

is obtained after taking the time derivative of

problem (2.22). Therefore,

= −





M(ξ)



ξ)η



ξ, 0), t





(ξ)



+ O(β). (2.25)

Now, φ

(x, 0, t) is a tangential derivative on the ﬂat upper boundary for prob-

lem (2.22), whose domain is a corrugated strip. Hence, it is also expressed as a

Hilbert transform acting on Neumann data. Since

(x, 0, t) =



ξ(x, 0)





ξ)η



ξ, 0), t



coth



πL



ξ − ξ(x, 0)





ξ,

a Hilbert-like transform on the corrugated strip has been identiﬁed as:

[f](x) =



ξ(x, 0)





ξ) f



ξ, 0)



coth



πL



ξ − ξ(x, 0)





ξ,

which is not a convolution operator, unlike Eq. (2.24).

Finally, substituting the expression for P

obtained in Eq. (2.25) in the upper

layer averaged equations (2.18) gives











−



(1 − η)



= 0,



1 −







ξ(x, 0)





ξ)η



ξ, 0), t



coth



πL



ξ − ξ(x, 0)





ξ + O(β).

In a compact notation this becomes











−



(1 − η)



= 0,



1 −





M(ξ)



M(·)η



x(·, 0), t



(ξ) + O(β),

where the dot indicates the variable on which the operator T is applied.

It remains to make a few manipulations with this set of equations: elimi-

nate the second order derivative in time and write all spatial derivatives in the

ξ-variable.

Note that the ﬁrst equation is exact. According to it η



(1 − η)u



only the ﬁrst time derivative of

needs to enter the right-hand side of the second

equation.

In conclusion, the reduced one-dimensional internal wave model is:











−



(1 − η)



= 0,



1 −





M(ξ)



M(·)



(1 − η)





x(·, 0), t





(2.26)

The transform in the forcing term is in curvilinear coordinates. For practical

purposes both sides must be in the same coordinate system, which is readily ad-

justed via the conformal mapping: every x-derivative is equal to a ξ-derivative

divided by the Jacobian M(ξ). Therefore, system (2.26) in the terrain-following

coordinates reads











−

M(ξ)



(1 − η)



= 0,

M(ξ)



1 −





M(ξ)





(1 − η)



ξt



(2.27)

This is a Boussinesq-type system with variable (time independent) coeﬃcients

depending on M(ξ) for the perturbation of the interface η and the mean-layer

horizontal upper velocity

. We will show that this is a dispersive model, where

the dispersion term comes in through the Hilbert transform. Since no smallness

assumption was made on the wave amplitude up to now, the model derived is

strongly nonlinear. It involves a Hilbert transform on the strip characterizing the

presence of harmonic functions (hence the potential ﬂow) below the interface.

System (2.27) is a reduction of the original Euler equations constituted by

a pair of 2D-systems of PDEs to a single 1D-system of PDEs at the interface.

Instead of the integration of the Euler equations in the presence of a free interface,

a single 1D-system of PDEs is to be solved. Eﬃcient computational methods can

be produced for this accurate reduced model which governs, to leading order, a

complex two-dimensional problem.

Remarks:

1. If the bottom is ﬂat, M(ξ) = 1 and the same system derived in [8] is recov-

ered, which is a nice consistency check.

2. When the lower depth tends to inﬁnity (h

→ ∞) the limit for this model

is the same one obtained in [8] because the bottom is not seen anymore

(M(ξ) → 1 and x(

ξ, 0) →

ξ). Therefore

(x, 0, t) →





(1 − η)





˜x, t



˜x − x

d˜x = H





(1 − η)





(x),

where H is the usual Hilbert transform deﬁned as

H[f](x) =



f(˜x)

˜x − x

d˜x.

In this (shallow upper layer) inﬁnite lower layer regime, system (2.26) be-

comes











−



(1 − η)



= 0,



1 −









(1 − η)





(2.28)

3. The Fourier Transform (FT) of a Hilbert transform is easily computed. We

now make a comment regarding the use of FTs in numerical schemes. The

operators T[f] and H[f] have Fourier transforms



T[f] = icoth







H[f] = i sgn(k)

where the operator symbol multiplies the transform of f, which is

f. There-

fore in Eqs. (2.27) and (2.28) a pseudospectral scheme would apply a DFT

to the terms inside the square brackets. FFTs are only aplicable directly

when M(ξ) = 1 and the waves are weakly nonlinear.

2.3 Dispersion relation for the linearized model

Consider the ﬂat bottom case in system (2.26), that is, M ≡ 1:











−



(1 − η)



= 0,



1 −





β T





(1 − η)





+ O(β).

Its linearization around the undisturbed state η = 0, u

= 0 gives:











−

= 0,



1 −





β T





By diﬀerentiating once in t, η can be eliminated from the second equation:



1 −





β T



xtt



Let

= Ae

i(kx−ωt)

and substituting above,

i(kx−ωt)



−ω

−



1 −





β kω

−iωt



−ie

ikx



Since T[e

ikx

] = icoth





ikx



− 1



1 +

√

βk coth





, (2.29)

which is the correct approximation for the full dispersion relation



− 1



coth





√

βk coth





when kh

is near zero. A nonvanishing value of the parameter β in the dispersion

relation makes the phase velocity a function of the wave number k. Observe that

→ 0 as k → ∞, so bounded phase velocities are obtained as k becomes large.

This is a good property for numerical schemes.

In the limit h

→ ∞ the operator T becomes H. Since H[e

ikx

] = isgn(k)e

ikx

the dispersion relation for system (2.28) is



− 1



1 +

√

β |k|

and

→ 0 as k → ∞.

2.4 Unidirectional wave regime

For weakly nonlinear unidirectional waves and slowly varying topography, our

model reduces to a single ILW equation with variable coeﬃcients.

Consider again system (2.26), except for the Jacobian of the conformal map-

ping which now is

M(ξ

) = 1 −

∞



−∞



εx(ξ, −h

/L)



cosh



πL

(ξ − ξ

)



dξ,

since we assume a slowly varying bottom topography described in non-dimen-

sional coordinates as z = −

h(εx), with ε ≪ 1. The restriction to a slowly

varying topography is consistent with the objective of the present Section, which

is to ﬁnd equations to model the unidirectional wave propagation. Hence there will

be no reﬂection nor any backward evolution opposite to the propagation direction.

To study the weakly nonlinear regime, we introduce the typical amplitude a for

the perturbation of the interface and introduce the nonlinearity parameter α =

order O



√



. As usual in a weakly nonlinear theory we set η = αη

∗

. With this, the

original dimensional perturbation η is non-dimensionalized as η = αh

∗

= aη

∗

We also state that u

= αc

∗

, t =

∗

where c



− 1



. Depending on the root

chosen, there will be a right- or left-travelling wave.

Then, dropping the asterisks, Eq. (2.26) becomes











−



(1 − αη)



= 0,

+ α

− η



M(ξ)



ξ)



(1 − αη)





+ O(β).

(2.30)

Note that

+ O(α); η

+ O



α,





. (2.31)

As in [8], we look for a solution, up to a ﬁrst order correction in α and

√

β, in

the form

η = A

+ αA



β A

M(ξ)



ξ)



. (2.32)

Substituting in the system of Eqs. (2.30) up to order α,

√

β, two equations for

are obtained:

0 = A

+ 2αA

−



β A

M(ξ)



ξ)



and

0 =

+ α

−



+ 2α A



β A

M(ξ)



ξ)





−



M(ξ)



ξ)



+ O



, β, α





For compatibility A

= ±1 and to choose a right-going wave we take A

= −1.

Therefore

= −u

+ O



α,





, (2.33)

so u

= −

+ O



α,

√



and u

= −

+ O



α,

√



and the two equations are

consistent if A

= −1, A

= −

and A

= −

2ρ

. As a result, the evolution equation

for

−

√

M(ξ)



ξ)



= O



β, α

, α





For the elevation of the interface η a similar equation can be obtained through

asymptotic relations which permit (to leading order) to exchange derivatives in η

by derivatives in

, as well as time derivatives for spatial derivatives. This is a

consequence of (2.32). To begin with, use that

= −η

+ O



α,





√

M(ξ)



ξ)



= −

√

M(ξ)



ξ)η



+ O



, β, α





. (2.34)

In virtue of Eqs. (2.33) and (2.32)



−η

+ O



α,





−η + O



α,





αηη

+ O



, α





, (2.35)

and for similar reasons

+ η

= −

(

)

−

) −

√

M(ξ)



ξ)(

)



= −

(

)

+ O



, α



β, β



. (2.36)

Substituting all these expressions in the evolution equation for

we obtain

the evolution equation for the elevation of the interface,

+ η

−

αηη

−

√

M(ξ)



ξ)η



= O



β, α

, α





Finally, in curvilinear coordinates we have

M(ξ)

−

M(ξ)

ηη

−

√

M(ξ)

T[η

ξt

] = 0. (2.37)

This is an ILW equation with variable coeﬃ cients accounting for the slowly

varying bottom topography. Instead of the usual Hilbert transform on the half-

space, a Hilbert transform on the strip appears. The dispersion relation for the ﬂat

bottom case (M(ξ) = 1) is

ω =

1 +

√

k coth





. (2.38)

The equation reduces to a regularized dispersive model in analogy with the Benjamin-

Bona-Mahony equation (BBM), [4].

Remarks:

1. The constant coeﬃcient version of Eq. (2.37) diﬀers from the ILW consid-

ered in [8], Eq. (4.33), page 23, in that the latter has a dispersion term with

spatial derivatives only, as in the KdV equation. Both constant coeﬃcient

equations are equivalent up to the order considered since we can substi-

tute η

ξt

by −η

ξξ

in the dispersion term up to order O



β, α

√



. There are

two advantages for our choice. First, for every change from a Cartesian

x-derivative to a curvilinear ξ-derivative we need the presence of the met-

ric term M(ξ). Hence for the Camassa-Choi model with the second order

x-derivative we would end up with a variable coeﬃcient within the nonlo-

cal operator. Second, the regularized dispersive operator (namely with an

xt-derivative) leads to the stable dispersion relation (2.38) regarding numer-

ical schemes. Short waves have bounded propagation speeds. This does not

happen with the ILW equation considered in [8], whose dispersion relation

ω = k −

√

coth





2. One step remains to be explained in the substitution of Eq. (2.32) into the

second equation of system (2.30), namely why it is valid (for slowly varying

topography) that



√

M(ξ)



ξ)





√

M(ξ)



ξ)



+ O(β). (2.39)

This approximation was not presented in [8] since the present work con-

tains for the ﬁrst time the conformal mapping technique used for the lower

layer. The approximation (2.39) is justiﬁed through the construction of an

auxiliary PDE problem. See Appendix A for details.

3. For system (2.28) a similar unidirectional reduction can be obtained leading

+ η

−

αηη

−

√

H[η

] = 0,

which is a regularized Benjamin-Ono (BO) equation over an inﬁnite bottom

layer.

Hence Eq. (2.37) is a generalization of the BO equation for intermediate

depth and the presence of a topography. This equation is valid only when

backscattering is negligible.

2.5 Solitary wave solutions

The ILW equation we referred to in Section 2.4 and derived in [8, 14, 17] is of the

form

+ η

+ c

ηη

+ c

T[η

ξξ

] = 0 (2.40)

where c

= −

α and c

√

. Eq. (2.40) admits a family in the parameter θ of

solitary wave solutions [14, 8]

η(x) =

acos

cos

θ + sinh

(x/λ)

, x = ξ − ct, (2.41)

where

a =

θ tanθ

, λ =

, c = 1 −

θ cot(2θ ),

with 0 ≤ θ ≤ π/2. Alternatively, we consider a Regularized Intermediate Long

Wave equation

+ η

+ c

ηη

− c

T[η

ξt

] = 0 (2.42)

to the same order of approximation. It also admits the solitary wave solution

Eq. (2.41) with

a =

θ tanθ

1 +

θ cot(2θ )

, λ =

, c =

1 +

θ cot(2θ )

In the numerical section we will present a few experiments with solitary waves

over a ﬂat bottom. In the near future we intend to study solitary waves interacting

with highly varying topography and submarine structures.

Chapter 3

A higher-order reduced model

Since we want to study wave interaction with highly variable topographies and

submarine structures, we need to be able to account for higher order (vertical)

coupling terms between the two layers. Namely we want to investigate if these

higher order terms do indeed play a role in the dynamics. Hence in this chapter

we improve the model from the previous chapter regarding the pressure term by

allowing nonhydrostatic terms to come into play.

3.1 Higher-order upper layer equations

The purpose of this chapter is to improve the order of approximation of sys-

tem (2.27) by using higher precision approximations for the pressure term with

respect to the dispersion parameter β. Instead of order β, we seek order O(β

3/2

We start with the mean-layer equations (2.10) and (2.17) obtained in Section 2.1.

For convenience they are repeated here:

+ (η

)

= 0, (3.1)

· u

= −

+ O(β

). (3.2)

To approximate p

with order β

3/2

we need to expand p

(x, z, t) with one more

term:

(x, z, t) = p

(0)

+ βp

(1)

+ O(β

so its vertical derivative is expanded as

(x, z, t) = p

(0)

+ βp

(1)

+ O(β

). (3.3)

Again, from the vertical momentum equation

= −1 − β(w

+ u

+ w

)

we have that

(0)

= −1.

This is the hydrostatic contribution to the pressure. We also have that

(1)

= −(w

(0)

+ u

(0)

+ w

(0)

Even though the upper layer is shallow, through p

we can compute the leading

order nonhydrostatic correction.

Let D

= ∂t + u

(0)

∂x be the leading order material derivative. We rewrite the

expression above as

(1)

= −D

(0)

− w

(0)

. (3.4)

Let us express the quantities w

(0)

and w

(0)

in terms of η and

. We begin by

expanding the incompressibility equation to obtain

(0)

= −u

(0)

(x, t). (3.5)

Integrating Eq. (3.5) from η to z ≤ 1 and taking into account the z-independence

expressed by Eq. (2.12) we have that

(0)

(x, z, t) = −u

(0)

(x, t)(z − η) + w

(0)



z=η(x,t)

Now, from the kinematic condition in Eq. (2.2), to leading order

(0)



z=η(x,t)

= η

+ u

(0)

(x, z, t) = −u

(0)

(x, t)(z −η) + η

+ u

(0)

that is,

(0)

(x, z, t) = −u

(0)

(x, t)(z − η) + D

η. (3.6)

Substituting Eqs. (3.6) and (3.5) in Eq. (3.4):

(1)

= (z − η)



(0)

) − u

(0)



− D

(η).

Since

= u

(0)

+ O(β), (3.7)

it is also valid that

(1)

= (z − η)



(

) −



− D

(η) + O(β).

Recall that η

= 1 − η and now deﬁne

(x, t) =

η).

It can be shown

that

(x, t) = D

(

) −

+ O(β).

Therefore,

(1)

= (z − η)G

(x, t) −η

(x, t) + O(β)

Write the conservation of mass Eq. (2.10) in the form

(∂t +

∂x)η

+ η

= 0,

which together with the approximation (3.7) leads to

+ η

= O(β),

which is the same as D

η = η

+ O(β). Apply D

to it,

η = D

(η

) + O(β).

Expanding the right-hand side above we have

(η

) + O(β) = η

(

) +

(η

) + O(β),

= η

(

+ u

− u

) + O(β),

= η



(

) −



+ O(β).

Then D

η = η



(

) −



+ O(β) as desired.

and substituting in the asymptotic expansion Eq. (3.3) for p

we have

(x, z, t) = −1 + β(z −1)G

(x, t) + O(β

Integrating from z = η(x, t) to z ≤ 1 we obtain

(x, z, t) = P(x, t) − (z −η) + βG

(x, t)



(z − 1)

−

(η − 1)



+ O(β

Diﬀerentiating once in x,

= η

+ P

(x, t) + β



(x, t)



(z − 1)

−

(η − 1)



+ O(β

)

and taking means

= η

+ P

(x, t) −



(x, t)



+ O(β

Substituting in Eq. (3.2) we have

· u

= −



+ P

(x, t) −



(x, t)



+ O(β

). (3.8)

If the lower ﬂuid layer is neglected and P is regarded as the external pressure

applied to the free surface, Eqs. (3.1) and (3.8) are the complete set of evolu-

tion equations for one homogeneous layer derived by Su and Gardner in [27] and

independently by Green and Naghdi in [11].

3.2 Improved approximation for pressure at the in-

terface

Now we want an approximation of order O(β

) for P

(x, t) =





x, η(x, t), t





from the Euler equations for the lower ﬂuid layer. This order of approximation

is suﬃcient to make the nonhydrostatic order β terms explicit in the asymptotic

expansion. Again, the scale

√

β comes from the lower layer reduction.

From the Bernoulli law for the lower layer:

P(x, t) = −





βφ

(φ

+ φ

) + η + C(t)





√

βη(x,t)

Using a Taylor expansion about z = 0 we obtain that

P(x, t) = −



η +





z=0



βη φ

z=0







z=0

+ φ



z=0



+ C(t)



+ O(β

Since from Eq. (2.20) we have that φ

= η

√

βη

= η

+ O



√



at z =

√

βη(x, t), it follows that

z=0

= φ

√

βη

+ O







= η

+ O







and

z=0

= φ

√

βη

+ O







= η

+ O







Therefore

P(x, t) = −



η +



β φ

z=0

+ βηη





z=0

+ η



+ C(t)



+ O(β

)

and it is easy to take x-derivatives since all quantities are evaluated at z = 0:

(x, t) = −





β φ

z=0

+ β



ηη



z=0



+ O(β

Note from the previous chapter that φ

z=0

is already known up to order

√

β,

namely that



x(ξ, 0), 0, t



M(ξ)



ξ)η



ξ, 0), t





(ξ) + O







(3.9)

which is all we need to approximate



z=0

We will obtain φ

z=0

with order of approximation β via the Hilbert transform

on the corrugated strip. Assuming that our potential problem for the lower deep

layer is deﬁned in a region surrounding the physical domain containing the in-

terface at rest z = 0, we restrict our potential problem to the unperturbed region.

There the potential problem satisﬁes a certain upper Neumann boundary condition

(φ

z=0

) to be determined up to order β. This order of approximation is suﬃcient

to obtain an x-derivative of order β because we are solving a linear problem and

we have the Hilbert transform connecting these two derivatives. The calculation

is as follows.

To obtain φ

z=0

a Taylor expansion is used as before:





βη



= φ

(x, 0) +



βηφ

(x, 0) + O(β).

Since the boundary condition φ



√

βη



is known we have

(x, 0) = η



βη

(x, 0) −



βηφ

(x, 0) + O(β)

and due to the Laplace equation,

(x, 0) = η



βη

(x, 0) +



βηφ

(x, 0) + O(β)

which is the same as

(x, 0) = η





ηφ

(x, 0)



+ O(β).

Using Eq. (3.9) in the expression above,

(x, 0) = η





M(ξ)

T[M(

ξ)η

]



+ O(β),

= η

√

M(ξ)



M(ξ)

T[M(

ξ)η

]



+ O(β).

Thus by means of the Hilbert transform on the corrugated strip,

(x, 0) =

M(ξ)



ξ)φ

(x, 0)



M(ξ)







ξ)







√

ξ)



ξ)

T[M(ξ

′

)η

]















+ O(β),

M(ξ)







ξ)η





ξ)

T[M(ξ

′

)η

]









+ O(β).

It is easy to take a time-derivative of this expression since the coeﬃcient M is

independent of t:

(x, 0) =

M(ξ)







ξ)η





ξ)

T[M(ξ

′

)η

]









+ O(β),

which together with



(1 − η)





(1 − η)



ξt

/M(ξ)

leads to



β φ

z=0

√

M(ξ)





(1 − η)



ξt



M(ξ)



ξ)





(1 − η)







ξt

+ O(β

and as we saw





z=0











M(ξ)



ξ)η











+ O(β

2M(ξ)









M(ξ)



ξ)η











+ O(β

2M(ξ)









M(ξ)





(1 − η)













+ O(β

Summarizing, the higher order pressure term connecting the top and lower

layer is

(x, t) = −







√

M(ξ)





(1 − η)



ξt



M(ξ)



M(ξ)





(1 − η)







ξt

2M(ξ)









M(ξ)





(1 − η)













M(ξ)



ηη









+ O(β

This order of approximation is compatible with the nonhydrostatic corrections

added to the upper shallow water layer model since it makes the order β terms

explicit.

Notice that composition of the Hilbert operator T arises in this case. A sim-

ilar situation also occurred in the fully dispersive Boussinesq model obtained by

Matsuno in [19] and Artiles and Nachbin in [1] for surface gravity waves.

The above pressure term is substituted into the system of Eqs. (3.1) and (3.8)

in curvilinear coordinates, that is,











M(ξ)



(1 − η)



M(ξ)

= −

M(ξ)

(1 − η)

3M(ξ)



(1 − η)



− P

+ O(β

)

where

(ξ, t) =

M(ξ)

ξt

M(ξ)



M(ξ)



−

M(ξ)

As a result, the strongly nonlinear model of order β

is:











M(ξ)



(1 − η)



M(ξ)



1 −



√

M(ξ)



(1 − η)



ξt

1 − η

3M(ξ)



(1 − η)



M(ξ)









(1 − η)



ξt

M(ξ)



(1 − η)









M(ξ)



M(ξ)



(1 − η)





ξt

2M(ξ)









M(ξ)





(1 − η)













+ O(β

For the weakly nonlinear model of order β

introduce η

∗

such that

η = αη

∗

= αu

∗

with α = O



√



, a typical scaling used for solitary waves. After dropping the

asterisks we have











M(ξ)



(1 − αη)



M(ξ)



1 −



√

M(ξ)



(1 − αη)



ξt

3M(ξ)



M(ξ)

ξt



+ O(β

The higher-order weakly nonlinear model has exactly the same form as the

lower-order strongly nonlinear model when the last term, namely,

3M(ξ)



M(ξ)

ξt



of the weakly nonlinear model is neglected. This implies that the weakly nonlinear

higher-order model should serve as a good model for moderate amplitude inter-

nal waves in a deep water conﬁguration. Furthermore, when the additional term

from the upper layer is included, the linear dispersion relation for the higher-order

weakly nonlinear model becomes the closest to the exact linear dispersion relation

when compared to lower-order models, as will be shown in the next Section. In

other words, the weakly nonlinear higher-order model might have a large domain

of validity so that the model can be used for a moderate (although still large) ra-

tio of h

, where the eﬀects of bottom topography are more pronounced. This

might be a justiﬁcation for using a higher-order weakly nonlinear model and will

be thoroughly explored in the near future.

3.3 Dispersion relation for the higher-order model.

Comparison with the previous model

To obtain the dispersion relation for the improved model, consider its linearization

around the undisturbed state so that













1 −





β T[u

]

ξt

ξξt

in the presence of a ﬂat bottom.

Taking derivatives once in t, η can be eliminated from the second equation,



1 −



ξξ

−



β T[u

]

ξtt

−

ξξtt

= 0.

Let

= Ae

i(kx−ωt)

. Substituting above and using again that T[e

ikx

] = icoth





ikx

we get



1 +



β k coth







− 1



that is,



− 1



1 +

√

β k coth





. (3.10)

Again,

→ 0 as k → ∞, so bounded phase velocities are obtained as k becomes

large.

Now, let us make a comparison between the diﬀerent dispersion relations ob-

tained throughout this work. Initially, we have the dimensional full dispersion

relation

g(ρ

− ρ

kcoth(kh

) + ρ

kcoth(kh

)

(3.11)

that comes from the linearized Euler equations around the undisturbed state, see

for example [18]. To compare it with the dimensionless dispersion relations (2.29)

and (3.10) it must be taken into account that because of the non-dimensionalization,

k =

, ω =

ω.

Therefore, (2.29) in dimensional form becomes

g(ρ

− ρ

+ ρ

kcoth(kh

)

. (3.12)

As it has been stated, the reduced model (2.27) captures the dispersion relation for

the shallow water (long waves) regime in the upper layer since

kcoth(kh

) →

as kh

→ 0.

On the other side, the relation (3.10) in dimensional form is

g(ρ

− ρ

+ ρ

kcoth(kh

)

Notice that both approaches are fully dispersive regarding the bottom layer, since

the second coth is completely retained. For the shallow water upper layer regime

(kh

near zero) we obtain again the correct approximation of the full dispersion

relation, but here the term ρ

kcoth(kh

) in the denominator is expanded with one

more term, namely

coth(kh

) =



1 +

(kh

)

+ O



(kh

)





and consequently

= ω

+ O



(kh

)



while

= ω

+ O



(kh

)



This means that the linear dispersion relation from the higher-order nonlinear

model ω

is closer to the full (exact) linear dispersion relation ω

than the lin-

ear dispersion relation from the lower-order model ω

is. See Fig. 3.1 and a detail

in Fig. 3.2 where the corresponding phase velocities are depicted.

The inclusion of the higher order pressure term has improved the accuracy of

the phase speed over a much wider wavenumber band. This is very important in

reﬂection-transmission problems as shown by Mu

noz and Nachbin in [23].

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

1.2

1.4

1.6

1.8

2.2

2.4

2.6

ω/k

Figure 3.1: Phase velocities for ρ

= 1, ρ

= 2, h

= 1, h

= 2, β = 0.01. Dotted

line: full phase velocity, dashed line: phase velocity for the higher order model,

solid line: phase velocity for the lower order model.

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

2.05

2.1

2.15

2.2

2.25

ω/k

Figure 3.2: Detail from Fig. 3.1.

Chapter 4

Numerical results

4.1 Hierarchy of one-dimensional models

For numerical implementations, we ﬁrst normalize the shallowwater velocity c



− 1



of system (2.27) by setting η = η

∗

= c

∗

, t =

∗

. Dropping the

asterisks, the following Strongly Nonlinear Corrugated Bottom Model (SNCM) is

obtained,











−

M(ξ)



(1 − η)



= 0,

M(ξ)

−

M(ξ)



M(ξ)

[0,2ℓ]



(1 − η)



ξt

(4.1)

We will work on the periodic domain ξ ∈ Π[0, 2ℓ], so that instead of the

operator T its periodic version T

[0,2ℓ]

appears above. See Appendix C for the

deﬁnition of T

[0,2ℓ]

. The choice of a computational periodic domain was made to

be able to use spectral methods. To avoid the inﬂuence of the boundaries on the

evolution of the perturbation η interacting with the bottom proﬁle we added two

homogeneuous (ﬂat) regions at the extremes and keep the dynamics away from

these regions. Hence both proﬁle and initial disturbance will be deﬁned and kept

away from ξ = 0 = 2ℓ.

From system (4.1), a hierarchy of one-dimensional models can be derived by

considering the diﬀerent regimes (linear, weakly nonlinear or strongly nonlinear)

as well as the ﬂat or corrugated bottom cases. The Weakly Nonlinear Corrugated

Bottom Model (WNCM) was already obtained in (2.30). In curvilinear coordi-

nates for a periodic domain it reads











−

M(ξ)



(1 − αη)



= 0,

M(ξ)

−

M(ξ)



M(ξ)

[0,2ℓ]





ξt

(4.2)

Setting α = 0 we obtained the Linear Corrugated Bottom Model (LCM)











−

M(ξ)

= 0,

−

M(ξ)



M(ξ)

[0,2ℓ]





ξt

(4.3)

The ﬂat bottom versions are obtained by simply taking M(ξ) = 1 for all ξ ∈

Π[0, 2ℓ]. To ﬁx a notation, let us use the abbreviations in Table 4.1 to refer to each

model.

4.2 Method of lines

To ﬁnd the solution for the initial value problem of systems SNCM, WNCM,

SNFM, WNFM, LCM is a nontrivial task. That is why we resort to numerical

methods to ﬁnd approximate solutions.

Linear Weakly non-

linear

Strongly non-

linear

Flat bottom LFM WNFM SNFM

Rough

bottom

LCM WNCM SNCM

Table 4.1: Abbreviations for the six diﬀerent models.

In order to use the method of lines to solve numerically the systems of equa-

tions (4.1, 4.2, 4.3) and their ﬂat bottom versions, let us rewrite them in a more

convenient way,











= E(η,

= F(η,

(4.4)

where V is an auxiliary variable deﬁned for each model in Table 4.2. The corre-

sponding vector ﬁeld (E, F) is deﬁned in Table 4.3. The mean-layer horizontal

upper velocity

can be recovered from η and V by inverting the relations in Ta-

ble 4.2 in a way to be speciﬁed later on for each case. For the time being, let us

assume that

= Ψ(η, V), given a certain operator Ψ.

According to the method of lines, we can discretize in space and solve a cou-

pled system of ODEs by a ﬁnite diﬀerence formula in t like, for example, a Runge-

Kutta integration scheme or a predictor-corrector solver with an Adams-Bashforth

predictor and an Adams-Moulton corrector.

First, an approximation scheme for the ξ-derivatives involved in the right-hand

side of the systems above must be used for the discretization in space. A choice

V Flat bottom Rough bottom

Linear

−

√

[0,2ℓ]





−

√

M(ξ )

[0,2ℓ]





Weakly

Nonlinear

−

√

[0,2ℓ]





−

√

M(ξ )

[0,2ℓ]





Strongly

Nonlinear

−

√

[0,2ℓ]



(1 − η)



−

√

M(ξ )

[0,2ℓ]



(1 − η)



Table 4.2: The auxiliary variable V.

(E, F) Flat bottom Rough bottom

Linear



, η



M(ξ )



, η



Weakly

Nonlinear





(1 − αη)



, η

− α



M(ξ )





(1 − αη)



, η

− α



Strongly

Nonlinear





(1 − η)



, η

−



M(ξ )





(1 − η)



, η

−



Table 4.3: The ﬁeld (E, F).

is to approximate the ξ-derivative of a function f(ξ) by the fourth order, ﬁve point

formula

(ξ

) =

8(f

j+1

− f

j−1

) + f

j−2

− f

j+2

12∆ξ

+ O(∆ξ

), (4.5)

where f

= f(ξ

), ξ

= j∆ξ, ∆ξ = 2ℓ/N, j = 1, . . . , N. So the spatial discretization

for the LFM with β = 0 (which is just the wave equation) leads to the system of

Ordinary Diﬀerential Equations (ODEs)











= C

= Cη,

where η = [η

, . . . , η

]

, with η

= η(j∆ξ, t) and

= [u

, . . . ,

]

, with

(j∆ξ, t), j = 1, . . . , N. C is the skew-symmetric, circulant matrix,

C =

∆ξ







0 2/3 −1/12 0 ··· 0 1/12 −2/3

−2/3 0 2/3 −1/12

1/12

1/12 −2/3 0 2/3

0 1/12 −2/3 0

0 2/3 −1/12

−1/12

−2/3 0 2/3

2/3 −1/12 0 1/12 −2/3 0







Since it is a skew-symmetric matrix, it has imaginary eigenvalues. The same is

true for the block matrix













which has the same eigenvalues of C with double multiplicity.

Another choice to approximate the ξ-derivatives is to use the spectral ξ-derivative,

whose corresponding matrix is

D =







0 −

cot(

1∆ξ

)

−

cot(

1∆ξ

)

cot(

2∆ξ

)

cot(

2∆ξ

)

−

cot(

3∆ξ

)

−

cot(

3∆ξ

)

−

cot(

(N−1)∆ξ

)

−

cot(

(N−1)∆ξ

) 0







Here we also have a skew-symmetric, Toeplitz, circulant matrix with imaginary

eigenvalues ik, k = −N/2 + 1, . . . , N/2 − 1, with zero having multiplicity 2. With

it, the discretization of the non-dispersive LFM leads to an ODEs system for η and

involving the block matrix













This block matrix has the same imaginary eigenvalues of D with double multi-

plicity.

The rule of thumb for stability (valid for normal matrices) is [28]: the method

of lines is stable if the eigenvalues of the linearized spatial discretization operator,

scaled by ∆t, lie in the stability region of the time-discretization operator.

In Fig. 4.1 we depict the stability regions for the fourth order Runge-Kutta in-

tegration scheme (RK4), the ﬁfth order, four step Adams-Moulton scheme (AM4)

and the fourth order, three step Adams-Moulton scheme (AM3). Along the imag-

inary axis (where the eigenvalues of the linearized spatial discretization operators

lie), RK4 is less restrictive, allowing larger time steps and numerical evolution

over longer time intervals, as we will see in the experiments below. For a better

visualization, in Fig. 4.2 we compute the ampliﬁcation factor

R(z)

where z = λ∆t

and λ represents the largest (in the sense of absolute value) eigenvalue of the lin-

earized spatial discretization operator, see for example [2]. The AM3 scheme

has the smallest stability interval along the imaginary axis since the ampliﬁca-

tion factor becomes greater than one very quickly. This implies a more restrictive

stability condition when using the AM3 scheme to solve the bidirectional wave

equation, for example. However, for the regularized dispersive systems consid-

ered here, we obtained less instability since the phase velocity actually decreases

as the wavenumber grows accommodating high wave-numbers better than in the

hyperbolic case (see the dispersion relation in Section 2.3). Dispersion comes in

as a physical regularization in comparison with its underlying hyperbolic counter-

part. Still, the classical fourth order Runge-Kutta seems to be the best choice.

The spectral approximation of the ξ-derivative is more accurate than the ﬁve

point formula exhibited above. However, we cannot use numerical velocity one

(∆t = ∆ξ) with it, in correspondence with the theoretical velocity, because of

stability restrictions of the method of lines already discussed. The reason for

choosing the Courant number as one is to avoid a numerical delay of the travelling

wave speed that could interfere with the expected delay that may result from the

interaction of the wave with a rapidly-varying bottom proﬁle. This particular issue

will be investigated in the near future. For the time being, we choose the ﬁve point

formula (4.5) in all the numerical experiments presented in this chapter.

The time domain is discretized as t

= n∆t, n = 0, . . . , T

ﬁnal

with ∆t = ∆ξ.

−3 −2 −1 0 1 2 3

−3

−2

−1

RK4

AM3

AM4

Re(z)

Im(z)

Figure 4.1: Stability regions for RK4, AM4, AM3.

0 0.5 1 1.5 2 2.5 3

0.4

0.6

0.8

1.2

1.4

1.6

1.8

2.2

2.4

RK4

AM3

AM4

R(z)

Figure 4.2: Ampliﬁcation factor along the imaginary axis for RK4, AM3, AM4.

So the numerical velocity equals the theoretical one in the LFM above and the

Courant number σ = 1 ×∆t/∆ξ equals one. This choice is maintained throughout

the numerical experiments with the diﬀerent systems.

The time evolution step for η and V using RK4 is

n+1

= η

∆t



K1 + 2K2 + 2K3 + K4



n+1

= V

∆t



KK1 + 2KK2 + 2KK3 + KK4



where η

and V

are N × 1 vectors with components η

= η(j∆ξ, n∆t) and V

V(j∆ξ, n∆t) respectively. After each evolution step the N × 1 vector

with

components

= u

(j∆ξ, n∆t) must be recovered from η

and V

via u

Ψ(η

, V

Also, recall that

K1 = E(η

KK1 = F(η

= Ψ(η

+ 0.5∆tK1, V

+ 0.5∆tKK1),

K2 = E(η

+ 0.5∆tK1,

KK2 = F(η

+ 0.5∆tK1,

= Ψ(η

+ 0.5∆tK2, V

+ 0.5∆tKK2),

K3 = E(η

+ 0.5∆tK2,

KK3 = F(η

+ 0.5∆tK2,

= Ψ(η

+ ∆tK3, V

+ ∆tKK3),

K4 = E(η

+ ∆tK3,

KK4 = F(η

+ ∆tK3,

Notice that the projection from (η, V) back to

must be done after each stage.

These intermediate values are denoted by

, i = 1, 2, 3.

Now we will specify how to deal with the operator Ψ in order to recover

terms of V and η. Since the Fast Fourier Transform (FFT) of a Hilbert transform

is easily computed, as well as the FFT for a ξ-derivative, we can go to frequency

space and solve

in terms of V for the linear and constant coeﬃcients relation

that appears in Table 4.2 for both LFM and WNFM. In Fourier space we have



V(k) =



1 +



kπ

ℓ

coth



kπh

ℓL





(k)

for k = −N/2 + 1, . . . , N/2. Therefore,



(k) =



V(k)



1 +

√

kπ

ℓ

coth



kπh

ℓL



since the denominator is always diﬀerent from zero.

The remaining models are essentially diﬀerent from the previous ones, as well

as from the terrain following Boussinesq system obtained in [24] and solved nu-

merically in [22] and the extended Boussinesq equations derived by Nwogu and

solved numerically in [29]. This is due to the fact that in the present work, either

the dispersive term has a nonlinear dependence between η and

, or it has a vari-

able coeﬃcient accompanying it. The use of the FFT is no longer straightforward,

so we use a matrix formulation to ﬁnd

in terms of V

and η.

For the SNCM (4.1) we have

V =

−



M(ξ)

[0,2ℓ]



(1 − η)



. (4.6)

Note a nonlinear term with the T operator and also a variable coeﬃcient in

front of it. The discrete version in matrix form of Eq. (4.6) reads

V =



1 −



−1

TDA



, (4.7)

where the N × N matrices involved are

• M

= M(ξ

), and zero elsewhere,

• A

= 1 − η(ξ

), and zero elsewhere,

• D is the Fourier spectral spatial diﬀerentiation matrix,

• T is the convolution matrix for T

[0,2ℓ]

. The composition of convolution and

diﬀerentiation can be computed with the help of the Fourier transform ma-

trix F,

= w

(j−1)(l−1)

, w = e

2πi/N

leading to the operator with symbol matrix











−

jπ

ℓ

coth



jπh

ℓL



, i = j = 1, . . . , N/2,

−

(j−N)π

ℓ

coth



(j−N)πh

ℓL



, i = j = N/2 + 1, . . . , N − 1,

0 elsewhere,

where

TD = DT =

FΛ

Although the original expression (4.6) is nonlinear, the relation (4.7) is a linear

algebraic system to be solved for

since η and V are already known at the current

time step n + 1. So at this stage, by using a spectral matrix instead of an FFT, we

are only paying a price in complexity but not in accuracy. Table 4.4 sumarizes the

discretizations used for each model, including the discretizations for the LFM and

WNFM cases. These cases were implemented with the help of the FFT and also

using a matrix formulation as a way to validate the method for the other models.

We also implemented in Matlab two predictor-corrector schemes, one that uses

a third order, three step, explicit Adams-Bashforth solver (AB3) as predictor and

a fourth order, three step implicit Adams-Moulton (AM3) as corrector. The other

scheme uses a fourth order, four step Adams-Bashforth solver (AB4) as predictor

and a ﬁfth order, four step implicit Adams-Moulton (AM4) as corrector.

V Flat bottom Rough bottom

Linear



1 −

√





1 −

√

−1



Weakly

Nonlinear



1 −

√





1 −

√

−1



Strongly

Nonlinear



1 −

√

TDA





1 −

√

−1

TDA



Table 4.4: Discretizations for the relations in Table 4.2.

The AB3 scheme used is:

n+1

= η

∆t



23E

− 16E

n−1

+ 5E

n−2



n+1

∆t



23F

− 16F

n−1

+ 5F

n−2



where E

= E(j∆ξ, n∆t) and F

= F(j∆ξ, n∆t). The AM3 scheme used is:

n+1

= η

∆t



n+1

+ 19E

− 5E

n−1

+ E

n−2



n+1

∆t



n+1

+ 19F

− 5F

n−1

+ F

n−2



The AB4 scheme used is:

n+1

= η

∆t



55E

− 59E

n−1

+ 37E

n−2

− 9E

n−3



n+1

= u

∆t



55F

− 59F

n−1

+ 37F

n−2

− 9F

n−3



The AM4 scheme used is:

n+1

= η

∆t

720



251E

n+1

+ 646E

− 264E

n−1

+ 106E

n−2

− 19E

n−3



n+1

∆t

720



251F

n+1

+ 646F

− 264F

n−1

+ 106F

n−2

− 19F

n−3



To initialize both predictor-corrector schemes, the RK4 described above is

employed.

4.3 Flat bottom experiments

Example 4.1. We ﬁrst consider the LFM for the non dispersive case β = 0,











= η

(4.8)

This is just the bidirectional wave equation for η

− η

ξξ

= 0

with initial contidions











η(ξ, 0) = η

(ξ),

(ξ, 0) =

(ξ).

This model (whose exact solution is well known) is useful for validating and com-

paring the numerical schemes, since it is more demanding than the dispersive

models regarding numerical stability, as we commented in the previous Section.

Consider the initial condition

(ξ) = 0.5e

−a(ξ−2π)

/64

, ξ ∈ Π[0, 16π],

with a = 50, η

(ξ, 0) =

(ξ) = −η

(ξ), where Π[0, 16π] is the interval [0, 16π]

with periodic boundary conditions. If we set

initially to −η

(ξ) then the Rie-

mann invariants for the nondispersive LFM, A = η +

and B = η − u

, will be

A(ξ, t) = 0 and B(ξ, t) = η(ξ − t, 0). Therefore

η =

A + B

= η(ξ − t, 0)

will be a right-travelling wave as shown in Fig. 4.3, where the numerical solution

obtained with the AM4 scheme is shown until time t = 39.2699. The solution for

the ﬁnal time is also plotted in Fig. 4.4, the absolute error is 0.00098872, a little

less than with RK4 (0.0013). Nevertheless, with RK4 we are able to advancemuch

more in time, until t = 151.1891, while with AM4 instabilities set for t = 49.0874.

AM3 is unstable as early as t = 9.8175.

Example 4.2. Another interesting example from the wave equation is that of the

ﬁssion of the wave. Take the initial condition for η as

(ξ) = 0.5e

−50(ξ−8π)

/64

, ξ ∈ Π[0, 16π],

and

= 0. If we set

initially to zero then the Riemann invariants will be

A(ξ, t) = η(ξ + t, 0) and B(ξ, t) = η(ξ − t, 0). Therefore

η =

A + B

η(ξ + t, 0) + η(ξ − t, 0)

−0.2

0.2

0.4

0.6

Figure 4.3: Travelling wave on a periodic domain Π[0, 16π]. The numerical

solution was obtained using AM4 with N = 512, ∆ξ = 2π/N = 0.098175,

∆t = ∆ξ = 0.098175.

0 10 20 30 40 50 60

−0.1

0.1

0.2

0.3

0.4

0.5

0.6

Figure 4.4: Travelling wave. Dotted line: numerical solution for the nondispersive

LFM using AM4 for t = 39.2699 and N = 512, ∆ξ = 2l/N = 0.098175, ∆t =

∆ξ = 0.098175, dashed line: initial condition, solid line: exact solution.

−0.2

0.2

0.4

0.6

Figure 4.5: Fission of the wave on a periodic domain Π[0, 16π].

and the amplitude is half of the original one, see Fig. (4.5). This is consistent

with D’Alembert’s solution for the wave equation. When the two travelling waves

coincide (overlap) in space and time the initial condition is recovered with an

error of 0.00374 by the RK4 method with N = 512, ∆ξ = 2l/N = 0.098175,

∆t = ∆ξ = 0.098175, see Fig. 4.6. This is a nice consistency check for the

numerical conservation of mass of two waves colliding.

Example 4.3. Let us study the numerical solutions for the LFM for the dispersive

case β  0. The initial value problem for this case can be solved explicitly by

means of Fourier Series as follows.

Consider the initial value problem (IVP) for the LFM on the periodic domain

0 10 20 30 40 50 60

−0.1

0.1

0.2

0.3

0.4

0.5

0.6

Figure 4.6: Fission of the wave. Dotted line: numerical solution for the nondis-

persive LFM using RK4 for t = 50.2655 and N = 512, ∆ξ = 2l/N = 0.098175,

∆t = ∆ξ = 0.098175, dashed line: initial condition, solid line: exact solution.

ξ ∈ Π[0, 2ℓ ],











−

= 0,

− η



[0,2ℓ]





ξt

η(ξ, 0) = η

(ξ),

(ξ, 0) = u

(ξ).

(4.9)

First, let us perform the change of variables

ξ = πξ/ℓ, t = πt/ℓ to the standard

periodic domain ξ ∈ Π[0, 2π]. The time variable is modiﬁed in order to maintain

the wave velocity as one for the hyperbolic regime (β = 0). In these coordinates

the IVP reads











−

= 0,

− η



[0,2ℓ]





η(

ξ, 0) = η

(ξ),

(ξ, 0) = u

(

ξ).

(4.10)

The symbol of the operator T

[0,2ℓ]

[

]

(that is the composition of one spatial deriva-

tive with the Hilbert transform) in the new coordinates is

−

kπ

ℓ

coth



kπ

ℓ



see Appendix C. Applying Fourier Transform (see Appendix C, Eq. (C.2) for its

deﬁnition) to problem (4.10) we have











ˆη

= ik





1 +



kπ

ℓ

coth



kπ

ℓ



= ikˆη, for k  0

ˆη(k, 0) =



(k),

(k, 0) =



(k).

(4.11)

Substitute



from the ﬁrst equation into the second:

ˆη

t t

= −

1 +

√

kπ

ℓ

coth



kπ

ℓ



ˆη = −ω

(k)ˆη.

The general solution for this ODE is

ˆη(k, t) = c

exp



iω(k)t



+ c

exp



−iω(k)t



, k  0,

where c

and c

are two functions of the wavenumber k. Using Fourier Series we

can write the general solution for η, at least in a formal manner, as

η(

ξ, t) =

2π

∞



k=−∞

(k) exp



iω(k)t



2π

∞



k=−∞

(k) exp



−iω(k)t



Each term represents one wave mode, see [30]. Since the dispersion relation ω(k)

is odd, each wave propagates in one direction: the ﬁrst wave travels to the left, the

second one to the right. Returning to Fourier space, from the initial condition in

(4.11) we have,











(k) + c

(k) =



(k),

(k) − c

(k) =

ω(k)



(k).

(4.12)

Therefore,

(k) = 0.5





(k) +

ω(k)



(k)



and

(k) = 0.5





(k) −

ω(k)



(k)



For one propagation direction we set c

= 0, then



(k) + k



(k)/ω(k) = 0,

which implies the following relation between each amplitude of the initial condi-

tion for

and the amplitude of the initial condition for η:



(k) = −

ω(k)



(k), k  0.

We use this relation to provide the initial condition for

(by means of an FFT)

for η to propagate in only one direction. Moreover, c

(k) =



(k) = ˆη(k, 0). The

exact solution is

ˆη(k, t) =













(0), k = 0,



(k) exp



−iω(k)t



, k  0.

Although the amplitude of each mode is preserved as time advances, each com-

ponent of the wavetrain travels with its own phase velocity ω(k)/k. As a result,

an initial Gaussian shape will disperse into an oscillatory train as shown in the

numerical experiments. The exact solution can be employed to test the numerical

scheme precision and stability properties. We consider the following Gaussian

function as the initial perturbation of the interface,

(ξ) = 0.5e

−a(ξ−4π)

/128

with a = 50. The LFM parameters are ℓ = 16π, ρ

= 1, ρ

= 2, β = 0.05, α = 0.

The numerical parameters are N = 256, ∆ξ = 2ℓ/N = 0.3927, ∆t = ∆ξ = 0.3927 .

See Fig. 4.7 were the numerical solution for t = 78.5398 is depicted together with

the initial condition. The numerical solution shows very good agreement with the

exact solution. The same experiment with a ﬁner grid, N = 512, ∆ξ = 2ℓ/N =

0.19635, ∆t = ∆ξ = 0.19635 is depicted in Fig. 4.8. The result is the same as the

one obtained in the ﬁrst experiment, which indicates convergence of the numerical

method.

0 20 40 60 80 100 120

−0.2

−0.1

0.1

0.2

0.3

0.4

0.5

Figure 4.7: Pulse propagating over a ﬂat bottom in the linear dispersive regime.

Dotted line: numerical solution for the LFM using RK4 for t = 78.5398 and

N = 256, dashed line: initial condition, solid line: exact solution.

0 20 40 60 80 100 120

−0.2

−0.1

0.1

0.2

0.3

0.4

0.5

Figure 4.8: Pulse propagating over a ﬂat bottom in the linear dispersive regime.

Dotted line: numerical solution for the LFM using RK4 for t = 78.5398 and

N = 512, dashed line: initial condition, solid line: exact solution.

4.4 Periodic topography experiments

When the bottom is ﬂat, the topography dependent coeﬃcient M(ξ) is identically

one. For the time being, we avoid the computation of M(ξ) from the variable depth

bottom, which can be costly even using Driscoll’s package [10]. Let us assume

that it is a function of the form M(ξ) = 1+ n(ξ) where n(ξ) describes periodic ﬂuc-

tuations. This choice is not far from the real coeﬃcient that comes from mapping

a periodic piecewise linear topography, see [21, 24, 22]. This strategy will prove

useful for testing the models and observing the phenomena we are interested in.

Example 4.4. As a ﬁrst example of a rough bottom let us consider a periodic

slowly-varying coeﬃcient M(ξ) deﬁned on the domain [0, 16π] as

M(ξ) =











1 + 0.5sin(5ξ), for 6π ≤ ξ ≤ 12π,

1, elsewhere.

The bottom irregularities are located in the region 6π ≤ ξ ≤ 12π. There are 15

oscillations. The period of the bottom irregularities is l = 1.2566. The initial

perturbation of the interface is the Gaussian function

(ξ) = 0.5e

−a(ξ−π)

/64

with a = 200, therefore its eﬀective width is L = 2.4 and the ratio inhomo-

geneities/wavelength is about 0.5236. For the mean velocity

we choose the

corresponding initial condition that gives one propagation direction for the LFM

with β = 0.0001, α = 0, as done in Section 4.3. See Fig. 4.9 where the numerical

solution for t = 36.3247 is depicted together with the solution for the ﬂat bottom

0 10 20 30 40 50 60

−0.3

−0.2

−0.1

0.1

0.2

0.3

0.4

0.5

Figure 4.9: Pulse propagating over a synthetic periodic slowly-varying topogra-

phy. Dotted line: numerical solution for the LCM using RK4 for t = 36.3247 and

N = 1024, dashed line: initial condition, solid line: ﬂat bottom exact solution.

and the initial condition. The other parameters are ρ

= 1, ρ

= 2, N = 1024,

∆ξ = 2ℓ/N = 0.0491, ∆t = ∆ξ = 0.0491.

A detailed analysis of Fig. 4.9 shows that twice the period of the bottom os-

cillations (2.5133) is in very good agreement with the reﬂected wavelength, as

expected from Bragg’s phenomenon theory [12]. See Fig. 4.10 where vertical

bars marking spatial intervals of size 2.5133 fall together with the end of each

period of the reﬂected signal. A comparison between the solutions for the ﬂat and

periodic bottoms suggests that the attenuation in the wave amplitude is mainly

due to the dispersive term. It is also patent that the oscillations behind the pulse

correspond to the reﬂected wave due to the topography.

Example 4.5. The present example adds the nonlinearity ingredient to the pre-

vious example. We consider again the periodic slowly-varying coeﬃcient M(ξ)

0 10 20 30 40 50 60

−0.3

−0.2

−0.1

0.1

0.2

0.3

0.4

0.5

Figure 4.10: Pulse propagating over a synthetic periodic slowly-varying topogra-

phy. Dashed line: initial condition. Dotted line: numerical solution for the LCM

using RK4 for t = 36.3247 and N = 1024, vertical bars mark spatial intervals of

size 2.5133 that fall together with the end of each period of the reﬂected signal.

deﬁned on the domain [0, 16π] as

M(ξ) =











1 + 0.5sin(5ξ), for 6π ≤ ξ ≤ 12π,

1, elsewhere.

The initial perturbation of the interface is the Gaussian function

(ξ) = 0.5e

−a(ξ−π)

/64

with a = 200 and eﬀective width L = 2.4. The ratio inhomogeneities/wavelength

is about 0.5236. The physical parameters are ρ

= 1, ρ

= 2, β = 0.0001, α =

0.01. We employ N = 1024, ∆ξ = 2ℓ/N = 0.0491, ∆t = ∆ξ = 0.0491. In Fig. 4.11

the numerical solution for t = 32.3977 is depicted together with the exact solution

0 10 20 30 40 50 60

−0.3

−0.2

−0.1

0.1

0.2

0.3

0.4

0.5

Figure 4.11: Pulse propagating over a synthetic periodic slowly-varying topogra-

phy. Dotted line: numerical solution for the WNCM using RK4 for t = 32.3977

and N = 1024, dashed line: initial condition, solid line: exact solution for the

LFM.

for the LFM and the initial condition. Again, twice the period of the bottom

oscillations (2.5133) is in very good agreement with the reﬂected wavelength. In

Fig. 4.12 vertical bars marking spatial intervals of size 2.5133 fall together with

the end of each period of the reﬂected signal.

Example 4.6. Let us consider now a periodic rapidly-varying coeﬃcient M(ξ)

deﬁned on the domain [0, 16π] as

M(ξ) =











1 + 0.5sin(15ξ), for 6π ≤ ξ ≤ 12π,

1, elsewhere.

The bottom irregularities are located in the region 6π ≤ ξ ≤ 12π. The period of

the bottom irregularities is l = 0.4189. The initial perturbation of the interface is

0 10 20 30 40 50 60

−0.3

−0.2

−0.1

0.1

0.2

0.3

0.4

Figure 4.12: Pulse propagating over a synthetic periodic slowly-varying topogra-

phy. Dotted line: numerical solution for the WNCM using RK4 for t = 32.3977

and N = 1024, vertical bars mark spatial intervals of size 2.5133 that fall together

with the end of each period of the reﬂected signal.

now the Gaussian function

(ξ) = 0.5e

−a(ξ−2π)

/64

with a = 50, therefore its eﬀective width is L = 4.8 and the ratio inhomo-

geneities/wavelength is about 0.0873. For the mean velocity

we choose the

corresponding initial condition to ensure one propagation direction for the LFM

with β = 0.0001, α = 0, as done in Section 4.3. See Fig. 4.13 where the numerical

solution for t = 35.3429 is depicted together with the solution for the ﬂat bottom

and the initial condition. The other parameters are ρ

= 1, ρ

= 2, N = 1024,

∆ξ = 2ℓ/N = 0.049087, ∆t = ∆ξ = 0.049087. Note that the solution is very

similar to that of the ﬂat bottom case. The wave is not modiﬁed by the rapidly-

varying topography and no reﬂections are generated. Only the propagation speed

0 10 20 30 40 50 60

−0.1

0.1

0.2

0.3

0.4

0.5

0.6

Figure 4.13: Pulse propagating over a synthetic periodic rapid-varying topogra-

phy. Dotted line: numerical solution for the LCM using RK4 with N = 1024 for

t = 35.3429, dashed line: initial condition, solid line: ﬂat bottom exact solution.

should be slightly decreased as predicted in Rosales and Papanicolaou [26]. But

this change is only noticeable over very large distances.

Example 4.7. Let us add the nonlinearity ingredient (α = 0.005) to the previous

example. We consider again the periodic rapidly-varying topography deﬁned in

Example 4.6, together with the same Gaussian shape of eﬀective width L = 4.8.

Therefore the ratio inhomogeneities/wavelength is kept at 0.0873. In Fig. 4.14 the

numerical solution for t = 35.3429 is depicted together with the exact solution for

the LFM and the initial condition. The other parameters are β = 0.0001, ρ

= 1,

= 2, N = 1024, ∆ξ = 2ℓ/N = 0.049087, ∆t = ∆ξ = 0.049087. Again,

the solution is very similar to that of the LFM. The wave is not modiﬁed by the

rapidly-varying topography and no reﬂections are generated.

0 10 20 30 40 50 60

−0.1

0.1

0.2

0.3

0.4

0.5

0.6

Figure 4.14: Pulse propagating over a synthetic periodic rapid-varying topogra-

phy. Dotted line: numerical solution for the WNCM using RK4 with N = 1024 for

t = 35.3429, dashed line: initial condition, solid line: ﬂat bottom exact solution.

4.5 Computing solitary waves solutions

Now we present two examples of internal solitary waves from the Regularized

ILW equation evolving according to the WNFM. That is, we take as initial condi-

tion for the WNFM a solitary wave from its unidirectional reduction. We expect

the wave to behave almost like a solitary wave. In particular, the balance between

nonlinearity and dispersion should be maintained and the wave should travel with-

out a signiﬁcant change of shape. The velocity of propagation should be similar

to that in the ILW equation. The numerical solutions are obtained by the RK4

numerical solver for the WNFM with ρ

= 1, ρ

= 2, β = 0.0001, α = 0.01,

N = 256, ∆ξ = 2ℓ/N = 0.1963, ℓ = 8π, ∆t = ∆ξ = 0.1963.

Example 4.8. In Fig. 4.15 the evolution of an approximate solitary wave solution

is shown. As initial condition for η, the parameters θ = π/8, a = −0.09953

0 10 20 30 40 50 60

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

Figure 4.15: Numerical solution of the WNFM with β = 0.0001, α = 0.01 for the

propagation of a single solitary wave.

and λ = 0.3323 are selected in the solitary wave solution of Eq. (2.42);

taken to be the corresponding dispersive solution for one propagation direction,

see Section 4.3.

The expected behaviour of the wave is captured by the numerical method for

long times as shown in Fig. 4.16. The pulse propagates with an approximate

velocity of 0.9884 in conformity with its propagation velocity c = 0.9961 in the

Regularized ILW equation (2.42). The shape of the solitary wave is preserved for

long times as shown in Fig. 4.17. The error between the initial condition and the

solution that returns to the original position at approximate time t = 50.8545 is

0.0047. Taking into account that the choice of

is an approximation from the

linear case (α = 0), the result is satisfactory.

Example 4.9. In Fig. 4.18 the ﬁssion of a single approximate solitary wave so-

lution is simulated. As initial condition for η, the same parameters θ = π/8,

a = −0.099536 and λ = 0.3323 are used, while

= 0. We observe two waves

−0.05

0.05

0.1

−η

Figure 4.16: Propagation of a single solitary wave until t = 50.8545.

0 10 20 30 40 50 60

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

−η

Figure 4.17: A single solitary wave, dashed line: initial condition, dotted line:

numerical solution for t = 50.8545.

0.02

0.04

0.06

0.08

0.1

−η

Figure 4.18: Numerical solution of the WNFM for the ﬁssion of a single solitary

wave.

traveling in opposite directions with approximate speed 0.9846. When they coin-

cide (overlap) in space and time, the initial condition is recovered with an error

of 9.8510 × 10

−4

, see Fig. 4.19. This behaviour of the wave is observed for long

times as shown in Fig. 4.20.

The examples in this chapter suggest that the model proposed in Chapter 2

can be implemented numerically and that its basic qualitative properties are well

captured by the numerical solutions.

0 10 20 30 40 50 60

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

−η

Figure 4.19: Dashed line: initial condition, dotted line: initial wave recovered at

t = 51.0509, error = 9.8510 × 10

−4

100

120

−0.05

0.05

0.1

−η

Figure 4.20: Fission of a single solitary wave until t = 102.1018.

Conclusions and future work

In the present work, a one-dimensional strongly nonlinear variable coeﬃcient

Boussinesq-type model for the evolution of internal waves in a two-layer sys-

tem is derived. The regime considered is a shallow water conﬁguration for the

upper layer and an intermediate depth for the lower layer. The bottom has an arbi-

trary, not necessarily smooth nor single-valued proﬁle generalizing the ﬂat bottom

model derived in [8]. This arbitrary topography is dealt with by performing a con-

formal mapping as in [24]. In the unidirectional propagation regime the model

reduces to an ILW equation when a slowly varying topography is assumed. The

adjustment for the periodic wave case and its computational implementation are

also performed. We study the interaction of internal waves with periodic bottom

proﬁles and the evolution of approximate solitary wave solutions. The expected

qualitative behaviour is captured. A higher-order reduced one-dimensional model

is also obtained, though it has not being implemented computationally yet. Both

reduced models have dispersion relations that reproduce correctly the limit from

the Euler equations for the shallow water (long waves) regime in the upper layer.

The higher-ordermodel, by taking into account the nonhydrostatic correction term

for the pressure, approximates better the full dispersion relation. It will be inter-

esting to study numerically the behaviour of the weakly nonlinear higher-order

model regarding this property.

The work initiated here points out some lines of research. We intend to use the

strongly nonlinear model and the weakly nonlinear higher order model to study

the interaction of large amplitude internal waves with multiscale topography pro-

ﬁles. The refocusing and stabilization of solitary waves for the large levels of

nonlinearity allowed by these models is the goal of current research. A diﬃculty

arises to this end, due to the nonlocal dispersive term: it is not easy to ﬁnd ex-

plicit solitary wave solutions for the strongly nonlinear system, and the solitary

waves suggested by the weakly one-directional lower-order theory steepen with

higher nonlinearity. In [8] the strategy to overcome this diﬃculty is to use the

weakly one-directional solitary wave from the ILW equation as an initial guess

for the solitary wave proﬁle. Then ﬁnd a numerical solitary wave iteratively via

the Newton-Raphson method. The proﬁle obtained will serve as initial data for

the propagating model.

Appendix A

Approximation for the horizontal

derivatives at the unperturbed

interface

We want to justify the use of

√

M(ξ )



ξ)



instead of



√

M(ξ )



ξ)



the substitution of η

in system (2.30), in the case of slowly varying topography.

Hence we identify

M(ξ )



ξ)



as the tangential derivative of the solution of

the Laplace equation with Neumann conditions deﬁned in the auxiliary problem











+ Φ

= 0, on −

h(εx) ≤ z ≤ 0,

, at z = 0,

− ε

′

(εx)Φ

= 0, at z = −

h(εx),

(A.1)

where ε is the small parameter deﬁned at the beginning of this text as L/l. It was

already shown that

(x, 0, t) =

M(ξ)



ξ)



via conformal mapping. Now we seek an approximation of Φ

(x, 0, t), our term of

interest, where t is kept frozen. Let ω(x, z) = Φ

(x, z, t). The tangential derivative

of ω at z = 0 is our goal. From Eqs. (A.1), ω satisﬁes











+ ω

= 0, on −

h(εx) ≤ z ≤ 0,

, at z = 0,

− ε

′

(εx)ω

− ε

′′



εx



ω = 0, at z = −

h(εx).

A conformal mapping taking a ﬂat strip into the corrugated strip above trans-

forms this problem into











ξξ

+ ω

ζζ

= 0, on −

≤ ζ ≤ 0,

= M(ξ)



x(ξ, 0), t



, at ζ = 0,

− ε

′′



ε x(ξ, 0)



ω = 0, at ζ = −

Note that a mixed boundary condition (Robin condition) is set on ζ = −

instead of a Neumann condition as in all previous problems. Because of the Robin

condition, the tangential derivative ω

(ξ, 0, t) is no longer T



ξ)



ξ, 0), t





but it can be approximated by this term up to a certain order in ε. Following a

perturbation approach, consider

ω = ω

+ εω

+ ε

+ ···

then ω

satisﬁes the Neumann problem











ξξ

+ ω

ζζ

= 0, on −

≤ ζ ≤ 0,

= M(ξ)



x(ξ, 0), t



, at ζ = 0,

= 0, at ζ = −

with tangential derivative ω

(ξ, 0) = T



ξ)



ξ, 0), t





(ξ). Thus ω

/M(ξ). The subsequent term is ω

= 0 because both boundary conditions are

homogeneous to O(ε). Next ω

satisﬁes











ξξ

+ ω

ζζ

= 0, on −

≤ ζ ≤ 0,

= 0, at ζ = 0,

−

′′



ε x(ξ, 0)



= 0, at ζ = −

Therefore

ω = ω

+ O(ε

). (A.2)

If we establish a relation between α, β and ε of the type ε

= O(β

), with q ≥

Eq. (A.2) leads to



M(ξ)



ξ)





M(ξ)



ξ)



+ O







which justiﬁes the approximation done in Eq. (2.39).

Appendix B

The Dirichlet-to-Neumann operator

Due to the importance of the Dirichlet-to-Neumann operator (DtN

) in modelling

water waves we pointed here its relation with the Hilbert transform on the ﬂat strip

used througout this work.

Given the problem











υ υ

+ z

ω ω

= 0, −h ≤

ω ≤ 0,

= g(υ), ω = 0,

= 0,

ω = −h,

the operator T returns the tangential derivative z

(

υ, 0) from the Neumann data

(υ, 0) = g(υ), that is

T[g] = z

(

υ, 0).

Therefore, its inverse T satisﬁes

g(υ) = z

(

υ, 0) = T[z

(

υ, 0)]. (B.1)

The inverse T of T has a symbol −itanh

(

)

. Alternatively it can be deﬁned

by the principal value integral

T[f] = −



f(v)

sinh



π(v − υ)/2h



dv. (B.2)

The DtN

for the problem











υ υ

+ z

ω ω

= 0, −h ≤

ω ≤ 0,

υ, 0) = f(υ),

υ, −h) = 0,

returns the normal derivative (Neumann condition) from the Dirichlet data. In

light of (B.1) we have

DtN

[f] = T[f

′

It means that the Dirichlet-to-Neumann operator applies one diﬀerentiation plus

the inverse operator

T to the Dirichlet data in order to obtain the Neumann con-

dition.

Appendix C

The periodic counterpart of the

operator T

Due to the nonlocal deﬁnition of the Hilbert transform, it is necessary to redeﬁne

it when restricting our problem to the periodic domain for numerical implemen-

tations. To that end we keep in mind its geometrical interpretation, that is, the

operator that takes a harmonic function’s normal derivative at the boundary and

transforms it into its tangential derivative at the boundary.

Consider the following problem for periodic ξ ∈ Π[0, 2ℓ],











ξξ

+ φ

ζζ

= 0, −h

/L ≤ ζ ≤ 0, 0 ≤ ξ ≤ 2ℓ,

φ(0, ζ) = φ(2ℓ, ζ),

(ξ, 0) = g(ξ),

(ξ, −h

/L) = 0.

Set ξ = πξ/ℓ. In the new coordinates φ(ξ, ζ) = φ(ξ, ζ) satisﬁes















ξ ξ

ζζ

= 0, −h

/L ≤ ζ ≤ 0, 0 ≤ ξ ≤ 2π,

φ(0, ζ) = φ(2π, ζ),

(ξ, 0) = g(ξ),

(ξ, −h

/L) = 0.

(C.1)

We consider the Fourier Series in

ξ ∈ [0, 2π] with its coeﬃcients given by

f(k) =

2π



ξ)e

−ik

ξ,

ξ) =

2π

∞



k=−∞

f(k)e

(C.2)

The Discrete Fourier Transform (DFT) in

ξ is

f(k) = ∆



j=1

−ik

= j∆ξ, ∆ξ =

2π

exactly the same one used by Trefethen [28] together with the inverse

) =

2π

N/2



k=−N/2+1

f(k)e

, j = 1, . . . , N,

where k ∈ {−N/2 + 1, . . . , N/2} because in the discrete domain e

ikj∆

= e

ikj2π/N

there is no diﬀerence for k = k

mod N.

Apply the Fourier Transform to problem (C.1):











(−ik)



ℓ



φ +

ζζ

= 0, on −h

/L ≤ ζ ≤ 0,

(k, 0) =

g(k),

(k, −h

/L) = 0.

(C.3)

The solution for problem (C.3) is

φ(k, ζ) =

g(k)

πk/ℓ

cosh



kπ

ℓ



ζ +



sinh



kπ

ℓ



, k  0.

Therefore, for C

initial Neumann data g,

φ(ξ, ζ) =

2π

∞



k=−∞

φ(k, ζ)e

2π

∞



k=−∞

k0

g(k)

πk/ℓ

cosh



kπ

ℓ



ζ +



sinh



kπ

ℓ



φ(0)

2π

the convergence is uniform in 0 ≤

ξ ≤ 2π and also in −h

/L ≤ ζ ≤ 0 since



cosh



kπ

ℓ



ζ +



sinh



kπ

ℓ





≤



cosh



kπ

ℓ



sinh



kπ

ℓ





coth



ℓ





, (C.4)

for all integer k  0.

We want φ

; returning to the original variables

φ(ξ, ζ) =

2π

∞



k=−∞

k0

g(k)

πk/ℓ

cosh



kπ

ℓ



ζ +



sinh



kπ

ℓ



ikξπ/ℓ

φ(0)

2π

Taking ξ -derivatives,

(ξ, ζ) =

2π

∞



k=−∞

k0

g(k)

icosh



kπ

ℓ



ζ +



sinh



kπ

ℓ



ikξπ/ℓ

The tangential derivative at the boundary is obtain making ζ → 0:

(ξ, 0) =

2π

∞



k=−∞

k0

g(k)

icosh



kπ

ℓ



sinh



kπ

ℓ



ikξπ/ℓ

The convergence is still uniform because of Eq. (C.4).

Therefore,

[0,2ℓ]

[f](ξ) =

2π

∞



k=−∞

k0

icoth



kπ

ℓ



ikξπ/ℓ

f(k),

where

f(k) =

2π



f(ξ)e

−ik

ξ, ξ = πξ/ℓ,

that is, the Fourier coeﬃcients in Π[0, 2π].

It is also convenient to write the composition of one spatial derivative with

the Hilbert transform T

[0,2ℓ]

[

]

because they always come together in the models

considered here,

[0,2ℓ]

[f]

(ℓ

ξ/π) =

2π

∞



k=−∞

k0

−

kπ

ℓ

coth



kπ

ℓ



f(k).

Finally, for the discretization of the periodic domain (ignoring the aliasing

eﬀect) we have that

[0,2ℓ]

[f]

(

) ≈

2π

N/2



k=−N/2+1

k0

−

kπ

ℓ

coth



kπ

ℓ



f(k), j = 1, . . . , N,

where

f(k) = ∆ξ



j=1

f(ξ

−ik

In conclusion, the symbol of the operator T

[0,2ℓ]

[

]

(that is the composition of

one spatial derivative with the Hilbert transform) in the new coordinate

ξ is

−

kπ

ℓ

coth



kπ

ℓ



Bibliography

[1] Artiles, W. & Nachbin, A., 2004. “Nonlinear evolution of surface grav-

ity waves over highly variable depth,” Physical Review Letters, vol. 93,

pp. 234501–1–234501–4.

[2] Ascher, U. M. & Petzold, L. R., 1988. Computer Methods for Ordinary Dif-

ferential Equations and Diﬀerential–Algebraic Equations, SIAM.

[3] Benjamin, T. B., 1967. “Internal waves of permanent form of great depth,”

Journal of Fluid Mechanics, vol. 29, pp. 559–592.

[4] Benjamin, T. B., Bona, J. L. & Mahony, J. J., 1972. “Model equations for

long waves in nonlinear dispersive systems,” Philosophical Transactions of

the Royal Society of London. Series A, Mathematical and Physical Sciences,

vol. 272, No. 1220, pp. 47–78.

[5] Camassa, R. & Levermore, C. D., 1997. “Layer-mean quantities, local con-

servation laws, and vorticity,” Physical Review Letters, vol. 78, pp. 650–653.

[6] Choi, W., & Camassa, R., 1996. “Long internal waves of ﬁnite amplitude,”

Physical Review Letters, vol. 77, pp. 1759–1762.

[7] Choi, W., & Camassa, R., 1996. “Weakly nonlinear internal waves in a two-

ﬂuid system,” Journal of Fluid Mechanics, vol. 313, pp. 83–103.

[8] Choi, W., & Camassa, R., 1999. “Fully nonlinear internal waves in a two-

ﬂuid system,” Journal of Fluid Mechanics, vol. 396, pp. 1–36.

[9] Davis, R. E., & Acrivos, A., 1967. “Solitary internal waves in deep water,”

Journal of Fluid Mechanics, vol. 29, pp. 593–607.

[10] Driscoll, T., Schwarz-Christoﬀel toolbox for Matlab,

http://www.math.udel.edu/˜driscoll/software.

[11] Green, A. E. & Naghdi, P. M., 1976. “A derivation of equations for wave

propagation in water of variable depth,” Journal of Fluid Mechanics, vol. 78,

pp. 237–246.

[12] Guazzeli, E., Rey, V. & Belzons, M., 1992. “Higher order Bragg reﬂection

of gravity surface waves by periodic beds,” Journal of Fluid Mechanics,

vol. 245, pp. 301–317.

[13] Jo, T.-C. & Choi, W., 2002. “Dynamics of strongly nonlinear internal solitary

waves in shallow water,” Studies in Applied Mathematics, vol. 109, pp. 205–

227.

[14] Joseph, R. I., 1977. “Solitary waves in ﬁnite depth ﬂuid,” Journal of

Physics A, vol. 10, pp. L225–L227.

[15] Keener, J. P., 2000. Principles of Applied Mathematics: Transformation and

Approximation, Westview Press.

[16] Koop, C. G. & Butler, G., 1981. “An investigation of internal solitary waves

in a two-ﬂuid system,” Journal of Fluid Mechanics, vol. 112, pp. 225–251.

[17] Kubota, T., Ko, D., & Dobbs, L., 1978. “Propagation of weakly non linear

internal waves in a stratiﬁed ﬂuid of ﬁnite depth,” AIAA Journal Hydrody-

namics, vol. 12, pp. 157–165.

[18] Lamb, H., 1932. Hydrodynamics, Dover.

[19] Matsuno, Y., 1992. “Nonlinear evolutions of surface gravity waves on ﬂuid

of ﬁnite depth,” Physical Review Letters, vol. 69, pp. 609–611.

[20] Matsuno, Y., 1993. “A uniﬁed theory of nonlinear wave propagation in two-

ﬂuid systems,” Journal of the Physical Society of Japan, vol. 62, pp. 1902–

1916.

[21] Mu

noz, J. C. & Nachbin, A., 2004. “Dispersive wave attenuation due to

orographic forcing,” SIAM Journal of Applied Mathematics, vol. 64, Issue 3,

pp. 977–1001.

[22] Mu

noz, J. C. & Nachbin, A., 2005. “Stiﬀ microscale forcing and solitary

wave refocusing,” SIAM Multiscale Modeling and Simulation, vol.3, issue 3,

pp. 680–705.

[23] Mu

noz, J. C. & Nachbin, A., 2006. “Improved Boussinesq-type equations

for highly-variable depths,” IMA Journal of Applied Mathematics, vol. 71,

pp. 600–633.

[24] Nachbin, A., 2003. “A terrain-following Boussinesq system,” SIAM Journal

on Applied Mathematics, vol. 63, pp. 905–922.

100

[25] Ono, H., 1975. “Algebraic solitary waves in stratiﬁed ﬂuids,” Journal of the

Physical Society of Japan, vol. 39, pp. 1082–1091.

[26] Rosales, R. R. & Papanicolaou, G. C., 1983. “Gravity waves in a channel

with a rough bottom,” Studies in Applied Mathematics, vol. 68, pp. 89–102.

[27] Su, C. H. & Gardner, C. S., 1969. “Korteweg-de Vries Equation and Gen-

eralizations. III. Derivation of the Korteweg-de Vries Equation and Burgers

Equation,” Journal of Mathematical Physics, vol. 10, issue 3, pp. 536–539.

[28] Trefethen, L. N., 2000. Spectral Methods in Matlab, SIAM.

[29] Wei, G. & Kirby, J., 1995. “Time-dependent numerical code for extended

Boussinesq equations,” Journal of Waterway, Port, Coastal and Ocean En-

gineering, vol. 121 pp. 251–261.

[30] Whitham, G. B., 1974. Linear and nonlinear waves, John Wiley, New York–

London–Sidney.

[31] Wu, T. Y., 1981. “Long waves in ocean and coastal waters,” Journal of the

Engineering Mechanics Division ASCE, vol. 107, No. 3, pp. 501–522.

101

Livros Grátis
( http://www.livrosgratis.com.br )
 
Milhares de Livros para Download:
 
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas

Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
 
 