Lecture Notes

“Modeling Lab Course”

If you find any typos or other errors, please let me know!

Prof. Dr. Hendrik Ranocha

May 8, 2025
Contents
1 Introduction 3
1.1 Scalar nonlinear dispersive transport equations . . . . . . . . . . . . . . . 3
1.1.1 Korteweg-de Vries (KdV) equation . . . . . . . . . . . . . . . . . . 4
1.1.2 The Benjamin-Bona-Mahony (BBM) equation . . . . . . . . . . . . 4
1.1.3 The nonlinear Schrödinger (NLS) equation . . . . . . . . . . . . . 5
1.1.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.1.5 First projects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Dispersion relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2 Projects 8
2.1 Scalar models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Systems with flat bathymetry in 1D . . . . . . . . . . . . . . . . . . . . . . 9
2.2.1 BBM-BBM equations . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.2 SGN equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Systems with flat bathymetry in 1D: comparison with experimental data 10

3 Deriving shallow water models from first principles 12

3.1 The incompressible Euler equations . . . . . . . . . . . . . . . . . . . . . 12
3.1.1 Irrotational flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.2 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.3 Non-dimensionalization . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2 Deriving simplified models . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2.1 Handling the pressure term . . . . . . . . . . . . . . . . . . . . . . 19
3.2.2 Handling the integrals . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3 The Serre-Green-Naghdi (SGN) equations . . . . . . . . . . . . . . . . . . 21
3.4 Boussinesq approximations and the BBM-BBM system . . . . . . . . . . 22
3.5 The shallow water (Saint-Venant) equations . . . . . . . . . . . . . . . . . 25
3.6 Unidirectional models: KdV and BBM . . . . . . . . . . . . . . . . . . . . 25

2
1 Introduction
In this course, we will study study transport equations that can be used to model
water waves with nonlinear and dispersive effects. You have already seen nonlinear
(hyperbolic) transport equations in the course introducing numerical methods for PDE,
e.g., the classical 1D shallow water equations

ℎ ℎ𝑣
   
𝜕𝑡 + 𝜕𝑥 =0
ℎ𝑣 ℎ𝑣 + 21 𝑔 ℎ 2
2
|{z} | {z }
=𝑢 = 𝑓 (𝑢)

for a flat bathymetry, where 𝑡 is the time, 𝑥 the space coordinate, ℎ is the water height,
𝑣 its speed, ℎ𝑣 the discharge, and 𝑔 the gravitational acceleration (9.81 m/s2 ). In this
course, we will focus on extensions of this hyperbolic system including

• variable bottom topography (bathymetry) 𝑏

• dispersive effects

• multiple space dimensions

To do so, we will derive several dispersive water wave models from the incompressible
Euler equations to give you an idea about the physical modeling process. Before doing
so, we will first present some basic scalar models of nonlinear dispersive transport
equations that illustrate some main ideas and behaviors.

1.1 Scalar nonlinear dispersive transport equations

We will consider three propotypes of scalar nonlinear dispersive transport equations
in one space dimension:

• the Korteweg-de Vries (KdV) equation named after

D. J. Korteweg and G. De Vries (1895). “On the change of form of long
waves advancing in a rectangular canal, and on a new type of long
stationary waves.” In: The London, Edinburgh, and Dublin Philosophical
Magazine and Journal of Science 39.240, pp. 422–443. doi: 10. 1080 /
14786449508620739

• the Benjamin-Bona-Mahony (BBM) equation named after

T. B. Benjamin, J. L. Bona, and J. J. Mahony (1972). “Model equations
for long waves in nonlinear dispersive systems.” In: Philosophical Trans-
actions of the Royal Society of London. Series A, Mathematical and Physical
Sciences 272.1220, pp. 47–78. doi: 10.1098/rsta.1972.0032

• the nonlinear Schrödinger (NLS) equation

3
1.1.1 Korteweg-de Vries (KdV) equation
A well-known dispersive transport equation is the Korteweg-de Vries (KdV) equation
𝑢2
𝜕𝑡 𝑢 + 𝜕𝑥 𝑢 + 𝜕𝑥
+ 𝜕𝑥3 𝑢 = 0, (1.1)
2
where we scaled physical parameters to obtain a clean form of the equations. This
equation conserves the quadratic invariant
∫
1
𝒥 (𝑢) = 𝑢(𝑡, 𝑥)2 d𝑥
2
since
𝑢2
 
1 2
𝜕𝑡 𝑢 = 𝑢𝜕𝑡 𝑢 = −𝑢𝜕𝑥 𝑢 − 𝑢𝜕𝑥 − 𝑢𝜕𝑥3 𝑢
2 2
𝑢2 𝑢3
 
1
= −𝜕𝑥 − 𝜕𝑥 − 𝜕𝑥 𝑢𝜕𝑥2 𝑢 − (𝜕𝑥 𝑢)2
2 3 2
for smooth solutions, which is of the general form of a conservation law
𝜕𝑡 𝑃(𝑢) + 𝜕𝑥 𝑄(𝑢) = 0.
You can check that the KdV equation (1.1) has solitary wave solutions of the form
√ !2
𝑐
𝑢(𝑡, 𝑥) = 3𝑐 sech 𝑥 − (𝑐 + 1)𝑡 − 𝑥 0 ,


2

where 𝑐 > 0 is a parameter determining the wave speed, 𝑥 0 ∈ R determines the initial
position at time 𝑡 = 0, and sech(𝑦) = 1/cosh(𝑦) is the hyperbolic secant.

1.1.2 The Benjamin-Bona-Mahony (BBM) equation

The Benjamin-Bona-Mahony (BBM) equation
𝑢2
− 𝜕𝑡 𝜕𝑥2 𝑢 = 0
𝜕𝑡 𝑢 + 𝜕𝑥 𝑢 + 𝜕𝑥 (1.2)
2
is closely related to the KdV equation; the only difference is in the dispersive term,
where the third derivative in space has been replaced by a mixed space-time derivative.
We have again scaled the physical parameters to obtain a clean form of the equations.
This equation conserves the quadratic invariant
∫
1 2 
𝒥 (𝑢) = 𝑢 + (𝜕𝑥 𝑢) ,
2
2
which can be shown by multiplying the BBM equation (1.2) by 𝑢 and computing
0 = 𝑢𝑢𝑡 + 𝑢𝑢𝑥 + 𝑢(𝑢 2 /2)𝑥 − 𝑢𝑢𝑡𝑥𝑥
= (𝑢 2 /2)𝑡 + (𝑢 2 /2)𝑥 + (𝑢 3 /3)𝑥 + 𝑢𝑥 𝑢𝑡𝑥 − (𝑢𝑢𝑡𝑥 )𝑥
   
1 2 1 1 1
= 𝜕𝑡 𝑢 + (𝜕𝑥 𝑢)2 + 𝜕𝑥 𝑢 2 + 𝑢 3 − 𝑢𝑢𝑡𝑥
2 2 2 3
for smooth solutions, which is of the general form of a conservation law. You can check
that the BBM equation (1.2) has solitary wave solutions of the form
!2
𝑐
r
1
𝑢(𝑡, 𝑥) = 3𝑐 sech 𝑥 − (𝑐 + 1)𝑡 − 𝑥 0 .


2 1+𝑐

4
1.1.3 The nonlinear Schrödinger (NLS) equation
The (cubic) nonlinear Schrödinger (NLS) equation
i𝜕𝑡 𝑢 + 𝜕𝑥2 𝑢 + |𝑢|2 𝑢 = 0 (1.3)
has a different structure than the KdV and BBM equations. In particular, the solution
𝑢 is complex-valued, i.e., 𝑢(𝑡, 𝑥) ∈ C. This equation conserves the quadratic invariant
∫
𝒥 (𝑢) = |𝑢|2 d𝑥,

which can be shown as follows. First, we compute

𝜕𝑡 |𝑢|2 = 𝜕𝑡 (𝑢𝑢) = 𝑢𝑢𝑡 + 𝑢 𝑡 𝑢.
Next, we rearrange the NLS equation (1.3) to obtain
𝑢𝑡 = i𝑢𝑥𝑥 + i|𝑢|2 𝑢 and 𝑢 𝑡 = −i𝑢 𝑥𝑥 − i|𝑢|2 𝑢.
Thus, we obtain
𝜕𝑡 |𝑢|2 = i𝑢𝑢𝑥𝑥 + i|𝑢|4 − i𝑢𝑢 𝑥𝑥 − i|𝑢|4
= i𝜕𝑥 𝑢𝑢𝑥 − 𝑢 𝑥 𝑢 ,

which is of the general form of a conservation law. You can check that the NLS equation
(1.3) has solitary wave solutions of the form
√ 
 √ 
𝑢(𝑡, 𝑥) = 2𝛼 exp i 𝑐(𝑥 − 𝑥 0 )/2 − (𝑐 2 /4 − 𝛼)𝑡 sech 𝛼(𝑥 − 𝑥 0 − 𝑐𝑡) ,

where 𝑐, 𝛼 > 0 are real parameters determining the wave speed and amplitude, and
𝑥0 ∈ R determines the initial position at time 𝑡 = 0.

1.1.4 Discussion
There are some similarities between these simple scalar models, e.g., there are solitary
wave solutions and (quadratic) invariants. However, not all dispersive wave equations
have obvious and non-trivial invariants. If there are such invariants, they can often give
helpful estimates for analysis and for developing numerical methods. For example, the
quadratic invariant of the BBM equation is just the squared 𝐻 1 norm of the solution,
i.e., it guarantees a certain smoothness of the solution.
Moreover, many dispersive wave equations have solitary wave solutions, but there
are explicit formulae only in very few cases. As you will see in the projects, the
interaction of solitary wave solutions is also different in these models. Some models
allow for soliton solutions, i.e., solitary waves that can interact with each other and
preserve their shape and speed after the interaction.
Next, I have some questions to you to prepare the next introductory lectures:
• Do you know what a Hamiltonian ODE/PDE system is?
• Do you know what a Lagrangian ODE/PDE system is?
• Have you heard of ”calculus of variations”, “variational methods”, or a “varia-
tional principle”?
• Have you heard of the Legendre transform?
• Do you know what a dispersion relation is?
• Have you heard of the group and phase speeds?

5
1.1.5 First projects
Split into two groups, each group shall work on one of the equations described above.
Your tasks are:
Task!
• Find/derive a numerical method to solve the equation.

• Implement the method and test it by performing convergence experiments.

• Do you experience a CFL condition for explicit time integration methods?

• Check conservation of the invariants. Can you find efficient numerical methods
that conserve the invariants? Does this change the error growth in time for long
time simulations with periodic boundary conditions? Consider cases where the
error in time dominates the error in space.

• What happens if you initialize a simulation with multiple solitary waves that are
isolated from each other but start to interact after some time due to different wave
speeds?

You have to present your results for these tasks on Thursday, 2025-04-24 (roughly 10–15
minutes per group) and be prepared to answer questions. We will reserve some time
for questions/discussions/project work on Thursday, 2025-04-17.
After the presentations, the groups will cross-check each other’s results and make
sure that everything is reproducible.
Please remember to send me your GitHub account names so that I can add you to a
private GitHub repository we use for sharing the results.
→ end of lecture 2025-04-15

1.2 Dispersion relations

Next, we will have a look at dispersion relations for the KdV and BBM equations. This
concept comes from Fourier analysis and requires linearization. For example, assume
that we are interested in a small perturbation of a background state solving the KdV
equation. For simplicity, we assume that the background state vanishes. Thus, we
can linearize the KdV equation (1.1) around the background state 𝑢0 = 0 simply by
dropping the quadratic term. Thus, we obtain

𝜕𝑡 𝑢 + 𝜕𝑥 𝑢 + 𝜕𝑥3 𝑢 = 0.


i 𝑘𝑥−𝜔(𝑘)𝑡
A plane wave 𝑢(𝑡, 𝑥) = e solves this equation if and only if
 

i 𝑘𝑥−𝜔(𝑘)𝑡
−i𝜔(𝑘) + i𝑘 − i𝑘 3
e = 0,

i.e., if
𝜔(𝑘) = 𝑘 − 𝑘 3 .
For the BBM equation (1.2), we obtain the linearization around the background state
𝑢0 = 0 as
𝜕𝑡 𝑢 + 𝜕𝑥 𝑢 − 𝜕𝑡 𝜕𝑥2 𝑢 = 0.
A plane-wave ansatz yields
 

𝑘𝑥−𝜔(𝑘)𝑡
−i𝜔(𝑘) + i𝑘 − i𝜔(𝑘)𝑘 2 ei = 0,

6
i.e., the dispersion relation
𝑘
𝜔(𝑘) = .
1 + 𝑘2
In general, the phase velocity is

𝜔(𝑘)
𝑐phase (𝑘) = .
𝑘
To see why this is called phase velocity, consider a plane wave



i 𝑘𝑥−𝜔(𝑘)𝑡 i𝑘 𝑥−𝜔(𝑘)𝑡/𝑘
𝑢(𝑡, 𝑥) = e =e .

This is clearly of the form of a traveling harmonic wave with speed 𝜔(𝑘)/𝑘. The group
velocity is more difficult to understand. It is defined as

d𝜔(𝑘)
𝑐group (𝑘) = .
d𝑘
Thus, the KdV equation has a negative phase velocity 𝜔(𝑘)/𝑘 = 1 − 𝑘 2 for 𝑘 > 1.
Moreover, the group velocity d𝜔(𝑘)/d𝑘 = 1 − 3𝑘 2 is unbounded, i.e., fine-scale features
(high-wavenumber components) propagate arbitrarily fast into the negative direction.
In contrast, the BBM equation has a positive and bounded phase velocity 𝜔(𝑘)/𝑘 =
1/(1 + 𝑘 2 ) and a bounded group velocity d𝜔(𝑘)/d𝑘 = (1 − 𝑘 2 )/(1 + 𝑘 2 )2 .

7
2 Projects
This chapter collects the projects that students work on. Please note that the project
work overlaps with the accompanying lecture introducing the modeling process. Thus,
you will have to work on projects described here while we still discuss material from
the other chapters. Some general notes:
We use the private GitHub repository
[Link]

to share the codes. You can also use it to collaborate on your projects.
Your code has to be reproducible. Typically, after each milestone set by the sub-
projects, you will present your results in class. Afterward, the others have to check
your code for reproducibility and report the results in the next meeting/week.
You will work in two groups of three students each.

2.1 Scalar models

This is basically Section 1.1.5. One group works with the KdV equation (1.1), the other
with the BBM equation (1.2). Your tasks are:

• Find/derive a numerical method to solve the equation.

• Implement the method and test it by performing convergence experiments.

• Do you experience a CFL condition for explicit time integration methods?

• Check conservation of the invariants. Can you find efficient numerical methods
that conserve the invariants? Does this change the error growth in time for long
time simulations with periodic boundary conditions? Consider cases where the
error in time dominates the error in space.

• What happens if you initialize a simulation with multiple solitary waves that are
isolated from each other but start to interact after some time due to different wave
speeds?

You have to present your results for these tasks on Thursday, 2025-04-24 (roughly 10–15
minutes per group) and be prepared to answer questions. We will reserve some time
for questions/discussions/project work on Thursday, 2025-04-17.
After the presentations, the groups will cross-check each other’s results and make
sure that everything is reproducible.
Please remember to send me your GitHub account names so that I can add you to a
private GitHub repository we use for sharing the results.

8
2.2 Systems with flat bathymetry in 1D
One group focuses on the BBM-BBM equations
1
𝜂𝑡 + (𝜂 + 𝐷)𝑣 − 𝐷 2 𝜂𝑡𝑥𝑥 = 0,


6 𝑥
1 2
𝑣 𝑡 + 𝑔𝜂 𝑥 + 𝑣𝑣 𝑥 − 𝐷 𝑣 𝑡𝑥𝑥 = 0,
6
formulated equivalently (for 𝜂0 = 0, see below) as
1
ℎ 𝑡 + (ℎ𝑣)𝑥 − 𝑏 2 ℎ 𝑡𝑥𝑥 = 0,
6
1
𝑣 𝑡 + 𝑔𝜂 𝑥 + 𝑣𝑣 𝑥 − 𝑏 2 𝑣 𝑡𝑥𝑥 = 0,
6
the other group on the Serre-Green-Naghdi equations
ℎ 𝑡 + (ℎ𝑣)𝑥 = 0,
   
1
𝑣 𝑡 + 𝑔 ℎ 𝑥 + 𝑣𝑣 𝑥 − ℎ 𝑣 𝑡𝑥 + 𝑣𝑣 𝑥𝑥 − (𝑣 𝑥 )
3 2
= 0,
3ℎ 𝑥
where
• 𝜂(𝑡, 𝑥) is the total water height,
• 𝑣(𝑡, 𝑥) is the (somehow approximated) horizontal velocity of the water,
• 𝑔 is the gravitational acceleration,
• 𝜂0 is the constant still water height (typically assumed to be zero)
• 𝑏 = 𝜂0 − 𝐷 is the constant bathymetry so that 𝐷 is the still water depth,
• ℎ = 𝜂 − 𝜂0 + 𝐷 = 𝜂 − 𝑏 is the total water height above the flat bathymetry.
Additional details depending on the equations are given in the subsection below. Your
tasks are:
• Find/derive a numerical method to solve the equation, focus on periodic bound-
ary conditions for now.
• Implement the method and test it by performing convergence experiments.
• Do you experience a CFL condition for explicit time integration methods?
• Check conservation of the invariants. Can you find efficient numerical methods
that conserve the invariants? Does this change the error growth in time for long
time simulations with periodic boundary conditions? Consider cases where the
error in time dominates the error in space.
• What happens if you initialize a simulation with multiple solitary waves that are
isolated from each other but start to interact after some time due to different wave
speeds?
You have to present your results for these tasks on Thursday, 2025-05-08 (roughly 15
minutes per group) and be prepared to answer questions. After the presentations,
the groups will cross-check each other’s results and make sure that everything is
reproducible.
Remark. Do not (only) implement the methods available in [Link]
— I want to see how other methods behave, too. ⊳

9
2.2.1 BBM-BBM equations
The BBM-BBM equations conserve the energy functional
∫  
1 2 1
𝑔𝜂 + (𝜂 + 𝐷)𝑣 2
2 2
for periodic boundary conditions.
There are no physically relevant known solitary wave solutions in closed form (as far
as I know). A code generating solitary wave profiles numerically to high accuracy will
be provided in the repository.

2.2.2 SGN equations

The SGN equations conserve the (modified) energy functional
∫  
1 2 1 2 1 3
𝑔𝜂 + ℎ𝑣 + ℎ (𝑣 𝑥 )2
2 2 6
for periodic boundary conditions.
The SGN equations have solitary wave solutions of the form

2
ℎ(𝑡, 𝑥) = ℎ ∞ 1 + 𝜀 sech 𝜅(𝑥 − 𝑐𝑡) ,
ℎ∞
 
𝑣(𝑡, 𝑥) = 𝑐 1 − ,
ℎ(𝑡, 𝑥)
where 𝜀 = 𝐴/ℎ∞ , 𝐴 is the solitary wave amplitude,
3𝜀
𝜅2 = , 𝑐 2 = 𝑔 ℎ∞ (1 + 𝜀).
2
4ℎ∞ (1 + 𝜀)
A typical choice is ℎ∞ = 1.

2.3 Systems with flat bathymetry in 1D: comparison with

experimental data
Keep using the same model (BBM-BBM or SGN) in the groups as before. This time, you
have to compare your results to experimental/numerical data based on the following
setups.
• interaction of solitary waves described in Section 5.3 of D. Dutykh, D. Clamond, P.
Milewski, and D. Mitsotakis (2013). “Finite volume and pseudo-spectral schemes
for the fully nonlinear 1D Serre equations.” In: European Journal of Applied Mathe-
matics 24.5, pp. 761–787. doi: 10.1017/S0956792513000168
• dispersive shock wave described in Section 5 of J. P. Pitt, C. Zoppou, and S. G.
Roberts (2022). “Numerical scheme for the generalised Serre-Green-Naghdi
model.” In: Wave Motion 115, p. 103077. doi: 10.1016/[Link].2022.
103077
• dispersive shock wave described in Section 5.4 of J. A. do Carmo, J. Ferreira, L.
Pinto, and G. Romanazzi (2018). “An improved Serre model: Efficient simulation
and comparative evaluation.” In: Applied Mathematical Modelling 56, pp. 404–423.
doi: 10.1016/[Link].2017.12.005

10
 • dispersive shock wave described in Section 11.6 of H. Ranocha and M. Ricchiuto
(Aug. 2024). Structure-preserving approximations of the Serre-Green-Naghdi equations
in standard and hyperbolic form. arXiv: 2408.02665 [[Link]]

• Favre waves described in Section 11.8 (see Figures 12 and 17) of Ranocha and
Ricchiuto (2024)

You will likely have to be creative to obtain the experimental data in some cases —
either by writing emails to some authors or by trying to extract data from plots.
You have to present your results for these tasks on Tuesday, 2025-05-27 and be
prepared to answer questions. After the presentations, the groups will cross-check
each other’s results and make sure that everything is reproducible.

11
3 Deriving shallow water models from
first principles
In this section, we follow
D. Mitsotakis (2013). A simple introduction to water waves. url: https :
//[Link]/hal-00805080
to derive shallow water models from the incompressible Euler equations. The article
also includes a derivation of the incompressible Euler equations from first principles in
mechanics, which we will not discuss here. A more detailed derivation of additional
models and some related discussions can be found in
D. Lannes (2020). “Modeling shallow water waves.” In: Nonlinearity 33.5,
R1. doi: 10.1088/1361-6544/ab6c7c
and the monograph
D. Lannes (2013). The water waves problem. Vol. 188. Mathematical Surveys
and Monographs. Providence, Rhode Island, USA: American Mathematical
Society (AMS)

3.1 The incompressible Euler equations

First, we will derive the incompressible Euler equations for a homogeneous fluid with
gravitational forces following
A. J. Chorin and J. E. Marsden (1993). A mathematical introduction to fluid
mechanics. Vol. 4. Texts in Applied Mathematics. New York: Springer. doi:
10.1007/978-1-4612-0883-9
Assume we can model the motion of a fluid by its density field 𝜚(𝑡, 𝑥) and its velocity
field 𝑣(𝑡, 𝑥), where 𝑡 is the time and 𝑥 is the space coordinate.

Conservation of mass The rate of change of the total mass of the fluid in a volume 𝑉
is given by ∫ ∫
d
𝜚(𝑡, 𝑥) = 𝜕𝑡 𝜚(𝑡, 𝑥).
d𝑡 𝑉 𝑉
Since mass cannot be created or destroyed, we have the balance equation
∫ ∫ ∫
d
𝜚(𝑡, 𝑥) = − 𝜚𝑣 · 𝑛 = − div(𝜚𝑣).
d𝑡 𝑉 𝜕𝑉 𝑉

Thus, if everything is sufficiently smooth to justify the steps above, we obtain the
conservation of mass equation (continuity equation)

𝜕𝑡 𝜚 + div(𝜚𝑣) = 0

in differential form.

12
Conservation of momentum We follow
a fluid particle at position 𝑥(𝑡) with corre-
sponding velocity 𝑥(𝑡) = d𝑡 𝑥(𝑡) = 𝑣 𝑡, 𝑥(𝑡) . Thus, the acceleration of the fluid particle
¤ d

is
d
𝑥(𝑡) 𝑣 𝑡, 𝑥(𝑡) = (𝑣 𝑡 + 𝑣 · ∇𝑣) 𝑡, 𝑥(𝑡) .



¥ =
d𝑡
Furthermore, we assume to have an ideal fluid having a well-defined pressure 𝑝(𝑡, 𝑥) such
that if 𝑆 is any surface in the fluid with unit normal 𝑛, the force of the stress exerted
across the surface per unit area at time 𝑡 is 𝑝(𝑡, 𝑥)𝑛. Thus, the total force exerted on the
fluid in the volume 𝑉 by the pressure is given by
∫
− 𝑝𝑛.
𝜕𝑉

For any fixed vector 𝑒, we have

∫ ∫ ∫ ∫
−𝑒 · 𝑝𝑛 = − 𝑝𝑒 · 𝑛 = − div(𝑝𝑒) = − ∇𝑝 · 𝑒.
𝜕𝑉 𝜕𝑉 𝑉 𝑉

If 𝑏(𝑡, 𝑥) is an additional body force per unit mass such as the gravitational force
𝑏(𝑡, 𝑥) = −𝑔𝑒3 with 𝑔 being the gravitational acceleration and 𝑒3 = (0, 0, 1)𝑇 , the total
force exerted on any fluid particle is

−∇𝑝 + 𝜚𝑏.

By Newton’s second law, the total force on a fluid particle is equal to its mass times its
acceleration, i.e.,
𝜚(𝑣 𝑡 + 𝑣 · ∇𝑣) = −∇𝑝 − 𝜚 𝑔𝑒3 .

Incompressibility We assume that the fluid is incompressible. Following Chorin and

Marsden (1993), this can be translated into

∇ · 𝑣 = 0.

If you are interested in the steps to derive this equation, you can find them in the
book of Chorin and Marsden (1993) — or you can attend the course “Modellieren mit
partiellen Differentialgleichungen” by Ferdinand Thein in the summer term 2025.

Incompressible Euler equations of a homogeneous ideal fluid Summing up the re-

sults above, we obtain the Euler equations of an incompressible and homogeneous
fluid with constant density and velocity 𝑣 = (𝑣1 , 𝑣2 , 𝑣3 )𝑇 :

div 𝑣 = ∇ · 𝑣 = 0,
1
𝜕𝑡 𝑣 + (𝑣 · ∇)𝑣 + ∇𝑝 = −𝑔𝑒3 ,
𝜚

where 𝜚 is the constant density, 𝑝 is the pressure, and 𝑔 is the gravitational acceleration.
The velocity 𝑣(𝑡, 𝑥) depends on time 𝑡 and space 𝑥. The vector 𝑒3 = (0, 0, 1)𝑇 is the unit
vector in the vertical direction. The pressure is basically a Lagrange multiplier that
enforces the incompressibility condition and not a primary variable.

13
3.1.1 Irrotational flows
For the derivation of shallow water models, we will assume that the vorticity

𝜔 := curl 𝑣 = ∇ × 𝑣

vanishes. We can use this to simplify the momentum equation by using the vector
calculus identity
1
∇|𝑣|2 = (𝑣 · ∇)𝑣 + 𝑣 × (∇ × 𝑣),
2
which we can prove as follows using the Einstein summation convention. First, we
have  
1 1
∇|𝑣| = 𝜕𝑖 (𝑣 𝑗 𝑣 𝑗 ) = 𝑣 𝑗 𝜕𝑖 𝑣 𝑗 .
2
2 𝑖 2
Second, we have

(𝑣 · ∇)𝑣 + 𝑣 × (∇ × 𝑣) 𝑖 = 𝑣 𝑗 𝜕 𝑗 𝑣 𝑖 + 𝜀𝑖𝑗 𝑘 𝑣 𝑗 𝜀𝑘ℓ 𝑚 𝜕ℓ 𝑣 𝑚

= 𝑣 𝑗 𝜕 𝑗 𝑣 𝑖 + (𝛿 𝑖ℓ 𝛿 𝑗𝑚 − 𝛿 𝑖𝑚 𝛿 𝑗ℓ )𝑣 𝑗 𝜕ℓ 𝑣 𝑚
= 𝑣 𝑗 𝜕 𝑗 𝑣 𝑖 + 𝑣 𝑗 𝜕𝑖 𝑣 𝑗 − 𝑣 𝑗 𝜕 𝑗 𝑣 𝑖 = 𝑣 𝑗 𝜕𝑖 𝑣 𝑗 ,

where we used the well-known properties

𝜀𝑖𝑗 𝑘 = 𝜀𝑘𝑖𝑗 , 𝜀𝑘𝑖𝑗 𝜀𝑘ℓ 𝑚 = 𝛿 𝑖ℓ 𝛿 𝑗𝑚 − 𝛿 𝑖𝑚 𝛿 𝑗ℓ

of the Levi-Civita symbol 𝜀𝑖𝑗 𝑘 and the Kronecker delta 𝛿 𝑖𝑗 . Thus, the momentum
equation of the incompressible Euler equations becomes
1 1
𝜕𝑡 𝑣 + ∇|𝑣|2 + ∇𝑝 = −𝑔𝑒3 .
2 𝜚

3.1.2 Boundary conditions

Since we are interested in shallow water models, we will assume that the fluid is
bounded from below by the impermeable bottom topography. Thus, we have the
boundary condition
𝑣 · 𝑛 = 0 on the bottom boundary,
where 𝑛(𝑡, 𝑥) is the outward unit normal vector on the boundary.
At the top boundary, the situation is more complicated since we have a free surface.
Assume that this free surface is given as

{𝑥 ∈ R𝑑 | 𝐹(𝑡, 𝑥) = 0}.

Then, the unit normal vector on the free surface is given by

∇𝐹(𝑡, 𝑥)
𝑛(𝑡, 𝑥) = .
|∇𝐹(𝑡, 𝑥)|
Next, we consider a point 𝑥(𝑡) on the free surface. The velocity 𝑥(𝑡,¤ 𝑥(𝑡)) of this point is
given by
d
0= 𝐹(𝑡, 𝑥(𝑡)) = 𝜕𝑡 𝐹(𝑡, 𝑥(𝑡)) + ∇𝐹(𝑡, 𝑥(𝑡)) · 𝑥(𝑡,
¤ 𝑥(𝑡)),
d𝑡
i.e.,
𝜕𝑡 𝐹(𝑡, 𝑥(𝑡))
𝑥(𝑡,
¤ 𝑥(𝑡)) · 𝑛(𝑡, 𝑥(𝑡)) = − .
|∇𝐹(𝑡, 𝑥(𝑡))|

14
The kinematic boundary condition at the free surface states that this normal velocity of
the surface is equal to the normal velocity of the fluid itself, i.e.,

𝑥(𝑡,
¤ 𝑥(𝑡)) · 𝑛(𝑡, 𝑥(𝑡)) = 𝑣(𝑡, 𝑥(𝑡)) · 𝑛(𝑡, 𝑥(𝑡)).

Written in terms of the level-set function 𝐹, this condition reads

𝜕𝑡 𝐹 + 𝑣 · ∇𝐹 = 0.

When the free surface can be parameterized by a function 𝜂(𝑡, 𝑥), i.e.,

𝐹(𝑡, 𝑥) = 𝜂(𝑡, 𝑥1 , . . . , 𝑥 𝑑−1 ) − 𝑥 𝑑 ,

this becomes
𝑑−1
Õ
𝜕𝑡 𝜂 + 𝑣 𝑖 𝜕𝑥 𝑖 𝜂 = 𝑣 𝑑 on the free surface.
𝑖=1

Moreover, the dynamic boundary condition states that the stresses on the free surface
are equal. Neglecting surface tension and the motion of the much less dense air,
this condition states that the pressure on the free surface is equal to the atmospheric
pressure 𝑝 𝑎 :
𝑝 = 𝑝 𝑎 on the free surface.
Since the pressure acts as a Lagrange multiplier and is only defined up to a constant,
we can set 𝑝 𝑎 = 0 without loss of generality.
→ end of lecture 2025-04-17

3.1.3 Non-dimensionalization
From now on, we simplify the following derivations by focusing on the 2D case with
flat bathymetry at 𝑥 2 = −ℎ 0 , where ℎ 0 is the undisturbed water height. Thus, the
undisturbed water surface is at 𝜂 = 0. In this case, the incompressible Euler equations
with curl-free velocity become

𝑣 1,𝑥 + 𝑣 2,𝑦 = 0, (∇ · 𝑣 = 0)
𝑣 1,𝑦 = 𝑣2,𝑥 , (curl 𝑣 = 0)
1
𝑣 1,𝑡 + 𝑣 1 𝑣1,𝑥 + 𝑣 2 𝑣 1,𝑦 + 𝑝 𝑥 = 0,
𝜚
1
𝑣2,𝑡 + 𝑣 1 𝑣2,𝑥 + 𝑣 2 𝑣 2,𝑦 + 𝑝 𝑦 = −𝑔,
𝜚

where 𝑣1 is the horizontal velocity, 𝑣2 is the vertical velocity, and 𝑝 𝑥 /𝑝 𝑦 are the hori-
zontal/vertical components of the pressure gradient.
We will now non-dimensionalize the equations by introducing the following charac-
teristic scales:

• characteristic water height ℎ0

• wave amplitude 𝐴

• wavelength 𝜆

15
The classical shallow water assumption is the undisturbed water depth is much smaller
than the typical wavelengths, i.e.,

ℎ0
𝜎 := ≪ 1.
𝜆
We will see later that 𝜎 determines the strength of the dispersive effects. Moreover, we
sometimes assume only weakly nonlinear waves, i.e.,

𝐴
𝜀 := ≪ 1.
ℎ0
Then, we introduce the non-dimensionalized variables
𝜎𝑔 1 1
𝑡∗ = 𝑡, 𝑥∗ = 𝑥, 𝑦∗ = 𝑦,
𝑐0 𝜆 ℎ0

where 𝑐 0 = 𝑔 ℎ 0 is the characteristic wave speed. Furthermore, we non-dimensionalize

p
the velocity and pressure by

1 1 1
𝑣1∗ = 𝑣1 , 𝑣2∗ = 𝑣2 , 𝑝∗ = 𝑝.
𝜀𝑐 0 𝜀𝜎𝑐0 𝜚 𝑔 ℎ0

By introducing this scaling, we assume that the non-dimensionalized variables are

of order unity and that relative strengths of the different terms is determined by the
(small) parameters 𝜀 and 𝜎.
The next step is to substitute the non-dimensionalized variables into the incompress-
ible Euler equations. Thus, we have
∗ ∗
𝑣1,𝑡 = 𝜀𝑐0 𝑣1,𝑡 = 𝜀𝜎𝑔𝑣1,𝑡 ∗,

∗ 𝜀𝑐 0 ∗
𝑣1,𝑥 = 𝜀𝑐0 𝑣1,𝑥 = 𝑣 ∗,
𝜆 1,𝑥
∗ 𝜀𝑐 0 ∗
𝑣1,𝑦 = 𝜀𝑐0 𝑣1,𝑦 = 𝑣 ∗,
ℎ0 1,𝑦
∗ ∗
𝑣2,𝑡 = 𝜀𝜎𝑐0 𝑣2,𝑡 = 𝜀𝜎2 𝑔𝑣2,𝑡 ∗,

∗ 𝜀𝜎𝑐0 ∗
𝑣2,𝑥 = 𝜀𝜎𝑐0 𝑣2,𝑥 = 𝑣 ∗,
𝜆 2,𝑥
∗ 𝜀𝜎𝑐0 ∗
𝑣2,𝑦 = 𝜀𝜎𝑐0 𝑣2,𝑦 = 𝑣 ∗,
ℎ0 2,𝑦
𝜚 𝑔 ℎ0 ∗
𝑝 𝑥 = 𝜚 𝑔 ℎ0 𝑝 ∗𝑥 = 𝑝 ∗,
𝜆 𝑥
𝑝 𝑦 = 𝜚 𝑔 ℎ0 𝑝 ∗𝑦 = 𝜚 𝑔𝑝 ∗𝑦 ∗ .

Substituting these relations into the incompressible Euler equations, we obtain

𝜀𝑐 0 ∗ 𝜀𝜎𝑐0 ∗
𝑣1,𝑥 ∗ + 𝑣 ∗ = 0, (∇ · 𝑣 = 0)
𝜆 ℎ0 2,𝑦
𝜀𝑐 0 ∗ 𝜀𝜎𝑐0 ∗
𝑣1,𝑦 ∗ = 𝑣 ∗, (curl 𝑣 = 0)
ℎ0 𝜆 2,𝑥
∗ ∗ 𝜀𝑐 0 ∗ 𝜀𝑐0 ∗ 1 𝜚 𝑔 ℎ0 ∗
𝜀𝜎 𝑔𝑣 1,𝑡 ∗ + 𝜀𝑐 0 𝑣 1 𝑣1,𝑥 ∗ + 𝜀𝜎𝑐 0 𝑣2∗ 𝑣1,𝑦 ∗ + 𝑝 ∗ = 0,
𝜆 ℎ0 𝜚 𝜆 𝑥
∗ ∗ 𝜀𝜎𝑐 0 ∗ 𝜀𝜎𝑐0 ∗ 1 𝜚 𝑔 ℎ0 ∗
𝜀𝜎 2 𝑔𝑣2,𝑡 ∗ + 𝜀𝑐 0 𝑣 1 𝑣2,𝑥 ∗ + 𝜀𝜎𝑐 0 𝑣2∗ 𝑣2,𝑦 ∗ + 𝑝 ∗ = −𝑔,
𝜆 ℎ0 𝜚 ℎ0 𝑦

16
which we can reformulate using 𝜎 = ℎ0 /𝜆 and 𝑐 02 = 𝑔 ℎ0 as
∗ ∗
𝑣 1,𝑥 ∗ + 𝑣 2,𝑦 ∗ = 0,

∗ 2 ∗
𝑣 1,𝑦 ∗ = 𝜎 𝑣 2,𝑥 ∗ ,

∗ 2 ∗ ∗ 2 ∗ ∗ ∗ (3.1)
𝜀𝑣1,𝑡 ∗ + 𝜀 𝑣 1 𝑣 1,𝑥 ∗ + 𝜀 𝑣 2 𝑣 1,𝑦 ∗ + 𝑝 𝑥 ∗ = 0,

∗ 2 2 ∗ ∗ 2 2 ∗ ∗ ∗
𝜀𝜎 2 𝑣2,𝑡 ∗ + 𝜀 𝜎 𝑣 1 𝑣 2,𝑥 ∗ + 𝜀 𝜎 𝑣 2 𝑣 2,𝑦 ∗ + 𝑝 𝑦 ∗ = −1.

Finally, we non-dimensionalize the free surface elevation by

1
𝜂∗ = 𝜂.
𝜀ℎ 0
Then, the boundary conditions become

𝑣2∗ = 0 for 𝑦 ∗ = −1,

𝜂𝑡∗∗ + 𝜀𝑣 1∗ 𝜂∗𝑥 ∗ = 𝑣2∗ for 𝑦 ∗ = 𝜀𝜂∗ (𝑡 ∗ , 𝑥 ∗ ), (3.2)
𝑝∗ = 0 for 𝑦 ∗ = 𝜀𝜂∗ (𝑡 ∗ , 𝑥 ∗ ).

3.2 Deriving simplified models

From here on, we will omit the asterisks for the non-dimensionalized variables and drop
the superscript ∗ in (3.1) and (3.2). The free-boundary problem for the incompressible
Euler equations is difficult to solve. Thus, people prefer solving simplified models that
are derived by averaging over the depth. Thus, we introduce the vertically averaged
velocity ∫ 𝜀𝜂
1
𝑣 := 𝑣 1 d𝑦
1 + 𝜀𝜂 −1
and the total water height
ℎ := 1 + 𝜀𝜂.
Integrating the divergence constraint of the incompressible Euler equations (3.1) over
𝑦 yields ∫ 𝜀𝜂
𝑣 1,𝑥 d𝑦 + 𝑣 2 (𝑡, 𝑥, 𝜀𝜂) − 𝑣2 (𝑡, 𝑥, −1) = 0.
−1
By Leibniz’ rule
∫ 𝜀𝜂 ∫ 𝜀𝜂
𝜕𝑥 𝑣 1 (𝑡, 𝑥, 𝑦) d𝑦 = 𝑣1,𝑥 (𝑡, 𝑥, 𝑦) d𝑦 + 𝑣 1 (𝑡, 𝑥, 𝜀𝜂)𝜕𝑥 𝜀𝜂(𝑡, 𝑥)


−1 −1
| {z }
=𝜕𝑥 (ℎ𝑣)

for differentiating integrals with parameters, we obtain

0 = 𝜕𝑥 (ℎ𝑣)(𝑡, 𝑥) −𝑣 1 (𝑡, 𝑥, 𝜀𝜂)𝜕𝑥 𝜀𝜂(𝑡, 𝑥) +𝑣 2 (𝑡, 𝑥, 𝜀𝜂) − 𝑣2 (𝑡, 𝑥, −1) .

| {z } | {z }
=𝜂𝑡 −𝑣2 =0

Inserting the kinematic boundary condition 𝜂𝑡 + 𝜀𝑣 1 𝜂 𝑥 = 𝑣2 at the free surface and

the impermeability boundary condition 𝑣 2 = 0 at the bottom, we obtain the continuity
equation
𝜂𝑡 + (ℎ𝑣)𝑥 = 0.

17
 Next, we integrate the momentum equation over 𝑦 to obtain
∫ 𝜀𝜂 ∫ 𝜀𝜂 ∫ 𝜀𝜂 ∫ 𝜀𝜂
𝜀 𝑣 1,𝑡 d𝑦 + 𝜀 2
𝑣1 𝑣1,𝑥 d𝑦 + 𝜀 2
𝑣 2 𝑣1,𝑦 d𝑦 + 𝑝 𝑥 d𝑦 = 0.
−1 −1 −1 −1

By Leibniz’ rule, we have

∫ 𝜀𝜂 ∫ 𝜀𝜂
𝑣 1,𝑡 d𝑦 = 𝜕𝑡 𝑣1 d𝑦 −𝑣 1 (𝑡, 𝑥, 𝜀𝜂)𝜕𝑡 𝜀𝜂(𝑡, 𝑥)


−1 −1
| {z }
=ℎ𝑣

and ∫ 𝜀𝜂 ∫ 𝜀𝜂
1 1
𝑣1 𝑣1,𝑥 d𝑦 = 𝜕𝑥 𝑣12 d𝑦 − 𝑣12 (𝑡, 𝑥, 𝜀𝜂)𝜕𝑥 𝜀𝜂(𝑡, 𝑥) .


−1 2 −1 2
The third integral can be simplified using integration by parts and the divergenc con-
straint (i.e., 𝑣2,𝑦 = −𝑣1,𝑥 ) as
∫ 𝜀𝜂 ∫ 𝜀𝜂
𝑣2 𝑣1,𝑦 d𝑦 = 𝑣2 (𝑡, 𝑥, 𝜀𝜂)𝑣 1 (𝑡, 𝑥, 𝜀𝜂) − 𝑣2 (𝑡, 𝑥, −1)𝑣 1 (𝑡, 𝑥, −1) − 𝑣1 𝑣 2,𝑦 d𝑦.
−1 −1 |{z}
=−𝑣1 𝑣1,𝑥

Thus, we obtain
=ℎ 𝑡 𝑣+ℎ𝑣 𝑡 =𝜀𝜂𝑡 𝑣+ℎ𝑣 𝑡 =−(ℎ𝑣)𝑥 𝑣+ℎ𝑣 𝑡
z}|{ ∫ 𝜀𝜂
0 = 𝜀 (ℎ𝑣)𝑡 −𝜀𝑣 1 (𝑡, 𝑥, 𝜀𝜂)𝜕𝑡 𝜀𝜂(𝑡, 𝑥) + 𝜀 𝜕𝑥 𝑣12 d𝑦 − 𝜀2 𝑣12 (𝑡, 𝑥, 𝜀𝜂)𝜕𝑥 𝜀𝜂(𝑡, 𝑥)

2


−1
∫ 𝜀𝜂
+ 𝜀 𝑣2 (𝑡, 𝑥, 𝜀𝜂)𝑣 1 (𝑡, 𝑥, 𝜀𝜂) − 𝜀 𝑣2 (𝑡, 𝑥, −1)𝑣 1 (𝑡, 𝑥, −1) +
2 2
𝑝 𝑥 d𝑦
−1
 
= 𝜀ℎ𝑣 𝑡 − 𝜀(ℎ𝑣)𝑥 𝑣 − 𝜀2 𝑣1 (𝑡, 𝑥, 𝜀𝜂) 𝜕𝑡 𝜂(𝑡, 𝑥) + 𝑣 1 (𝑡, 𝑥, 𝜀𝜂)𝜕𝑥 𝜀𝜂(𝑡, 𝑥) − 𝑣 2 (𝑡, 𝑥, 𝜀𝜂)


∫ 𝜀𝜂 ∫ 𝜀𝜂
+ 𝜀 𝜕𝑥
2
𝑣12 d𝑦 − 𝜀 𝑣2 (𝑡, 𝑥, −1)𝑣 1 (𝑡, 𝑥, −1) +
2
𝑝 𝑥 d𝑦.
−1 −1

Next, we use again the kinematic boundary condition 𝜂𝑡 + 𝜀𝑣1 𝜂 𝑥 = 𝑣2 at the free surface
and the impermeability boundary condition 𝑣2 = 0 at the bottom to get rid of the
boundary terms. Next, we compute
∫ 𝜀𝜂
−𝜀(ℎ𝑣)𝑥 𝑣 = −𝜀(ℎ𝑣 )𝑥 + 𝜀(ℎ𝑣)𝑣 𝑥 = −𝜀𝜕𝑥
2
𝑣 2 d𝑦 + 𝜀(ℎ𝑣)𝑣 𝑥 .
−1

Finally, we use Leibniz rule again for the pressure term

∫ 𝜀𝜂 ∫ 𝜀𝜂
𝑝 𝑥 d𝑦 = 𝜕𝑥 𝑝 d𝑦 − 𝑝(𝑡, 𝑥, 𝜀𝜂) 𝜕𝑥 𝜀𝜂(𝑡, 𝑥) =: 𝜕𝑥 (ℎ𝑝).


−1 −1 | {z }
=0

To sum up, we have obtaind the (exact!) system

ℎ 𝑡 + 𝜀(ℎ𝑣)𝑥 = 0,
∫ 𝜀𝜂   (3.3)
𝜀ℎ𝑣 𝑡 + 𝜀 ℎ𝑣 𝑣 𝑥 + 𝜀 𝜕𝑥
2 2
𝑣12 −𝑣 2
d𝑦 + 𝜕𝑥 (ℎ𝑝) = 0.
−1

Next, we will derive simplified shallow water models from this formulation by approx-
imating the remaining integrals and ignoring higher-order terms in 𝜀 and 𝜎.
→ end of lecture 2025-04-22

18
3.2.1 Handling the pressure term
First, we rewrite the second momentum equation of the incompressible Euler equations
(3.1) as
𝑝 𝑦 (𝑡, 𝑥, 𝑦) = −1 − 𝜀𝜎 2 Γ(𝑡, 𝑥, 𝑦),
where we have defined
Γ = 𝑣2,𝑡 + 𝜀𝑣 1 𝑣2,𝑥 + 𝜀𝑣 2 𝑣2,𝑦 .
Thus, we have
=0
∫ 𝜀𝜂 z }| { ∫ 𝜀𝜂
𝑝(𝑡, 𝑥, 𝑦) = − 𝑝 𝑦 (𝑡, 𝑥, 𝑠) d𝑠 − 𝑝(𝑡, 𝑥, 𝜀𝜂) = (𝜀𝜂 − 𝑦) + 𝜀𝜎 2
Γ(𝑡, 𝑥, 𝑠) d𝑠
𝑦 𝑦

due to the dynamic boundary condition 𝑝(𝑡, 𝑥, 𝜀𝜂) = 0 at the free surface. Thus, we
obtain the mean value
!
∫ 𝜀𝜂 ∫ 𝜀𝜂 ∫ 𝜀𝜂
ℎ𝑝 = 𝑝(𝑡, 𝑥, 𝑦) d𝑦 = 𝜀𝜂 − 𝑦 + 𝜀𝜎2 Γ(𝑡, 𝑥, 𝑠) d𝑠 d𝑦
−1 −1 𝑦
  𝜀𝜂 ∫ 𝜀𝜂 ∫ 𝜀𝜂
1
= 𝜀𝜂𝑦 − 𝑦 2 +𝜀𝜎 2
Γ(𝑡, 𝑥, 𝑠) d𝑠 d𝑦.
2 −1 −1 𝑦
| {z }
= 12 𝜀𝜂2 +𝜀𝜂+ 21 = 12 (1+𝜀𝜂)2 = 21 ℎ 2

Thus, we have
∫ 𝜀𝜂 ∫ 𝜀𝜂
1
𝜕𝑥 (ℎ𝑝) = 𝜕𝑥 (ℎ 2 ) + 𝜀𝜎 2 𝜕𝑥 Γ(𝑡, 𝑥, 𝑠) d𝑠 d𝑦
2 −1 𝑦

with
1
𝜕𝑥 (ℎ 2 ) = ℎ ℎ 𝑥 = 𝜀ℎ𝜂 𝑥 .
2
Hence, the second equation of (3.3) becomes
∫ 𝜀𝜂   ∫ 𝜀𝜂 ∫ 𝜀𝜂
𝜀ℎ𝑣 𝑡 + 𝜀ℎ𝜂 𝑥 + 𝜀 ℎ𝑣 𝑣 𝑥 + 𝜀 𝜕𝑥
2 2
𝑣12 −𝑣 2
d𝑦 + 𝜀𝜎 𝜕𝑥 2
Γ(𝑡, 𝑥, 𝑠) d𝑠 d𝑦 = 0.
−1 −1 𝑦

3.2.2 Handling the integrals

Next, we approximate the remaining integrals.

Velocity integralFirst, we approximate the horizontal velocity 𝑣 1 by its Taylor poly-

nomial around the bottom at 𝑦 = 1, i.e.,

1
𝑣1 (𝑡, 𝑥, 𝑦) = 𝑣1 (𝑡, 𝑥, −1) + (𝑦 + 1) 𝑣1,𝑦 (𝑡, 𝑥, −1) + (𝑦 + 1)2 𝑣1,𝑦 𝑦 (𝑡, 𝑥, −1) + . . .
| {z } 2
=𝜎2 𝑣2,𝑥 (𝑡,𝑥,−1)

Next, we use the irrotational condition 𝑣1,𝑦 = 𝜎2 𝑣2,𝑥 and the incompressibility condition
𝑣1,𝑥 + 𝑣 2,𝑦 = 0 from (3.1) to rewrite

𝑣 1,𝑦 𝑦 = 𝜎2 𝑣 2,𝑥 𝑦 = 𝜎2 𝑣2,𝑦𝑥 = −𝜎2 𝑣1,𝑥𝑥 .

19
Hence, we also see that higher-order derivatives like 𝑣1,𝑦 𝑦 𝑦 and 𝑣 1,𝑦 𝑦 𝑦 𝑦 are 𝑂(𝜎 4 ).
Moreover, we use the boundary condition 𝑣2 ≡ 0 at the bottom boundary to deduce
𝑣2,𝑥 (𝑡, 𝑥, −1) = 0. Thus, we obtain
1
𝑣1 (𝑡, 𝑥, 𝑦) = 𝑣1 (𝑡, 𝑥, −1) − 𝜎 2 (𝑦 + 1)2 𝑣1,𝑥𝑥 (𝑡, 𝑥, −1) + 𝑂(𝜎4 ). (3.4)
2
Integrating yields
∫ 𝜀𝜂
ℎ𝑣(𝑡, 𝑥) = 𝑣 1 (𝑡, 𝑥, 𝑦) d𝑦
−1
∫ 𝜀𝜂 ∫ 𝜀𝜂
1
= 𝑣1 (𝑡, 𝑥, −1) d𝑦 − 𝜎2 𝑣1,𝑥𝑥 (𝑡, 𝑥, −1) (𝑦 + 1)2 d𝑦 + 𝑂(𝜎 4 )
−1 2 −1
1
= ℎ𝑣1 (𝑡, 𝑥, −1) − 𝜎 2 ℎ 3 𝑣1,𝑥𝑥 (𝑡, 𝑥, −1) + 𝑂(𝜎4 ).
6
Canceling ℎ, we obtain

𝑣 𝑥𝑥 (𝑡, 𝑥) = 𝑣 1,𝑥𝑥 (𝑡, 𝑥, −1) + 𝒪(𝜎2 ).

Thus, we get
1
𝑣1 (𝑡, 𝑥, −1) = 𝑣(𝑡, 𝑥) + 𝜎 2 ℎ 2 𝑣 𝑥𝑥 (𝑡, 𝑥) + 𝑂(𝜎 4 ).
6
Hence, we obtain
1 1
𝑣 1 (𝑡, 𝑥, 𝑦) = 𝑣(𝑡, 𝑥) + 𝜎 2 ℎ 2 𝑣 𝑥𝑥 (𝑡, 𝑥) − 𝜎2 (𝑦 + 1)2 𝑣 𝑥𝑥 (𝑡, 𝑥) + 𝑂(𝜎 4 ). (3.5)
6 2
Squaring this equation results in
1
𝑣12 = 𝑣 2 + 𝜎2 ℎ 2 𝑣 𝑣 𝑥𝑥 − 𝜎 2 (𝑦 + 1)2 𝑣 𝑣 𝑥𝑥 + 𝑂(𝜎 4 ).
3
Integrating this and using
∫ 𝜀𝜂  
 𝜀𝜂 1  2 
1 2 1 2
ℎ − (𝑦 + 1) d𝑦 =
2
ℎ 𝑦 − (𝑦 + 1)3 −1 = ℎ 𝜀𝜂 + ℎ 2 − ℎ 3 = 0
−1 3 3 3
yields ∫ 𝜀𝜂
(𝑣12 − 𝑣 2 ) d𝑦 = 𝑂(𝜎 4 ).
−1
These terms are also referred to as “turbulent terms” since they are associated with
fine-scale structures of the velocity deviating from the mean velocity.

Pressure integral Next, we analyze the pressure term involving the integral of

Γ = 𝑣2,𝑡 + 𝜀𝑣 1 𝑣2,𝑥 + 𝜀𝑣 2 𝑣2,𝑦 .

We compute the Taylor expansion of the vertical velocity 𝑣2 around the bottom:

=0 =−𝑣 1,𝑥 (𝑡,𝑥,−1) =−𝜎2 𝑣2,𝑥𝑥 (𝑡,𝑥,−1)=0

z }| { z }| { z }| {
1
𝑣2 (𝑡, 𝑥, 𝑦) = 𝑣2 (𝑡, 𝑥, −1) +(𝑦 + 1) 𝑣2,𝑦 (𝑡, 𝑥, −1) + (𝑦 + 1)2 𝑣 2,𝑦 𝑦 (𝑡, 𝑥, −1)
2
1
+ (𝑦 + 1)3 𝑣2,𝑦 𝑦 𝑦 (𝑡, 𝑥, −1) +𝑂(𝜎 4 ),
6 | {z }
=𝜎2 𝑣 1,𝑥𝑥𝑥 (𝑡,𝑥,−1)

20
where we have used the boundary condition, the irrotational condition 𝑣 2,𝑦 = −𝑣1,𝑥 ,
and the incompressibility condition 𝜎2 𝑣 2,𝑥 = 𝑣1,𝑦 . Inserting

𝑣2 (𝑡, 𝑥, 𝑦) = −(𝑦 + 1)𝑣 1,𝑥 (𝑡, 𝑥, −1) + 𝑂(𝜎2 ) = −(𝑦 + 1)𝑣 𝑥 (𝑡, 𝑥) + 𝑂(𝜎 2 )

and
𝑣 1 (𝑡, 𝑥, 𝑦) = 𝑣(𝑡, 𝑥) + 𝑂(𝜎 2 )
into the equation defining Γ yields
 
Γ = −(𝑦 + 1) 𝑣 𝑡𝑥 + 𝜀𝑣 𝑣 𝑥𝑥 − 𝜀(𝑣 𝑥 )2 + 𝒪(𝜎 2 ).

Thus, we have
∫ 𝜀𝜂 ∫ 𝜀𝜂  ∫ 𝜀𝜂 ∫ 𝜀𝜂
Γ(𝑡, 𝑥, 𝑠) d𝑠 d𝑦 = 𝑣 𝑡𝑥 + 𝜀𝑣 𝑣 𝑥𝑥 − 𝜀(𝑣 𝑥 ) 2
−(𝑠 + 1) d𝑠 d𝑦 + 𝒪(𝜎 2 ),
−1 𝑦 −1 𝑦

where we compute the integral as

∫ 𝜀𝜂 ∫ 𝜀𝜂 ∫ 𝜀𝜂   𝜀𝜂
1 2
− (𝑠 + 1) d𝑠 d𝑦 = − 𝑠 +𝑠 d𝑦
−1 𝑦 −1 2 𝑦
∫ 𝜀𝜂  
1 2 2 1
=− 𝜀 𝜂 + 𝜀𝜂 − 𝑦 2 − 𝑦 d𝑦
−1 2 2
  𝜀𝜂
1 1 1
= − 𝜀2 𝜂2 𝑦 + 𝜀𝜂𝑦 − 𝑦 3 − 𝑦 2
2 6 2 −1
1 3 3 1 3 3 1 2 2 1 2 2 1 1
= − 𝜀 𝜂 − 𝜀 𝜂 + 𝜀 𝜂 + 𝜀 𝜂 − 𝜀 𝜂 − 𝜀𝜂 − +
2 2
2 6 2 2 6 2
1 3 3 1 1 1
= − 𝜀 𝜂 − 𝜀2 𝜂2 − 𝜀𝜂 − = − (1 + 𝜀𝜂)3 = − ℎ 3 .
3 3 3 3
Hence, we obtain
∫ 𝜀𝜂 ∫ 𝜀𝜂
1 3 
Γ(𝑡, 𝑥, 𝑠) d𝑠 d𝑦 = − ℎ 𝑣 𝑡𝑥 + 𝜀𝑣 𝑣 𝑥𝑥 − 𝜀(𝑣 𝑥 ) + 𝑂(𝜎 2 ).
2
−1 𝑦 3

These terms are also referred to as “non-hydrostatic terms”, since they are associated
to the deviation of the pressure from hydrostatic equilibrium.

Combining the results Combining the results, we obtain the system

𝜕𝑡 ℎ + 𝜀(ℎ𝑣)𝑥 = 0,
   
1 2 (3.6)
𝑣 𝑡 + 𝜂 𝑥 + 𝜀𝑣 𝑣 𝑥 − 𝜎 𝜕𝑥 ℎ 3 𝑣 𝑡𝑥 + 𝜀𝑣 𝑣 𝑥𝑥 − 𝜀(𝑣 𝑥 )2 = 𝑂(𝜀𝜎4 ).
3ℎ

3.3 The Serre-Green-Naghdi (SGN) equations

The Serre-Green-Naghdi (SGN) equations are derived for the nonlinear (𝜀 = 𝑂(1)) but
weakly dispersive regime, i.e., higher-order terms 𝑂(𝜎4 ) are neglected in (3.6), resulting
in the non-dimensionalized system

ℎ 𝑡∗∗ + 𝜀(ℎ ∗ 𝑣 ∗ )𝑥 ∗ = 0,
   
1
𝑣 ∗𝑡 ∗ + 𝜂∗𝑥 ∗ + 𝜀𝑣 ∗
𝑣 ∗𝑥 ∗ − ∗ 𝜎 2 𝜕𝑥 ∗ (ℎ ∗ )3 𝑣 ∗𝑡 ∗ 𝑥 ∗ + 𝜀𝑣 ∗ 𝑣 ∗𝑥 ∗ 𝑥 ∗ − 𝜀(𝑣 ∗𝑥 ∗ )2 = 0,
3ℎ

21
where we have re-introduced the notation ·∗ for non-dimensionalized variables. Recall
the nondimensionalization
𝜎𝑔 1 1 1 1 1
𝑡∗ = 𝑡, 𝑥∗ = 𝑥, 𝜂∗ = 𝜂, ℎ ∗ = 1 + 𝜀𝜂∗ = (ℎ0 + 𝜂) = ℎ, 𝑣∗ = 𝑣,
𝑐0 𝜆 𝜀ℎ 0 ℎ0 ℎ0 𝜀𝑐0

with 𝑐 0 = 𝑔 ℎ 0 and 𝜎 = ℎ0 /𝜆. Turning back to the dimensional form, we transform

p
the terms as
1 𝑐0
ℎ 𝑡∗∗ = ℎ𝑡 ,
ℎ0 𝜎𝑔
1 1
𝜀(ℎ ∗ 𝑣 ∗ )𝑥 ∗ = 𝜀 𝜆(ℎ𝑣)𝑥 ,
ℎ0 𝜀𝑐 0
1 𝑐0
𝑣 ∗𝑡 ∗ = 𝑣𝑡 ,
𝜀𝑐0 𝜎𝑔
1
𝜂∗𝑥 ∗ = 𝜆𝜂 𝑥 ,
𝜀ℎ 0
1
𝜀𝑣 ∗ 𝑣 ∗𝑥 ∗ = 𝜀 𝜆𝑣 𝑣 𝑥 ,
𝜀2 𝑐02
and
1 𝜎 2 𝜆3
       
1 2 ∗ 3 ∗ ∗ ∗ ∗ 2
𝜎 𝜕𝑥 (ℎ ) 𝑣 𝑡 ∗ 𝑥 ∗ + 𝜀𝑣 𝑣 𝑥 ∗ 𝑥 ∗ − 𝜀(𝑣 𝑥 ∗ )
∗ = 𝜕𝑥 ℎ 𝑣 𝑡𝑥 + 𝜀𝑣 𝑣 𝑥𝑥 − 𝜀(𝑣 𝑥 )
3 2
.
3ℎ ∗ 3ℎ 𝜀𝑐02 ℎ 02

Thus, we get the classical SGN equations

𝜕𝑡 ℎ + (ℎ𝑣)𝑥 = 0,
   
1 (3.7)
𝑣𝑡 + 𝑔 ℎ𝑥 + 𝑣 𝑣 𝑥 − 𝜕𝑥 ℎ 𝑣 𝑡𝑥 + 𝑣 𝑣 𝑥𝑥 − (𝑣 𝑥 )
3 2
= 0,
3ℎ

where ℎ = ℎ0 + 𝜂 is the total water height above the flat bathymetry at 𝑦 = −ℎ0 .

3.4 Boussinesq approximations and the BBM-BBM system

In contrast to the SGN equations, Boussinesq approximations also required small am-
plitudes, i.e., not only 𝜎 ≪ 1 but also 𝜀 ≪ 1. Writing (3.6) in terms of 1 + 𝜀𝜂 instead of
ℎ, we obtain

𝜂𝑡 + (1 + 𝜀𝜂)𝑣


𝑥 = 0,
  
1
𝑣 𝑡 + 𝜂 𝑥 + 𝜀𝑣 𝑣 𝑥 − 𝜎 2 𝜕𝑥 (1 + 𝜀𝜂)3 𝑣 𝑡𝑥 + 𝜀𝑣 𝑣 𝑥𝑥 − 𝜀(𝑣 𝑥 )2 = 𝑂(𝜀𝜎4 ).
3 + 3𝜀𝜂

Thus, we have

𝜂𝑡 + (1 + 𝜀𝜂)𝑣


= 0, 𝑥
1 (3.8)
𝑣 𝑡 + 𝜂 𝑥 + 𝜀𝑣 𝑣 𝑥 − 𝜎 2 𝑣 𝑡𝑥𝑥 = 𝑂(𝜀2 ) + 𝑂(𝜀𝜎 2 ) + 𝑂(𝜎 4 ).
3
Neglecting terms of order 𝑂(𝜀2 ), 𝑂(𝜀𝜎2 ), and 𝑂(𝜎 4 ) yields the so-called classical Boussi-
nesq system.
Next, we derive the BBM-BMM system from the classical Boussinesq system by
choosing another velocity variable following Bona, Chen, and Saut (2002). “Boussinesq
equations and other systems for small-amplitude long waves in nonlinear dispersive

22
media. I: Derivation and linear theory.” In: Journal of Nonlinear Science 12, pp. 283–318.
doi: 10.1007/s00332-002-0466-4. So far, we have used the averaged velocity
∫ 𝜀𝜂
1
𝑣= 𝑣 1 (𝑡, 𝑥, 𝑦) d𝑦.
ℎ −1

For the BBM system (and the general 𝑎𝑏𝑐𝑑-Boussinesq system), we use the velocity

𝑣 (𝑡, 𝑥) := 𝑣1 (𝑡, 𝑥, −1 + 𝜃).

e

Recall the Taylor expansion (3.4)

1
𝑣1 (𝑡, 𝑥, 𝑦) = 𝑣1 (𝑡, 𝑥, −1) − 𝜎 2 (𝑦 + 1)2 𝑣1,𝑥𝑥 (𝑡, 𝑥, −1) + 𝑂(𝜎4 ).
2
Thus, we have
1
𝑣 (𝑡, 𝑥) = 𝑣 1 (𝑡, 𝑥, −1) − 𝜎 2 𝜃2 𝑣1,𝑥𝑥 (𝑡, 𝑥, −1) + 𝑂(𝜎4 ).
e
2
By using the Taylor expansion of the inverse of the operator
  −1
1 1
I − 𝜎 2 𝜃2 𝜕𝑥2 = I + 𝜎2 𝜃2 𝜕𝑥2 + 𝑂(𝜎 4 )
2 2

(which can be made rigorous in Fourier space), we obtain

1
𝑣 1 (𝑡, 𝑥, −1) = e
𝑣 (𝑡, 𝑥) + 𝜎 2 𝜃2e
𝑣 𝑥𝑥 (𝑡, 𝑥) + 𝑂(𝜎 4 ).
2
Inserting this into the Taylor expansion (3.4), we obtain
1
𝑣 1 (𝑡, 𝑥, 𝑦) = 𝑣1 (𝑡, 𝑥, −1) − 𝜎2 (𝑦 + 1)2 𝑣 1,𝑥𝑥 (𝑡, 𝑥, −1) + 𝑂(𝜎4 )
2
1 2 2 1
=e𝑣 (𝑡, 𝑥) + 𝜎 𝜃 e 𝑣 𝑥𝑥 (𝑡, 𝑥) − 𝜎2 (𝑦 + 1)2e 𝑣 𝑥𝑥 (𝑡, 𝑥) + 𝑂(𝜎 4 ).
2 2
Thus, we can write the averaged velocity as
𝜀𝜂
∫
1
𝑣= 𝑣1 (𝑡, 𝑥, 𝑦) d𝑦
ℎ −1
∫ 𝜀𝜂  
1 1 2 2 1 2
= 𝑣 (𝑡, 𝑥) + 𝜎 𝜃 e 𝑣 𝑥𝑥 (𝑡, 𝑥) − 𝜎 (𝑦 + 1) e
2
𝑣 𝑥𝑥 (𝑡, 𝑥) d𝑦 + 𝑂(𝜎 4 )
ℎ −1
e
2 2
  𝜀𝜂
1 1 1 1
𝑣 (𝑡, 𝑥) + 𝜎2 𝜃2e
=e 𝑣 𝑥𝑥 (𝑡, 𝑥) − 𝜎 2e
𝑣 𝑥𝑥 (𝑡, 𝑥) (𝑦 + 1)3 + 𝑂(𝜎 4 )
2 2 ℎ 3 −1
 
1 1
𝑣 (𝑡, 𝑥) + 𝜎2e
=e 𝑣 𝑥𝑥 (𝑡, 𝑥) 𝜃2 − ℎ 2 + 𝑂(𝜎 4 )
2 3
 
1 1
𝑣 (𝑡, 𝑥) + 𝜎2e
=e 𝑣 𝑥𝑥 (𝑡, 𝑥) 𝜃2 − + 𝑂(𝜀𝜎 2 ) + 𝑂(𝜎 4 ).
2 3

Inserting this into the classical Boussinesq system (3.8) yields

 
1 2 2 1
𝜂𝑡 + (1 + 𝜀𝜂)e 𝑣 𝑥+ 𝜎 𝜃 − 𝑣 𝑥𝑥𝑥 = 𝑂(𝜀2 ) + 𝑂(𝜀𝜎 2 ) + 𝑂(𝜎 4 ),


e
2 3
1  
𝑣 𝑡 + 𝜂 𝑥 + 𝜀e
e 𝑣e𝑣 𝑥 + 𝜎2 𝜃2 − 1 e𝑣 𝑡𝑥𝑥 = 𝑂(𝜀2 ) + 𝑂(𝜀𝜎 2 ) + 𝑂(𝜎 4 ).
2

23
Next, we observe

𝜂𝑡 + e
𝑣 𝑥 = 𝑂(𝜀) + 𝑂(𝜎 2 ), 𝑣 𝑡 + 𝜂 𝑥 = 𝑂(𝜀) + 𝑂(𝜎 2 ).
e (3.9)

𝑣 𝑥𝑥𝑥 by −𝜂𝑡𝑥𝑥 without changing the

Thus, we can replace the higher-order derivative e
order of approximation, resulting in
 
1 2 2 1
𝜂𝑡 + (1 + 𝜀𝜂)e 𝑣 𝑥− 𝜎 𝜃 − 𝜂𝑡𝑥𝑥 = 𝑂(𝜀2 ) + 𝑂(𝜀𝜎 2 ) + 𝑂(𝜎 4 ),


2 3
1  
𝑣 𝑡 + 𝜂 𝑥 + 𝜀e
e 𝑣e𝑣 𝑥 + 𝜎2 𝜃2 − 1 e𝑣 𝑡𝑥𝑥 = 𝑂(𝜀2 ) + 𝑂(𝜀𝜎 2 ) + 𝑂(𝜎 4 ).
2
Finally, we choose1
2
𝜃2 =
3
to obtain
1
𝜂𝑡 + (1 + 𝜀𝜂)e
𝑣
− 𝜎 2 𝜂𝑡𝑥𝑥 = 𝑂(𝜀2 ) + 𝑂(𝜀𝜎 2 ) + 𝑂(𝜎 4 ),


6 𝑥
1 2
𝑣 𝑡 + 𝜂 𝑥 + 𝜀e
e 𝑣e𝑣𝑥 − 𝜎 e𝑣 𝑡𝑥𝑥 = 𝑂(𝜀2 ) + 𝑂(𝜀𝜎 2 ) + 𝑂(𝜎 4 ).
6
Neglecting higher-order terms 𝑂(𝜀2 ) + 𝑂(𝜀𝜎2 ) + 𝑂(𝜎4 ) and writing this explicitly in
terms of non-dimensionalized variables, we have
1
𝜂𝑡∗∗ + (1 + 𝜀𝜂∗ )e
𝑣∗ − 𝜎 2 𝜂𝑡∗∗ 𝑥 ∗ 𝑥 ∗ = 0,


6 𝑥∗
1
𝑣 𝑡∗∗ + 𝜂∗𝑥 ∗ + 𝜀e
e 𝑣∗ e
𝑣 ∗𝑥 ∗ − 𝜎 2e
𝑣 𝑡∗∗ 𝑥 ∗ 𝑥 ∗ = 0.
6
Going back to dimensional variables using the transformations
𝜎𝑔 1 1 1
𝑡∗ = 𝑡, 𝑥∗ = 𝑥, 𝜂∗ = 𝜂, 𝑣∗ = 𝑣,
𝑐0 𝜆 𝜀ℎ 0 𝜀𝑐0
e e

with 𝑐0 = 𝑔 ℎ 0 and 𝜎 = ℎ0 /𝜆, we compute

p

1 𝑐0 1
𝜂𝑡∗∗ = 𝜂𝑡 = 𝜂𝑡 ,
𝜀ℎ 0 𝜎𝑔 𝑐0 𝜀𝜎
1 1
𝑣 ∗𝑥 ∗ = 𝜆e𝑣𝑡 = ℎ0e𝑣𝑥 ,
𝜀𝑐0 𝑐0 𝜀𝜎
e
1 1 1
𝜀(𝜂∗e
𝑣 ∗ )𝑥 ∗ = 𝜀 𝜆(𝜂e 𝑣 )𝑥 = 𝑣 )𝑥 ,
(𝜂e
𝜀ℎ 0 𝜀𝑐 0 𝑐0 𝜀𝜎
1 𝑐0 2 1 2
𝜎 2 𝜂𝑡∗∗ 𝑥 ∗ 𝑥 ∗ = 𝜎2 𝜆 𝜂𝑡𝑥𝑥 = ℎ 𝜂𝑡𝑥𝑥 ,
𝜀ℎ 0 𝜎𝑔 𝑐0 𝜀𝜎 0
1 𝑐0 1
𝑣 𝑡∗∗ = 𝑣𝑡 = 𝑣𝑡 ,
𝜀𝑐0 𝜎𝑔 𝜀𝜎𝑔
e e e
1 1
𝜂∗𝑥 ∗ = 𝜆𝜂 𝑥 = 𝑔𝜂 𝑥 ,
𝜀ℎ 0 𝜀𝜎𝑔
1 1
𝑣∗ e
𝜀e 𝑣 ∗𝑥 ∗ = 𝜀 𝜆 𝑣 𝑣 𝑥 = 𝑣e 𝑣𝑥 ,
𝜀2 𝑐02 𝜀𝜎𝑔
e e e

1 𝑐0 2 1 2
𝑣 𝑡∗∗ 𝑥 ∗ 𝑥 ∗ =
𝜎 2e 𝜎2 𝜆e 𝑣 𝑡𝑥𝑥 = ℎe 𝑣 𝑡𝑥𝑥 ,
𝜀𝑐0 𝜎𝑔 𝜀𝜎𝑔 0
1 For a general 𝜃 ∈ [0, 1], we obtain the so-called 𝑎𝑏𝑐𝑑-Boussinesq system.

24
we finally obtain the BBM-BBM system

1
𝜂𝑡 + (𝜂 + ℎ 0 )e
𝑣 − ℎ02 𝜂𝑡𝑥𝑥 = 0,


6 𝑥
1 2
𝑣 𝑡 + 𝑔𝜂 𝑥 + e
e 𝑣e𝑣 𝑥 − ℎ 0e
𝑣 𝑡𝑥𝑥 = 0,
6
where 𝜂pis the total water height above the flat bathymetry at 𝑦 = −ℎ 0 , e
𝑣 is the velocity
at 𝑦 = ( 2/3 − 1)ℎ 0 , and the still water level is at 𝑦 = 0.

3.5 The shallow water (Saint-Venant) equations

The classical shallow water equations (or Saint-Venant equations, in particular in the
French literature) are derived from (3.6) by neglecting terms of order 𝑂(𝜎2 ). Thus, we
obtain
𝜕𝑡 ℎ + (ℎ𝑣)𝑥 = 0,
(3.10)
𝑣 𝑡 + 𝑔 ℎ 𝑥 + 𝑣 𝑣 𝑥 = 0,
in dimensional variables, which is the same as the SGN equations (3.7) without the
dispersive term. We can also write the shallow water equations in the classical conser-
vative form by observing

(ℎ𝑣)𝑡 = ℎ 𝑡 𝑣 + ℎ𝑣 𝑡
= −ℎ 𝑥 𝑣 𝑣 − ℎ𝑣 𝑥 𝑣 − 𝑔 ℎ ℎ 𝑥 − ℎ𝑣 𝑣 𝑥
 
1
= − ℎ𝑣 + 𝑔 ℎ 2
2
.
2 𝑥

Thus, we obtain the more familiar conservative form

𝜕𝑡 ℎ + 𝜕𝑥 ℎ𝑣 = 0,


 
1 (3.11)
𝜕𝑡 ℎ𝑣 + 𝜕𝑥 ℎ𝑣 + 𝑔 ℎ 2 = 0,

2
2

of the 1D shallow water equations with flat bathymetry.

→ end of lecture 2025-05-08

3.6 Unidirectional models: KdV and BBM

Recall the intermediate results

𝜂𝑡 + 𝑣 𝑥 = 𝑂(𝜀) + 𝑂(𝜎 2 ), 𝑣 𝑡 + 𝜂 𝑥 = 𝑂(𝜀) + 𝑂(𝜎 2 ),

we obtained when deriving the BBM-BBM system. Differentiating the first one with
respect to time and inserting the second equation, we obtain the wave equation

𝜂𝑡𝑡 = 𝜂 𝑥𝑥 + 𝑂(𝜀) + 𝑂(𝜎2 ).

Since
(𝜕𝑡2 − 𝜕𝑥2 ) = (𝜕𝑡 − 𝜕𝑥 )(𝜕𝑡 + 𝜕𝑥 ),
we take the solution of the wave equation traveling to the right, i.e.,

𝜂𝑡 + 𝜂 𝑥 = 𝑂(𝜀) + 𝑂(𝜎 2 ).

25
We make the ansatz
𝑣 = 𝜂 + 𝜀𝐴 + 𝜎 2 𝐵 + 𝑂(𝜀2 ) + 𝑂(𝜎 4 )
and insert it into the classical Boussinesq system (3.8)

𝜂𝑡 + 𝑣 𝑥 + 𝜀(𝜂𝑣)𝑥 = 0,
1
𝑣 𝑡 + 𝜂 𝑥 + 𝜀𝑣 𝑣 𝑥 − 𝜎 2 𝑣 𝑡𝑥𝑥 = 𝑂(𝜀2 ) + 𝑂(𝜀𝜎 2 ) + 𝑂(𝜎 4 ),
3
to obtain

𝜂𝑡 + 𝜂 𝑥 + 𝜀(𝐴 𝑥 + 2𝜂𝜂 𝑥 ) + 𝜎 2 𝐵 𝑥 = 𝑂(𝜀2 ) + 𝑂(𝜎 4 ),

 
1
𝜂𝑡 + 𝜂 𝑥 + 𝜀(𝐴𝑡 + 𝜂𝜂 𝑥 ) + 𝜎 2
𝐵𝑡 − 𝜎 2 𝜂𝑡𝑥𝑥 = 𝑂(𝜀2 ) + 𝑂(𝜀𝜎 2 ) + 𝑂(𝜎 4 ).
3

Equating the terms of the same orders to make both equations equivalent and using
𝐴𝑡 + 𝐴 𝑥 = 𝑂(𝜀) + 𝑂(𝜎 2 ), we obtain
1 1
𝐴 𝑥 = − 𝜂𝜂 𝑥 i.e., 𝐴 = − 𝜂2
2 4
and analogously
1
𝐵 = 𝜂 𝑥𝑥 .
6
Thus, ignoring higher-order terms, the ansatz becomes

1 1
𝑣 = 𝜂 − 𝜀𝜂2 + 𝜎 2 𝜂 𝑥𝑥 ,
4 6
leading to
1 1
𝜂𝑡 + 𝜂 𝑥 − 𝜀(𝜂2 )𝑥 + 𝜎 2 𝜂 𝑥𝑥𝑥 + 𝜀(𝜂2 )𝑥 = 𝑂(𝜀2 ) + 𝑂(𝜎 4 ).
4 6
Ignoring higher-order terms and going back to dimensional variables, we obtain the
dimensional KdV equation

𝑔 ℎ3 𝑔
r r
p 3
𝜂𝑡 + 𝑔 ℎ0𝜂𝑥 + 𝜂𝜂 𝑥 + 0 𝜂 𝑥𝑥𝑥 = 0. (3.12)
2 ℎ0 6 ℎ0

Similarly, we can again use 𝜂𝑡 + 𝜂 𝑥 = 𝑂(𝜀) + 𝑂(𝜎 2 ) to turn the dispersive third-
derivative term into a mixed derivative, leading to the dimensional BBM equation

𝑔 ℎ 02
r
p 3
𝜂𝑡 + 𝑔 ℎ0𝜂𝑥 + 𝜂𝜂 𝑥 − 𝜂𝑡𝑥𝑥 = 0. (3.13)
2 ℎ0 6

26