As we are all now very much aware, tsunamis are water waves that start in the deep ocean, usually because of an underwater earthquake (though tsunamis can also be caused by underwater landslides or volcanoes), and then propagate towards shore. Initially, tsunamis have relatively small amplitude (a metre or so is typical), which would seem to render them as harmless as wind waves. And indeed, tsunamis often pass by ships in deep ocean without anyone on board even noticing.

However, being generated by an event as large as an earthquake, the wavelength of the tsunami is huge – 200 kilometres is typical (in contrast with wind waves, whose wavelengths are typically closer to 100 metres). In particular, the wavelength of the tsunami is far greater than the depth of the ocean (which is typically 2-3 kilometres). As such, even in the deep ocean, the dynamics of tsunamis are essentially governed by the shallow water equations. One consequence of these equations is that the speed of propagation {v} of a tsunami can be approximated by the formula

\displaystyle v \approx \sqrt{g b} \ \ \ \ \ (1)

where {b} is the depth of the ocean, and {g \approx 9.8 ms^{-2}} is the force of gravity. As such, tsunamis in deep water move very fast – speeds such as 500 kilometres per hour (300 miles per hour) are quite typical; enough to travel from Japan to the US, for instance, in less than a day. Ultimately, this is due to the incompressibility of water (and conservation of mass); the massive net pressure (or more precisely, spatial variations in this pressure) of a very broad and deep wave of water forces the profile of the wave to move horizontally at vast speeds. (Note though that this is the phase velocity of the tsunami wave, and not the velocity of the water molecues themselves, which are far slower.)

As the tsunami approaches shore, the depth {b} of course decreases, causing the tsunami to slow down, at a rate proportional to the square root of the depth, as per (1). Unfortunately, wave shoaling then forces the amplitude {A} to increase at an inverse rate governed by Green’s law,

\displaystyle A \propto \frac{1}{b^{1/4}} \ \ \ \ \ (2)

at least until the amplitude becomes comparable to the water depth (at which point the assumptions that underlie the above approximate results break down; also, in two (horizontal) spatial dimensions there will be some decay of amplitude as the tsunami spreads outwards). If one starts with a tsunami whose initial amplitude was {A_0} at depth {b_0} and computes the point at which the amplitude {A} and depth {b} become comparable using the proportionality relationship (2), some high school algebra then reveals that at this point, amplitude of a tsunami (and the depth of the water) is about {A_0^{4/5} b_0^{1/5}}. Thus, for instance, a tsunami with initial amplitude of one metre at a depth of 2 kilometres can end up with a final amplitude of about 5 metres near shore, while still traveling at about ten metres per second (35 kilometres per hour, or 22 miles per hour), and we have all now seen the impact that can have when it hits shore.

While tsunamis are far too massive of an event to be able to control (at least in the deep ocean), we can at least model them mathematically, allowing one to predict their impact at various places along the coast with high accuracy. (For instance, here is a video of the NOAA’s model of the March 11 tsunami, which has matched up very well with subsequent measurements.) The full equations and numerical methods used to perform such models are somewhat sophisticated, but by making a large number of simplifying assumptions, it is relatively easy to come up with a rough model that already predicts the basic features of tsunami propagation, such as the velocity formula (1) and the amplitude proportionality law (2). I give this (standard) derivation below the fold. The argument will largely be heuristic in nature; there are very interesting analytic issues in actually justifying many of the steps below rigorously, but I will not discuss these matters here.

— 1. The shallow water wave equation —

The ocean is, of course, a three-dimensional fluid, but to simplify the analysis we will consider a two-dimensional model in which the only spatial variables are the horizontal variable {x} and the vertical variable {z}, with {z=0} being equilibrium sea level. We model the ocean floor by a curve

\displaystyle z = - b(x),

thus {b} measures the depth of the ocean at position {x}. At any time {t} and position {x}, the height of the water (compared to sea level {z=0}) will be given by an unknown height function {h(t,x)}; thus, at any time {t}, the ocean occupies the region

\displaystyle \Omega_t := \{ (x,z): -b(x) < z < h(t,x) \}.

Now we model the motion of water inside the ocean by assigning at each time {t} and each point {(x,z) \in \Omega_t} in the ocean, a velocity vector

\displaystyle \vec u(t,x,z) = (u_x(t,x,z), u_z(t,x,z)).

We make the basic assumption of incompressibility, so that the density {\rho} of water is constant throughout {\Omega_t}.

The velocity changes over time according to Newton’s second law {F=ma}. To apply this law to fluids, we consider an infinitesimal amount of water as it flows along the velocity field {\vec u}. Thus, at time {t}, we assume that this amount of water occupies some infinitesimal area {dA} and some position {\vec x(t) = (x(t), z(t))}, where we have

\displaystyle \frac{d}{dt} \vec x(t) = \vec u( t, \vec x(t) ).

Because of incompressibility, the area {dA} stays constant, and the mass of this infinitesimal portion of water is {m = \rho dA}. There will be two forces on this body of water; the force of gravity, which is {(0, -mg) = (0, -\rho) dA}, and the force of the pressure field {p(t,x,z)}, which is given by {-\nabla p dA}. At the length and time scales of a tsunami, we can neglect the effect of other forces such as viscosity or surface tension. Newton’s law {m \frac{du}{dt} = F} then gives

\displaystyle m \frac{d}{dt} \vec u(t, \vec x(t) ) = - \nabla p dA + (0,-mg)

which simplifies to the incompressible Euler equation

\displaystyle \frac{\partial}{\partial t} \vec u + (\vec u \cdot \nabla) \vec u = - \frac{1}{\rho} \nabla p + (0,-g).

At present, the pressure is not given. However, we can simplify things by making the assumption of (vertical) hydrostatic equilibrium, i.e. the vertical effect { - \frac{1}{\rho} \frac{\partial}{\partial z} p} of pressure cancels out the effect {-g} of gravity. We also assume that the pressure is zero on the surface {z=h(t,x)} of the water. Together, these two assumptions force the pressure to be the hydrostatic pressure

\displaystyle p = \rho g ( h(t,x) - z ). \ \ \ \ \ (3)

This reflects the intuitively plausible fact that the pressure at a point under the ocean should be determined by the weight of the water above that point.

The incompressible Euler equation now simplifies to

\displaystyle \frac{\partial}{\partial t} \vec u + (\vec u \cdot \nabla) \vec u = - g (\frac{\partial}{\partial x} h, 0). \ \ \ \ \ (4)

We next make the shallow water approximation that the wavelength of the water is far greater than the depth of the water. In particular, we do not expect significant changes in the velocity field in the {z} variable, and thus make the ansatz

\displaystyle \vec u(t,x,z) \approx \vec u(t,x). \ \ \ \ \ (5)

(This ansatz should be taken with a grain of salt, particularly when applied to the {z} component {u_z} of the velocity, which does actually have to fluctuate a little bit to accomodate changes in ocean depth and in the height function. However, the primary component of the velocity is the horizontal component {u_x}, and this does behave in a fairly vertically insensitive fashion in actual tsunamis.)

Taking the {x} component of (4), and abbreviating {u_x} as {u}, we obtain the first shallow water wave equation

\displaystyle \frac{\partial}{\partial t} u + u \frac{\partial}{\partial x} u = - g \frac{\partial}{\partial x} h. \ \ \ \ \ (6)

The next step is to play off the incompressibility of water against the finite depth of the ocean. Consider an infinitesimal slice

\displaystyle \{ (x,z) \in \Omega_t: x_0 \leq x \leq x_0+dx \}

of the ocean at some time {t} and position {x_0}. The total mass of this slice is roughly

\displaystyle \rho (h(t,x_0) + b(x_0)) dx

and so the rate of change of mass of this slice over time is

\displaystyle \rho \frac{\partial h}{\partial t}(t,x_0) dx.

On the other hand, the rate of mass entering this slice on the left {x=x_0} is

\displaystyle \rho u(t,x_0) (h(t,x_0) + b(x_0))

and the rate of mass exiting on the right {x=x_0+dx} is

\displaystyle \rho u(t,x_0+dx) (h(t,x_0+dx) + b(x_0+dx)).

Putting these three facts together, we obtain the equation

\displaystyle \rho \frac{\partial h}{\partial t}(t,x_0) dx = \rho u(t,x_0) (h(t,x_0) + b(x_0))

\displaystyle - \rho u(t,x_0+dx) (h(t,x_0+dx) + b(x_0+dx))

which simplifies after Taylor expansion to the second shallow water wave equation

\displaystyle \frac{\partial}{\partial t} h + \frac{\partial}{\partial x}( u (h+b) ) = 0. \ \ \ \ \ (7)

Remark 1 Another way to derive (7) is to use a more familiar form of the incompressibility, namely the divergence-free equation

\displaystyle \frac{\partial}{\partial x} u_x + \frac{\partial}{\partial z} u_z = 0. \ \ \ \ \ (8)

(Here we will refrain from applying (5) to the vertical component of the velocity {u_z}, as the approximation (5) is not particularly accurate for this component.) Also, by considering the trajectory of a particle {(x(t),h(t,x(t)))} at the surface of the ocean, we have the formulae

\displaystyle \frac{d}{dt} x(t) = u_x(x(t),h(t,x(t)))

and

\displaystyle \frac{d}{dt} h(t,x(t)) = u_z(x(t),h(t,x(t)))

which after application of the chain rule gives the equation

\displaystyle \frac{\partial}{\partial t} h(t,x) + (\frac{\partial}{\partial x} h( x )) u_x( x, h(t,x) ) = u_z( x, h(t,x) ). \ \ \ \ \ (9)

A similar analysis at the ocean floor (which does not vary in time) gives

\displaystyle -\frac{\partial}{\partial x} b( x ) u_x( x, -b(x) ) = u_z( x, -b(x) ). \ \ \ \ \ (10)

We apply these equations to the evaluation of the expression

\displaystyle \frac{\partial}{\partial x} \int_{-b(x)}^{h(t,x)} u_x( t, x, z )\ dz.

which is the spatial rate of change of the velocity flux through a vertical slice of the ocean. On the one hand, using the ansatz (5), we expect this expression to be approximately

\displaystyle \frac{\partial}{\partial x} ( u (h+b) ).

On the other hand, by differentiation under the integral sign, we can evaluate this expression instead as

\displaystyle \int_{-b(x)}^{h(t,x)} \frac{\partial}{\partial x} u_x( t, x, z)\ dz

\displaystyle + (\frac{\partial}{\partial x} h(t,x)) u_x(x, h(t,x))

\displaystyle + (\frac{\partial}{\partial x} b(x)) u_x(x, -b(x)).

If we then substitute in (8), (9), (10) and apply the fundamental theorem of calculus, one ends up with {-\frac{\partial}{\partial t}h(t,x)}, and the claim (7) follows.

The equations (6), (7) are nonlinear in the unknowns {u, h}. However, one can approximately linearise them by making the hypothesis that the amplitude of the wave is small compared to the depth of the water:

\displaystyle |h| \ll b. \ \ \ \ \ (11)

This hypothesis is fairly accurate for tsunamis in the deep ocean, and even for medium depths, but of course is not reasonable once the tsunami has reached shore (where the dynamics are far more difficult to model).

The hypothesis (11) already simplifies (7) to (approximately)

\displaystyle \frac{\partial}{\partial t} h + \frac{\partial}{\partial x}( u b ) = 0. \ \ \ \ \ (12)

As for (6), we argue that the second term on the left-hand side is negligible, leading to

\displaystyle \frac{\partial}{\partial t} u = - g \frac{\partial}{\partial x} h. \ \ \ \ \ (13)

To explain heuristically why we expect this to be the case, let us make the ansatz that {h} and {u} have amplitude {A}, {V} respectively, and propagate at some phase velocity {v} and wavelength {\lambda}; let us also make the (reasonable) assumption that {b} varies much slower in space than {u} does (i.e. that {b} is roughly constant at the scale of the wavelength {\lambda}), so we may (for a first approximation) replace {\frac{\partial}{\partial x}(ub)} by {b \frac{\partial}{\partial x} u}. Heuristically, we then have

\displaystyle \frac{\partial}{\partial x} u = O(V/\lambda)

\displaystyle \frac{\partial}{\partial x} h = O(A/\lambda)

\displaystyle \frac{\partial}{\partial t} u = O(vV/\lambda)

\displaystyle \frac{\partial}{\partial t} h = O(vA/\lambda)

and equation (12) then suggests

\displaystyle vA/\lambda \approx V b/\lambda. \ \ \ \ \ (14)

From (11) we expect {A \ll b}, and thus {v \gg V}; the wave propagates much faster than the velocity of the fluid. In particular, we expect {u \frac{\partial}{\partial x} u = O(V^2 / \lambda)} to be much smaller than {\frac{\partial}{\partial t} u = O( v V/\lambda)}, which explains why we expect to drop the second term in (6) to obtain (13).

If we now insert the above ansatz into (13), we obtain

\displaystyle v V/\lambda \approx g A / \lambda;

combining this with (14), we already get the velocity relationship (1).

Remark 2 One can also obtain (1) more quickly (up to a loss of a constant factor) by dimensional analysis, together with some additional physical arguments. Indeed, it is clear from a superficial scan of the above discussion that the velocity {v} is only going to depend on the quantities {\rho, g, b, A, V, \lambda}. As the density {\rho} is the only input that involves mass in its units, dimensional analysis already rules out any role for {\rho}. As we are in the small amplitude regime (11), we expect the dynamics to be linearised, and thus not dependent on amplitude; this rules out {A} (and similarly {V}, which is the amplitude of the velocity field, and which is negligible when compared against the phase velocity {V}). Finally, in the long wavelength regime {\lambda \gg b}, we expect the wavelength to be physically undetectable at local scales (it requires not only knowledge of the slope of the height function at one’s location, but also the second derivative of that function (i.e. the curvature of the ocean surface), which is lower order). So we rule out dependence on {\lambda} also, leaving only {g} and {b}, and at this point dimensional analysis forces the relationship (1) up to constants. (Unfortunately, I do not know of an analogous dimensional analysis argument that gives (2).)

To get the relation (2), we have to analyse the ansatz a bit more carefully. First, we combine (13) and (12) into a single equation for the height function {h}. Indeed, differentiating (12) in time and then substituting in (13) and (1) gives

\displaystyle \frac{\partial^2}{\partial t^2} h - \frac{\partial}{\partial x}( v^{2} \frac{\partial}{\partial x} h ) = 0.

To solve this wave equation, we use a standard sinusoidal ansatz

\displaystyle h(t,x) = A(t,x) \sin( \phi( t, x ) / \epsilon )

where {A, \phi} are slowly varying functions, and {\epsilon > 0} is a small parameter. Inserting this ansatz and extracting the top order terms in {\epsilon}, we conclude the eikonal equation

\displaystyle \phi_t^2 - v^{2} \phi_x^2 = 0

and the Hamilton-Jacobi equation

\displaystyle 2 A_t \phi_t + A \phi_{tt} - v^2 (2A_x \phi_x + A \phi_{xx}) - 2vv_x A \phi_x = 0.

From the eikonal equation we see that {\phi} propagates at speed {v}. Assuming rightward propagation, we thus have

\displaystyle \phi_t = -v \phi_x. \ \ \ \ \ (15)

As for the Hamilton-Jacobi equation, we solve it using the method of characteristics. Multiplying the equation by {A}, we obtain

\displaystyle (A^2 \phi_t)_t - v^2 (A^2 \phi_x)_x - 2vv_x A^2 \phi_x = 0.

Inserting (15) and writing {F := A^2 \phi_x}, one obtains

\displaystyle - v F_t - v^2 F_x - 2 v v_x F = 0

which simplifies to

\displaystyle (\partial_t + v \partial_x) (v^2 F) = 0.

Thus we see that {v^2 F} is constant along characteristics. On the other hand, differentiating (15) in {x} we see (after some rearranging) that

\displaystyle (\partial_t + v \partial_x) ( v \phi_x ) = 0

so {v \phi_x} is also constant along characteristics. Dividing, we see that {A^2 v} is constant along characteristics, leading to the proportionality relationship

\displaystyle A \propto \frac{1}{\sqrt{v}}

which gives (2).

Remark 3 It becomes difficult to retain the sinusoidal ansatz once the amplitude exceeds the depth, as it leads to the absurd conclusion that the troughs of the wave lie below the ocean floor. However, a remnant of this effect can actually be seen in real-life tsunamis, namely that if the tsunami starts with a trough rather than a crest, then the water at the shore draws back at first (sometimes for hundreds of metres), before the crest of the tsunami hits. As such, the sudden withdrawal of water of a shore is an important warning sign of an immediate tsunami.