## Abstract

Within the framework of linear theory, applicable far from the shore, we investigate the two-dimensional propagation of waves generated in the ocean by a sudden seabed deformation.

## 1. Introduction

Several recent investigations are devoted to the generation of tsunami waves caused by a sudden seabed deformation owing to an earthquake [1–5]. In contrast to earlier approaches, which assume that the free-surface deformation at initiation is identical to the seabed deformation after the occurrence of earthquake, these studies take into account the dynamics of seafloor displacement over a short period of time and do not consider the initial wave form as instantly generated. Instead, the water is still prior to the earthquake and a transient coupling occurs between the free-surface displacement and the seabed deformation. Within the framework of linear irrotational theory, explicit integral formulae were obtained using Fourier transforms (or combined Fourier and Laplace transforms) and asymptotic analysis, and/or numerical computations were performed. However, while the setting of linear theory is appropriate for waves of small amplitude, the hypothesis of irrotational flow is justified only if the length scales of the displacements in the horizontal and vertical directions are comparable. While such a scenario is definitely realistic, it is of interest to investigate the more general setting in which these length scales are different. The aim of this paper is to describe an approach that is applicable for two-dimensional waves, with emphasis on the wave propagation after initiation. The restriction to two-dimensional flows not only offers the advantage of mathematical simplicity but is also anchored in practical considerations. Namely, many tsunami source regions—in particular, those for the major three tsunami events since modern record-keeping began in the twentieth century—are elliptical, with a major axis in excess of hundreds of kilometres corresponding to the active part of the fault. The fact that most of the tsunami energy is transmitted at right angles to the major axis allows us to regard in this context the propagation of the tsunami in the open sea as being two-dimensional.

## 2. The governing equations

The fluid domain, representing the sea near the fault area, is bounded above by the water's free surface *Y* =*d*+*F*(*X*,*T*), *d* being the average depth of the sea, and below by the rigid seabed *Y* =*H*(*X*,*T*), and is considered to be unbounded in the horizontal *X*-direction. At time *T*=0, the seabed is flat, given by *Y* =0. For some short time *T*_{0}>0, the seabed moves in the region *X*∈[−*L*,*L*], so that *H*(*X*,*T*)=0 for *T*≥*T*_{0} and *X*∉(−*L*,*L*). Our aim is to describe the deformation of the free surface *Y* =*d* that is induced by this sudden motion of the seabed. We assume that the water is inviscid and the resulting flow is a two-dimensional irrotational flow. In rectangular Cartesian coordinates (*X*,*Y*), the equations of motion are the incompressible Euler equations
2.1where *ρ* is the constant density, *g* is the constant acceleration of gravity, (*U*,*V*) is the velocity field of the flow and *P* is the pressure. As the effects of surface tension are negligible for wavelengths greater than a few centimetres [6,7], the major factor governing the wave motion is the balance between gravity and the inertia of the system. The corresponding kinematic boundary conditions,
2.2and
2.3ensure that both fluid boundaries are interfaces: particles on these boundaries are confined to them at all times. In addition, we also have the dynamic boundary condition
2.4where *P*_{atm} is the (constant) atmospheric pressure, which decouples the motion of the water from that of the air above it [8]. The above-mentioned system (2.1)–(2.4) is to be solved with the initial conditions expressing the fact that the water is at rest at time *T*=0, that is,
2.5the surface waves *F*(*X*,*T*) being generated by the prescribed motion *H*(*X*,*T*) of the bed. Because an incompressible inviscid flow that is irrotational initially remains so at later times, the initial conditions (2.5) ensure that throughout the wave propagation the vorticity of the flow vanishes, that is,
2.6For a discussion of the local well-posedness (existence, uniqueness and continuous dependence on initial data for sufficiently small time intervals) of the governing equations, we refer to Alvarez-Samaniego & Lannes [9] (see also [10,11]).

Let us now briefly describe the process of non-dimensionalizing and scaling the governing equations (2.1)–(2.6). In this process, we introduce the two fundamental parameters *ε* and *δ* that describe the regimes used to identify different general types of water wave propagation. To measure the change in pressure as a wave propagates at the water's surface, it is convenient to introduce the non-dimensional excess pressure, *p*, relative to the hydrostatic pressure distribution by
To generate the non-dimensional version of the governing equations, we introduce the average or typical wavelength *λ* of the wave, using as a scale for the wave speed and as a time scale. The change of variables
2.7coupled with the scaling
2.8where *a* is a typical—perhaps maximum—amplitude of the wave, introduces non-dimensional dependent and independent variables. Note that (2.8) should be interpreted as expressing the fact that the variations of the wave and of the seabed are of comparable size. Introducing the amplitude and shallowness parameters by
2.9which measure the relative size of the amplitude to average water depth and of the average water depth to wavelength, respectively, we obtain the non-dimensional version of the governing equations in the physical variables as
2.10The magnitudes of *ε* and *δ* can be used to identify different general types of wave problems. For example, the limits *δ*→0 and produce the shallow water and deep water regime, respectively, whereas *ε*→0 corresponds to the linear regime of waves of small amplitude. Before proceeding with the linearization approach, we would like to point out that the different non-dimensionalization of the horizontal and vertical fluid velocity components in (2.7) is consistent with the equation of mass conservation (and, equivalently, the existence of a stream function). Indeed, we may write *U*=*Ψ*_{Y} and *V* =−*Ψ*_{X} for a stream function *Ψ*, which leads to the transformation introduced in (2.7). The preservation of the irrotational character of the flow in addition to incompressibility in the non-dimensionalization process requires a regime in which the length scales in the horizontal and vertical directions are of comparable size (i.e. *δ*=1). The investigations in earlier studies [1–5] are performed within the framework of a linear theory in this type of regime. Nevertheless, some simple considerations show that it is possible to dispense with this limitation. To start with, we consider the linear regime *ε*→0 with *δ* fixed. A glance at the system (2.10) shows first that *v* and *p* are both of order *ε*, and consequently also *u*. We now perform the scaling
2.11avoiding the introduction of a new notation. The dimensionless, scaled governing equations are then
2.12We now linearize the system by taking the limit as *ε*→0. The linearized problem in non-dimensional scaled variables is
2.13While the fluid velocity in the system (2.13) does not originate from a velocity potential if *δ*≠1, there exists a stream function defined up to an additive function of time by
2.14We can re-formulate (2.13) as
2.15On *y*=1, we have *f*_{x}=*p*_{x}=−*ψ*_{yt} from the second and fourth relation in (2.15), so that *f*_{xt}=−*ψ*_{ytt} on *y*=1. In combination with the fifth relation in (2.15), this yields *ψ*_{xx}=*ψ*_{ytt} on *y*=1. We obtain the linear system
2.16To find *f* and *p*, note that, in terms of the solution *ψ* of (2.16), we have
2.17where *p* is determined by the second, third and fourth relation in (2.15).

Before proceeding with the analysis of the system (2.16) in §3, let us present some considerations related to its derivation. Switching from (2.13) back to the physical variables yields the linear system
2.18The first two equations in (2.18) ensure the existence of a velocity potential *Φ*(*X*,*Y*,*T*) that is harmonic in the spatial variables (*X*,*Y*) throughout the fluid domain, and such that *U*=*Φ*_{X}, *V* =*Φ*_{Y} throughout the domain 0<*Y* <*d*. We can therefore write (2.18) as
2.19The second and third equations in (2.19) are equivalent to the fact that the expression *Φ*_{T}+(1/*ρ*)(*P*−*P*_{atm})+*g*(*Y* −*d*) is, at any fixed instant *T*≥0, constant throughout the fluid domain. Absorbing an additive function of time in the definition of *Φ*, we define *Φ* unambiguously by requiring that
at every fixed *T*≥0. This, in combination with the fourth relation in (2.19), yields *Φ*_{T}+*gF*(*X*,*T*)=0 on *Y* =*d*. Taking into account the fifth relation in (2.19), we obtain *Φ*_{TT}+*gΦ*_{Y}=0 on *Y* =*d*. Thus
2.20This is precisely the linearization of the governing equations obtained by Dutykh & Dias [1] under the assumption that the horizontal and vertical length scales of the motion are comparable. The approach by Dutykh & Dias [1] was different, starting from a direct linearization of the governing equations, under the assumption that the horizontal and vertical length scales are comparable. For the sake of completeness, we briefly sketch this approach. One starts from the governing equations (2.1)–(2.6) in the physical variables. The first equation in (2.1) and (2.6) ensures the existence of a velocity potential *Φ*(*X*,*Y*,*T*) that is harmonic in the spatial variables (*X*,*Y*) throughout the fluid domain and such that *U*=*Φ*_{X}, *V* =*Φ*_{Y}. Then, the Euler equations in (2.1) simply mean that at each instant *T* the expression
is constant throughout the fluid (this is Bernoulli's law). In the absence of waves, that is, for *Φ*≡0 and for the hydrostatic pressure *P*=*P*_{atm}−*ρg*(*Y* −*d*), this expression is zero. Even in the presence of waves, absorbing the constant in the definition of *Φ* we may consider that the expression displayed above is always identically zero in the fluid. Evaluating the expression on the free surface *Y* =*d*+*F*(*X*,*T*), where (2.4) holds, yields
The linearization of the above-mentioned relation reads
2.21On the other hand, the linearization of (2.1)–(2.6), obtained by neglecting all nonlinear terms, yields
Because *U*=*Φ*_{X} and *V* =*Φ*_{Y}, we can write (2.21) and the above-mentioned relations as the system
Consequently,
and we recover the system (2.20).

Rather than working with the linear system (2.20) in the physical variables, we prefer to analyse the system (2.16) in non-dimensional scaled variables. This approach has the advantage that, once we obtain the formula for the solution for a fixed *δ*, because we are in non-dimensional variables, it is permissible to neglect small terms and we can consider the limits *δ*→0 (shallow water waves) and (deep water waves). Another advantage of the formulation (2.13) is that it can be easily adapted to deal with seabed disturbances that do not act on still water but on a water flow that might even be rotational (such flows with vorticity describe non-uniform currents—see [12,13] for aspects of wave–current interactions related to tsunami). The reason is that, in contrast to (2.20), a velocity potential is not required for (2.13) and stream functions are known to exist for two-dimensional water flows with vorticity (see [14,15]).

## 3. The general solution formula for linear waves

We solve the problem (2.16) using the space–time Fourier transform. More precisely, we will denote in the following
3.1This enables us to define on the Schwartz class (functions that are smooth and for which all derivatives decay faster than polynomials at infinity (see [16])) the Fourier multiplier *m*(*D*) by
where *D*=−*i*∂_{x}. With these notations,
3.2The general solution of the differential equation is of the form
and *A*,*B* are determined from the boundary conditions. Because, by (2.17), we have
a straightforward calculation yields
3.3Because the function is even on and an increasing diffeomorphism of , for every fixed *ω*>0, there are precisely two roots *ξ*_{±}(*ω*) with *ξ*_{−}(*ω*)=−*ξ*_{+}(*ω*) of the equation . In the limit *ω*↓0, the two roots coalesce to *ξ*_{±}(0)=0, whereas *ξ*_{+}(*ω*)*ω*^{2}*δ* for because for . It is convenient to apply to (3.3) the inverse Fourier transform, obtaining that
3.4Let us now assume that the up- and downward motions of the sea bottom compensate at each instant,
3.5This will cancel a singularity in the next formula. The considerations in appendix A (see (A7) and (A8)) yield the representation
3.6of the solution, where .

## 4. Asymptotic behaviour of linear waves generated by a moving bed

### (a) The shallow water regime

As , equation (3.4) becomes the inhomogeneous linear wave equation
The initial conditions are *f*(*x*,0)=*f*_{t}(*x*,0)=0 owing to (2.5) and (2.17), so that Duhamel's formula (see [17]) yields
4.1We will now see that this model gives surprisingly accurate predictions for the propagation of some tsunamis. Let *h*(*x*,*t*)=*a*(*t*) *b*(*x*) with , such that *a*(*t*)=0 for *t*≤0, *a*(*t*)=1 for *t*>*t*_{0} (with *t*_{0}>0 representing the duration of the earthquake in the case study of the generation of tsunami waves by submarine earthquakes), and with *b*(*x*)=0 for *x*∉(−*L*,*L*) modelling a localized tsunami source. In this case, (4.1) takes the simpler form
4.2In the limiting case *t*_{0}↓0, *a*(*t*) becomes the Heaviside step function (modelling an instantaneous upward thrust of the seabed near to the earthquake's epicentre) and (4.2) simplifies to
4.3In the special case (4.3), some simple but insightful conclusions can be drawn as follows:

— at each instant

*t*>0 after initiation, the generated wave is localized as*f*(*x*,*t*)=0 for |*x*|≥*L*+*t*;— the surface wave consists of two travelling waves, one moving to the right and the other moving to the left; and

— each of the above travelling waves moves with non-dimensionalized speed 1 (corresponding to the speed in the physical variables) and with a shape that remains unchanged and is precisely that of the bed deformation at half-scale—in particular, the amplitude is half that of the bed deformation.

These conclusions are realistic. Indeed, let us discuss the two largest tsunamis for which records are available—the December 2004 and the May 1960 tsunamis. As for the March 2011 tsunami, the wave propagation in the open ocean is of secondary importance. The reason for this is that, for this particular event, the tsunami waves were generated about 70 km offshore Japan, and the waves that hit Japan propagated in water of rapidly diminishing depth. As for the tsunami waves towards Hawaii, their propagation also occurred in an ocean region with considerable changes in the ocean bathymetry (see the discussion in Constantin [18]).

The December 2004 tsunami was triggered by a submarine earthquake at the interface of the India and Burma plates, off Sumatra. The tsunami initiation region was roughly elliptical, about 100 km wide and more than 1000 km long, with the disturbance subsequently moving eastward towards Sumatra and Thailand, and westward towards Sri Lanka and India (see [19]). The generating earthquake raised the ocean floor by up to 2 m over a distance of about 100 km to the west of the epicentre and lowered it to the east (over 900 km). The Indian Ocean/Bay of Bengal and the Andaman Basin through which the tsunami propagated are relatively flat, with average depths *d*=4 km and *d*=1 km, respectively. As an indication of the applicability of the above simple considerations, we notice that, for India and Sri Lanka, the first tsunami wave reaching the shore was a wave of elevation, whereas in southern Thailand a wave of depression signalled the arrival of the tsunami. This is consistent with our predictions based on (4.3) because the bed deformation presents a depression eastward and elevation westward of the earthquake's epicentre. Also, accurate measurements, provided by a radar altimeter on board a satellite travelling along a track traversing the Indian Ocean/Bay of Bengal, indicate that out in the open sea the tsunami amplitude was about 1 m (see [13]). This is consistent with the halving of the amplitude encompassed in (4.3). Furthermore, the above-mentioned theory predicts a wave speed of approximately 712 km h^{−1}, corresponding to with *g*=9.8 m s^{−2} and *d*=4 km. This would mean that in 2 h 12 min the tsunami waves would cover a distance of about 1566 km, and this is the real time the tsunami needed for the 1550 km across the Indian Ocean/Bay of Bengal, between the epicentre of the earthquake and the first affected coast of Sri Lanka (see [20]). For the Andaman Basin, a predicted speed of 356 km h^{−1} would correspond to with *d*=1 km. This is also quite accurate because the Thai resort at Hat Ray Leh, about 700 km from the epicentre of the earthquake, was hit by the tsunami 2 h after initiation (see [13]). While the field data for the May 1960 tsunami are scarcer (see, however, the reference lists of Constantin & Henry [21] and Stuhlmeier [22]), it is known that it was caused by the largest earthquake ever recorded. Several earthquakes in fast succession occurred along 1000 km of fault parallel to the Chilean coastline, with the epicentre within 200 km off the coast of Central Chile, and a submarine lift of 1 m and subsidence of 1.6 m ensued over a stretch of 300 km (see [20]). Tsunami waves propagated in the northwest direction across the Pacific Ocean, reaching Hawaii, at about 10 500 km from the epicentre of the earthquake, after 14 h 48 min. The ocean floor of the Central Pacific Basin between Chile and Hawaii is relatively uniform, with a mean depth of approximately *d*=4 km. The above-mentioned simple theory would, therefore, predict a wave speed km h^{−1}, so that the waves would travel 10 535 km in 14 h 48 min, which is about right.

### (b) The regime of finite, non-vanishing *δ*

The asymptotic formula (A13) for the solution to equation (3.4) is derived in appendix A. It gives
where
For tsunamis propagating in open sea, typical values of the physical parameters introduced above are (for instance, see Constantin [20])
This leads to *ε*=0.00025 and *λ*=0.02, which point towards the linear wave approximation. The stationary-phase analysis conducted above is justified if
4.4Using the discussion in appendix A, can be written as *F*(*δξ*), for an odd function *F*, decreasing on , with *F*(0+)=1 and . It makes sense to restrict away from 0: very low frequencies travel very slowly, and correspond thus to waves which will barely reach the coast. With this restriction (say ), the stationary-phase point *ξ*>0 where is of order 1/*δ*. Now, for *δξ*=*O*(1), (A12) yields so that from (4.4) we get the criterion
Coming back to physical variables, this corresponds to
Thus, with the choice of the above-mentioned parameters, stationary-phase analysis does not seem to be justified. However, under different physical conditions, stationary-phase analysis could be the right choice.

## Acknowledgements

The support by the ERC advanced grant ‘Nonlinear studies of water flows with vorticity’ (NWFV) is gratefully acknowledged. It is a pleasure to thank the referees for very helpful comments.

## Appendix A. Asymptotic formulae

(*a*) *A general approach*

We want to derive an asymptotic formula, as , for the solution *f*(*x*,*t*) of
A1where *τ* is a given dispersion relation, such that *τ*(*ξ*)>0 for any *ξ*≠0, and *θ* is a given forcing term. We assume that *θ* and *f* are both of Schwartz class , the limit in (A1) being also interpreted in the Schwartz class. If solves (A1), then the function solves . Filtering *v* by the group gives that solves . This last equation can be solved by integrating in time; coming back to *f*, this means
A2Spelling out the Fourier multipliers in the above formula, and writing the formula in order to highlight the oscillatory nature of the integral, we see that, with ,
The stationary-phase principle gives
for *t* large, where *ξ*_{0} is the function of such that , and *σ* is the sign of . Performing the integration over in the above-mentioned integral gives the asymptotics: for *t* large,
A3The above-mentioned derivation has been rigorous, except for two points: we assumed implicitly the uniqueness of *ξ*_{0}; and we divided by , which is problematic where this function vanishes. We will address these two issues in the subsequent paragraph, where we specialize the above-mentioned analysis to the case of linear water waves.

(*b*) *The case of linear water waves*

We apply here the analysis of the previous paragraph to the equation A4where is a test function (smooth and with compact support). The setting of the previous paragraph corresponds to the notations A5and A6It is easily seen that

— is an even function,

*τ*(−*ξ*)=*τ*(*ξ*);—

*τ*(*ξ*)>0 for*ξ*≠0;— ; and

—

*τ*has linear growth at infinity.

Because the first property ensures that *τ*(*D*) *φ* is real-valued if , (A2) yields
so that
A7Let us now assume that is such that (3.5) holds—this hypothesis expresses that at each instant the upward motion of the seabed is compensated by the downward motion. This hypothesis cancels the singularity at *ξ*=0 which is present in (A7). Indeed, as (3.5) ensures for all and
A8we get for all . To make the cancellation of the singularity more obvious, we can rewrite (A7) as
A9As the Fourier transform acts on the Schwartz class as an isomorphism, the absence of singularities on the right of (A9) ensures . For the evaluation of the asymptotic behaviour of the solution (A9) by means of the stationary-phase principle, let us write (A9) in the equivalent form
A10where . A computation gives
A11and
A12Consequently, is strictly decreasing on from the asymptotic value 1 towards the asymptotic value 0, and on from the asymptotic value 0 towards the asymptotic value −1. Thus, stationary points that possibly contribute towards the leading order exist only for the range , and in this range each fixed produces precisely one stationary point . As
the formula (A3) becomes
A13(notice that the singular factor present in (A3) has been cancelled). For *ξ*>0 fixed, by expanding in (A11) in power series (in terms of *δξ*), we see that is a smooth () and even function of *δ*. As , we must have for some function *α*(*ξ*). To determine *α*(*ξ*), notice that 2*α*(*ξ*) is the coefficient of *δ*^{2} in the expansion of with respect to the parameter *δ*. As
expanding in power series with respect to *δξ* yields . Therefore, . For we deduce that the stationary-phase point *ξ*_{0}>0 where is of order 1/*δ*. Now, for *δξ*_{0}=*O*(1), note that (A12) yields .

## Footnotes

One contribution of 13 to a Theme Issue ‘Nonlinear water waves’.

- This journal is © 2012 The Royal Society