![MathJax Logo](/templates/jsp/_style2/_tandf/pb2/images/math-jax.gif)
Abstract
We study the propagation of two-dimensional tsunami waves triggered by a seaquake in the open sea in the presence of underlying wind-generated currents, corresponding to background flows of constant vorticity. A suitable scaling of the governing equations introduces dimensionless parameters, of particular interest being the setting of linear waves that only depend on the vertical movement of the sea bed. We use Fourier analysis methods to extract formulae for the function f which describes the vertical displacement of the water's free surface. We show that the results are particularly useful in the physically relevant shallow-water regime: in the irrotational case the predictions fit well with the observed behaviour of some historical tsunamis. In other situations, the stationary-phase principle gives insight into the asymptotic behaviour of f.
Keywords:
1. Introduction
The aim of this paper is to pursue a mathematical analysis of tsunami waves triggered by an earthquake in the open sea.
Tsunamis are water waves which are generated by an impulsive disturbance that vertically displaces the water column, most commonly by tectonic activity, but sometimes by other means as well, such as underwater landslides or asteroid impact. Hundreds of tsunamis were confirmed since 1900, about 80% of which occurred in the Pacific Ocean. Most of these tsunami waves have small amplitude, being nondestructive and barely detectable in the ocean (less than 0.3 m high, the average wave height in the Pacific Ocean being about 2.5 m). Large, destructive tsunamis, causing loss of life and major coastal destruction, take place typically once every several decades. The highest tsunami wave that was ever measured exceeded 500 m and occurred on 9 July 1958, in the Lituya Bay inlet of the Gulf of Alaska, caused by a massive rock dislocated by an earthquake and plunging 900 m down into the water (see [Citation1]). Let us also mention that 66 million years ago, an asteroid with a diameter of about 10 km struck the waters of the Gulf of Mexico near the Yucatán Peninsula, and it is estimated that the impact generated tsunami waves more than 1.5 km high that caused sea levels to rise in all corners of the Earth (see [Citation2]); the subsequent cataclysmic events – plumes of aerosol, soot and dust filling the air, and wildfires starting as flaming pieces of material blasted from the impact re-entered the atmosphere and rained down – ended the era of the dinosaurs and triggered a mass extinction of about 75% of animal and plant life on Earth. Since it appears [Citation3] that up to 75% of all recorded tsunamis were generated by undersea earthquakes, in this paper we focus on this type. Three events stand out in the last 100 years: the 22 May 1960 Chile tsunami (caused by the largest earthquake ever recorded and with waves propagating across the Pacific Ocean, causing havoc in Hawaii and in Japan), the 26 December 2004 tsunami (that killed more than 200,000 people around the shores of the Indian Ocean), and the tsunami off the coast of Japan on 11 March 2011 (that killed thousands of people and triggered a nuclear accident). In all three cases the earthquake was undersea, along a fault line, in which case the generated waves are typically two-dimensional and with long wavelengths, of the order of 100 km (see the discussions in [Citation1,Citation3,Citation4]). Given that the average ocean depth is about 4 km, this means that tsunami waves in the open ocean are shallow water waves. While near the shore, where the water depth gradually diminishes, nonlinear effects become dominant (and induce wave-breaking), a long-standing issue was whether in the open sea the dominant linear behaviour of the flow is greatly affected by the accumulation of nonlinear factors. Some authors (see [Citation5–7]) advocated that KdV-like models apply, but even in the case of propagation over the maximal distance that is possible (across the Pacific Ocean, which was the case for the 1960 Chile tsunami), the scale analysis reveals that weakly nonlinear theory is not adequate (see the discussion in [Citation1,Citation3,Citation8–10]): the effects remain negligible in the open sea. Advances in the modelling of the tsunami wave propagation in the open sea, triggered by an undersea earthquake, were obtained in [Citation4,Citation11] for the setting of a still sea prior to the perturbation. In this paper we extend these considerations to accommodate the presence of underlying currents. These are typically modelled by constant-vorticity flows (see the discussion in [Citation12,Citation13]) but even the irrotational setting (flows with zero vorticity) might admit uniform tidal currents.
2. The governing equations
We consider an analogous setup as in [Citation11]: we consider two-dimensional surface waves with the unbounded X-direction corresponding to the direction of wave propagation. The water's free surface is , where d is the average depth of the sea, the seabed being
. At time T = 0, we assume that
and
. Then the seabed moves for a short time in some compact region
and is flat after that again. We are ultimately interested in the deformation of the free surface, i.e. in
. To do this, we study a partial differential equation that describes its evolution in time. Consider
the velocity field of the flow in rectangular Cartesian coordinates
. We assume that the water is inviscid and that the resulting flow is with constant vorticity. Note that the presence of non-uniform currents makes the hypothesis of irrotational flow inappropriate and, since in a setting in which the waves are long compared to the water depth, the existence of a non-zero mean vorticity is important rather than its specific distribution [cf. Citation12]. This is not merely a mathematical simplification: wind-generated currents in shallow waters with nearly flat beds have been shown to be accurately described as flows with constant vorticity [Citation13]. We now present the governing equations for water flows [Citation1]. We have the incompressible Euler equations
(1)
(1)
in
, where ρ is the constant density, g is the constant acceleration of gravity and P is the pressure. The effects of surface tension are negligible for wavelengths greater than a few centimetres [Citation14,Citation15], hence the major factor governing the wave motion is the balance between gravity and the inertia of the system. We also have the kinematic boundary conditions
(2)
(2)
and
(3)
(3)
which ensure that particles on these boundaries are confined to them at all times. We also have the dynamic boundary condition
(4)
(4)
where
is the constant atmospheric pressure. This decouples the motion of the water from the motion of the air above it [Citation16]. The system (Equation1
(1)
(1) )–(Equation4
(4)
(4) ) is to be solved with the following initial conditions
(5)
(5)
These initial conditions express the assumption that initially, at time T = 0, there is only a horizontal flow whose speed depends linearly on the depth. This is a generalisation of [Citation11], where it was assumed that the water is completely at rest at T = 0, a case which corresponds to
. A two-dimensional water flow that has constant vorticity at some instant will have this feature at all other times (see the discussion in Chapter 1 of [Citation1]). By the initial condition (Equation5
(5)
(5) ) this implies
(6)
(6) for
. Figure illustrates our set-up in the case
.
The first step in the analysis of the system (Equation1(1)
(1) )–(Equation6
(6)
(6) ) will be to non-dimensionalize and scale the equations. First, it will be convenient to introduce the non-dimensional excess pressure p relative to the hydrostatic pressure distribution by
(7)
(7) To obtain meaningful scales, we introduce the average or typical wavelength λ of the wave, we use
as a scale for the wave speed and
as a time scale. We use the change of variables
(8)
(8) with
(9)
(9) where a is a typical, perhaps maximal, amplitude of the wave. Note that (Equation9
(9)
(9) ) should be interpreted as ensuring that the variations of the wave and of the seabed are of comparable size.
We introduce the parameters
(10)
(10) where ε measures the relative size of the amplitude to the average water depth and δ measures the average water depth to the wavelength. Now we get the non-dimensional equations
(11)
(11) The magnitudes of ε and δ correspond to the different general types of water wave problem: The limits
and
produce the shallow water and deep water regime, respectively. On the other hand,
corresponds to regime of waves of small amplitude (see the discussion in [Citation16,Citation17]).
We want to linearise the problem and in this case we will want to let and keep δ fixed. But the system (Equation11
(11)
(11) ) shows that v and p are of order ε and hence also u. Hence we use the following scaling:
(12)
(12) where we avoid a new notation. We then obtain the system
(13)
(13) in
and
(14)
(14) We now linearise this system by taking
and obtain
(15)
(15) By the first equation in (Equation15
(15)
(15) ), there exists a stream function ψ such that
(16)
(16) If we plug this in (Equation15
(15)
(15) ), we get
(17)
(17) Note that we want to eventually find f. But p and ψ are also unknown, only h is given. We will now derive equations only involving ψ and h.
We first differentiate the second equation in (Equation17(17)
(17) ) with respect to x and t, respectively, to obtain
(18)
(18) and
(19)
(19) Furthermore, we can differentiate the fifth equation in (Equation17
(17)
(17) ) with respect to x to obtain
(20)
(20) Observe that the fourth equation in (Equation17
(17)
(17) ) implies that we can exchange the derivatives of f in (Equation20
(20)
(20) ) by the corresponding derivatives of p. But then we can insert Equations (Equation18
(18)
(18) ) and (Equation19
(19)
(19) ) in (Equation20
(20)
(20) ) to obtain
on y = 1. We introduce the constant
(21)
(21) which is the non-dimensional horizontal velocity of the flow at the surface. This allows us to rewrite the previous equation as
on y = 1. By observing the quadratic expression on the right-hand side and introducing the differential operator
(22)
(22) we can rewrite this again as
So we obtain
(23)
(23) which only involves ψ and h. We can then finally recover f by the following procedure: We know by the fourth and fifth equation in (Equation17
(17)
(17) ) that
on y = 1 and hence
(24)
(24) on y = 1. Together with the initial condition
we can recover f by
and p can be reconstructed similarly.
3. General solution formulae for linear waves
We will consider the space and space-time Fourier transform with notations
and
We also define Fourier multipliers: Let
be some function. We define
by
or equivalently
where
, whenever this is defined. Note that
maps real-valued functions to real-valued functions whenever the condition
(25)
(25) is satisfied. In particular, this is satisfied if m is real-valued and even. We can apply the space-time Fourier transform to the first three equations of system (Equation23
(23)
(23) ) to obtain
(26)
(26) where
(27)
(27) with the relation
(28)
(28) Since (Equation26
(26)
(26) ) is just a linear ODE of second order in y for fixed ξ and ω, we can, somewhat explicitly, write down the solution:
(29)
(29) where
are determined by the boundary conditions. We can then reconstruct
by applying the Fourier transform on the first equation in (Equation24
(24)
(24) ) to obtain
and hence
(30)
(30) We ignore the potential singularity for the time being. We know that the functions
and
are uniquely determined. Once we have formulae for them, we get a formula for
by inserting in (Equation29
(29)
(29) ). Then we can insert this formula in (Equation30
(30)
(30) ) to obtain an expression for
depending on
.
We insert (Equation29(29)
(29) ) in (Equation26
(26)
(26) ) to obtain
and
This linear equation in
and
is solved easily and by (Equation30
(30)
(30) ) we have
which we will write as
(31)
(31) We will also write this as
(32)
(32) One would need to consider the roots of the denominator in (Equation31
(31)
(31) ) in order to justify (Equation32
(32)
(32) ) rigorously. However for the case A = 0, it turns out that the singularity cancels.
In the case A = 0 formula (Equation32(32)
(32) ) has the simpler form
(33)
(33) We will extract a formula for f whenever it solves the following type of equation. Let
be an even function such that
for all
. Note that in particular (Equation25
(25)
(25) ) is satisfied. We try to derive a formula for the solution f of
(34)
(34) where θ is a given forcing term. We will assume f and θ to be in
and the limit in (Equation34
(34)
(34) ) is adequate for functions in the Schwartz class [Citation18]. We can decompose the operator
into
(35)
(35) since S and
commute. If we now set
then we see by (Equation35
(35)
(35) ) and the fact that we can commute the two factors that
We can now filter
by the group
: Let
. A calculation shows that
This is an inhomogeneous transport equation with solution
This yields
If we write out the Fourier multipliers explicitly, we get
(36)
(36) We will deal with the possible singularity in the for us relevant case: We have
and
. We see that τ is even and
for
. We have furthermore that τ is smooth with derivatives of at most polynomial growth at infinity. We will use that
(37)
(37) Now observe that the function
is real-valued and hence we can subtract
in (Equation36
(36)
(36) ) to obtain
(38)
(38) Clearly,
remains bounded around 0. For the quotient, note that we can rewrite it as follows:
The second factor converges to 1 by (Equation37
(37)
(37) ). For the first factor, note that the function
is left and right differentiable in 0 and hence the first factor remains bounded as
for any s. We conclude that (Equation38
(38)
(38) ) holds (at least pointwise).
Note that we cannot do the same for the case : The operator on the left hand side of (Equation32
(32)
(32) ) is given by
and there does not seem to be an apparent way to obtain a root for
In fact, it is unclear whether this operator is even positive.
4. Behaviour in different regimes
We want to investigate some regimes for which our considerations give useful predictions. We will consider on the one hand the case where δ is small, which is justified considering the model as . On the other hand, we will discuss the case where δ is finite and non-vanishing. The case where δ is large (corresponding to the formal limit
) is not relevant for us since deep-water tsunamis are rare.
4.1. The shallow water regime
For simplicity, we use the model in the case A = 0: As , (Equation33
(33)
(33) ) becomes
(39)
(39) The initial conditions are on the one hand
from (Equation15
(15)
(15) ). On the other hand, this gives
and inserting this in the sixth relation in (Equation15
(15)
(15) ) for t = 0 yields
. By using Duhamel's principle (similar to [Citation19, p. 80–81]), we get the formula
(40)
(40) We will assume that h can be separated in the following way:
for
,
such that
for
and
for
(where
represents the duration of the earthquake), and with
for
modelling a localised tsunami source. Inserting this into (Equation40
(40)
(40) ) yields the formula
which is a lot simpler. In the limiting case
, a becomes the Heaviside step function (modelling an instantaneous upward thrust of the seabed near to the earthquake's epicentre) and the formula further simplifies to
(41)
(41) for
and t>0. We can draw some insightful conclusions:
at each instant t>0 after initiation, the generated wave is localised as
for
;
the surface wave consists of one stationary part (which can be disregarded, since
) and two travelling waves, one moving to the right and the other moving to the left;
the wave travelling to the left moves with non-dimensionalized speed 1−C (corresponding to
in the physical variables) and the wave travelling to the right moves with non-dimensionalized speed 1 + C (corresponding to
in the physical variables);
the shapes of the waves travelling to the left and to the right remain unchanged and are precisely that of the bed deformation at a scale of
and
, respectively.
This behaviour is illustrated in Figures and . On the one hand, the wave travelling to the right is faster for bigger C, but the maximal amplitude decreases. On the other hand, the wave travelling to the left becomes slower for bigger C, but the maximal amplitude increases. For B = 0 and hence C = 0, this model gives good predictions for the maximal amplitude and speed of the propagating wave for the two largest historical tsunamis: the December 2004 and the May 1960 tsunamis, see [Citation11]. Since , i.e.
, the model for
will give a very similar prediction, which might be slightly more precise, if one were in a situation, where B is known somewhat precisely, say due to a current.
4.2. The finite, non-vanishing regime
Here we cannot expect to be able to write down f in any explicit fashion. However, we will discuss a method to analyse the asymptotic behaviour of f in (Equation38(38)
(38) ). Ultimately, we will see that this approach does not work, since the method only works for too large t, see (Equation48
(48)
(48) ).
We cannot do the same for B = 0 or the general case, since we lack a formula which would need to correspond to (Equation36(36)
(36) ) or (Equation38
(38)
(38) ). If such a formula were available, stationary-phase analysis might then be appropriate.
4.2.1. Principle of stationary-phase
We want to analyse the asymptotic behaviour of (Equation38(38)
(38) ). For this we try to use the stationary-phase principle. Recall that we set
and
. We first rewrite (Equation38
(38)
(38) ) as
where
. The derivative of the phase factor
is given by
and one can check that
(42)
(42) and
(43)
(43) One computes the limits
and concludes that
is strictly decreasing on
from the asymptotic value 1 towards the asymptotic value 0 and on
from the asymptotic value 0 to the asymptotic value
. So we conclude that
is possible if
for exactly one
for given
. So as long as
is in the corresponding range, we can write
. The stationary-phase principle gives
for large t, where σ is the sign of
, which is just
by (Equation43
(43)
(43) ). We see that this is a Fourier transform (with respect to time) and hence get
We have
and hence get
(44)
(44) Finally, we try to extract the asymptotics with respect to δ: We fix a
. A look at (Equation42
(42)
(42) ) gives that
is
and clearly even in δ. Since we have
, the Taylor expansion with respect to δ around 0 has the form
for some coefficient
. One could find this coefficient through tedious calculations, but one can instead notice that
is the coefficient of
of the expansion of
. We have
and hence we conclude
. We conclude that for
the stationary-phase point
with
satisfies
. Inserting this into (Equation43
(43)
(43) ) yields
(45)
(45)
4.2.2. Some typical physical parameters
We want to see whether stationary-phase analysis, using the formula
could be justified here. Since the stationary-phase principle applies for large t we insert typical values of the physical parameters for tsunamis propagating at open sea (see [Citation8]):
This leads to
and
which point to linear wave approximation. Applying the stationary-phase principle is justified if
(46)
(46) We have by (Equation45
(45)
(45) ),
, which translates to
(47)
(47) In physical variables this means
(48)
(48) which takes way too long to be applicable. If one were to relax (Equation47
(47)
(47) ) to a condition of the form
one would arrive at
which would barely be applicable in the case where the epicentre of the seaquake is far off the shore. Nevertheless, relaxing (Equation47
(47)
(47) ) is quite dubious. We conclude that using stationary-phase analysis is not an effective tool in the case A = 0.
Disclosure statement
No potential conflict of interest was reported by the author(s).
References
- Constantin A. Nonlinear water waves with applications to wave-current interactions and tsunamis. In: CBMS-NSF regional conference series in applied mathematics, vol. 81. Philadelphia (PA): SIAM; 2011.
- Range MM, Arbiv BK, Johnson BC, et al. The Chicxulub impact produced a powerful global tsunami. AGU Adv. 2022;3:Article ID e2021AV00627.
- Arcas D, Segur H. Seismically generated tsunamis. Phil Trans Roy Soc London A. 2012;370:1505–1542.
- Dias F, Dutykh D. Water waves generated by a moving bottom. In: Kundu A, editor. Tsunami and nonlinear waves. Berlin: Springer; 2007. p. 65–95.
- Craig W. Surface water waves and tsunamis. J Dynam Diff Eq. 2006;18:525–549.
- Lakshmanan M. Integrable nonlinear wave equations and possible connections to tsunami dynamics. In: Kundu A, editor. Tsunami and nonlinear waves. Springer; 2007. p. 31–49.
- Segur H. Integrable nonlinear wave equations and possible connections to tsunami dynamics. In: Kundu A, editor. Tsunami and nonlinear waves. Springer; 2007. p. 3–29.
- Constantin A. On the relevance of soliton theory to tsunami modelling. Wave Motion. 2009;46:420–426.
- Constantin A, Johnson RS. Propagation of very long water waves, with vorticity, over variable depth, with applications to tsunamis. Fluid Dynam Res. 2008;40:175–211.
- Stuhlmeier R. KdV theory and the Chilean tsunami of 1960. Discrete Cont Dyn Syst Ser B. 2009;12:623–632.
- Constantin A, Germain P. On the open sea propagation of water waves generated by a moving bed. Phil Trans Roy Soc London A. 2012;370:1587–1601.
- da Silva AFT, Peregrine DH. Steep, steady surface waves on water of finite depth with constant vorticity. J Fluid Mech. 1988;195:281–302.
- Ewing JA. Wind, wave and current data for the design of ships and offshore structures. Mar Struct. 1990;3:421–459.
- Barnett TP, Kenyon KE. Recent advances in the study of wind waves. Rep Progr Phys. 1975;38:667–729.
- Lighthill MJ. Waves in fluids. Cambridge, UK: Cambridge University Press; 1978.
- Johnson RS. A modern introduction to the mathematical theory of water waves. Cambridge, UK: Cambridge University Press; 1997.
- Constantin A, Johnson RS. On the non-dimensionalisation, scaling and resulting interpretation of the classical governing equations for water waves. J Nonlinear Math Phys. 2008;15:58–73.
- Strichartz RS. A guide to distribution theory and Fourier transforms. Boca Raton (FL): CRC Press; 1994.
- Evans LC. Partial differential equations. Providence (RI): Amer Math Soc; 2010.