On the open sea propagation of two-dimensional rotational water waves generated by a moving bed

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:

• Water waves
• free boundary
• vorticity

AMS Classifications:

• 45G99
• 58J32
• 76B03

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 [1]). 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 [2]); 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 [3] 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 [1,3,4]). 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 [5–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 [1,3,8–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 [4,11] 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 [12,13]) 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 [11]: we consider two-dimensional surface waves with the unbounded X-direction corresponding to the direction of wave propagation. The water's free surface is $Y=d+F(X,T)$, where d is the average depth of the sea, the seabed being $Y=H(X,T)$. At time T = 0, we assume that $F(X,0)=0$ and $H(X,0)=0$. Then the seabed moves for a short time in some compact region $X\in [-L,L]$ and is flat after that again. We are ultimately interested in the deformation of the free surface, i.e. in $F(X,T)$. To do this, we study a partial differential equation that describes its evolution in time. Consider $(U,V)$ the velocity field of the flow in rectangular Cartesian coordinates $(X,Y)$. 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. 12]. 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 [13]. We now present the governing equations for water flows [1]. We have the incompressible Euler equations (1) $\begin{array}{r}\begin{array}{rl}{U}_X+V_Y& =0,\\ \rho (U_T+U{U}_X+V{U}_Y)& =-P_X,\\ \rho (V_T+U{V}_X+V{V}_Y)& =-P_Y-\rho g,\end{array}\}\end{array}$(1) in $H(X,T)<Y<d+F(X,T)$, 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 [14,15], 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) $\begin{array}{r}V=F_T+U{F}_X\text{on the free surface}Y=d+F(X,T)\end{array}$(2) and (3) $\begin{array}{r}V=H_T+U{H}_X\text{on the rigid bed}Y=H(X,T),\end{array}$(3) which ensure that particles on these boundaries are confined to them at all times. We also have the dynamic boundary condition (4) $\begin{array}{r}P-P_{\mathrm{a}\mathrm{t}\mathrm{m}}=0\text{on the free surface}Y=d+F(X,T),\end{array}$(4) where $P_{\mathrm{a}\mathrm{t}\mathrm{m}}$ is the constant atmospheric pressure. This decouples the motion of the water from the motion of the air above it [16]. The system (1(1) $\begin{array}{r}\begin{array}{rl}{U}_X+V_Y& =0,\\ \rho (U_T+U{U}_X+V{U}_Y)& =-P_X,\\ \rho (V_T+U{V}_X+V{V}_Y)& =-P_Y-\rho g,\end{array}\}\end{array}$(1) )–(4(4) $\begin{array}{r}P-P_{\mathrm{a}\mathrm{t}\mathrm{m}}=0\text{on the free surface}Y=d+F(X,T),\end{array}$(4) ) is to be solved with the following initial conditions (5) $\begin{array}{r}F(X,0)=0,U(X,Y,0)=AY+B,V(X,Y,0)=0.\end{array}$(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 [11], where it was assumed that the water is completely at rest at T = 0, a case which corresponds to $A=B=0$. 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 [1]). By the initial condition (5(5) $\begin{array}{r}F(X,0)=0,U(X,Y,0)=AY+B,V(X,Y,0)=0.\end{array}$(5) ) this implies (6) $\begin{array}{r}{U}_Y-V_X=A\end{array}$(6) for $H(X,T)<Y<d+F(X,T)$. Figure 1 illustrates our set-up in the case $A,B>0$.

Figure 1. Pure current flow (in the absence of waves) at the initial time T = 0.

The first step in the analysis of the system (1(1) $\begin{array}{r}\begin{array}{rl}{U}_X+V_Y& =0,\\ \rho (U_T+U{U}_X+V{U}_Y)& =-P_X,\\ \rho (V_T+U{V}_X+V{V}_Y)& =-P_Y-\rho g,\end{array}\}\end{array}$(1) )–(6(6) $\begin{array}{r}{U}_Y-V_X=A\end{array}$(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) $\begin{array}{r}P=P_{\mathrm{a}\mathrm{t}\mathrm{m}}+\rho g(d-Y)+\rho gdp.\end{array}$(7) To obtain meaningful scales, we introduce the average or typical wavelength λ of the wave, we use $\sqrt{gd}$ as a scale for the wave speed and $\lambda /\sqrt{gd}$ as a time scale. We use the change of variables (8) $\begin{array}{r}X=\lambda x,Y=dy,T=\frac{\lambda }{\sqrt{gd}}t,U=\sqrt{gd}u,V=\frac{d\sqrt{gd}}{\lambda }v\end{array}$(8) with (9) $\begin{array}{r}F=af,H=ah,\end{array}$(9) where a is a typical, perhaps maximal, amplitude of the wave. Note that (9(9) $\begin{array}{r}F=af,H=ah,\end{array}$(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) $\begin{array}{r}ϵ=\frac{a}{d},\delta =\frac{d}{\lambda },\end{array}$(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) $\begin{array}{rl}& \begin{array}{rl}& \{\begin{array}{l}{u}_x+v_y=0\\ u_t+u{u}_x+v{u}_y=-p_x\\ {\delta }^2(v_t+u{v}_x+v{v}_y)=-p_y\\ u_y-{\delta }^2{v}_x=\sqrt{\frac{d}{g}}A\end{array}\mathrm{i}\mathrm{n}ϵh(x,t)<y<1+ϵf(x,t)\\ & \begin{array}{ll}p=ϵf& \mathrm{o}\mathrm{n}y=1+ϵf(x,t)\\ v=ϵ(f_t+u{f}_x)& \mathrm{o}\mathrm{n}y=1+ϵf(x,t)\\ v=ϵ(h_t+u{h}_x)& \mathrm{o}\mathrm{n}y=ϵh(x,t)\end{array}\\ & f(x,0)=0,u(x,y,0)=\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}},v(x,y,0)=0\end{array}\}.\end{array}$(11) The magnitudes of ε and δ correspond to the different general types of water wave problem: The limits $\delta \to 0$ and $\delta \to \mathrm{\infty }$ produce the shallow water and deep water regime, respectively. On the other hand, $ϵ\to 0$ corresponds to regime of waves of small amplitude (see the discussion in [16,17]).

We want to linearise the problem and in this case we will want to let $ϵ\to 0$ and keep δ fixed. But the system (11(11) $\begin{array}{rl}& \begin{array}{rl}& \{\begin{array}{l}{u}_x+v_y=0\\ u_t+u{u}_x+v{u}_y=-p_x\\ {\delta }^2(v_t+u{v}_x+v{v}_y)=-p_y\\ u_y-{\delta }^2{v}_x=\sqrt{\frac{d}{g}}A\end{array}\mathrm{i}\mathrm{n}ϵh(x,t)<y<1+ϵf(x,t)\\ & \begin{array}{ll}p=ϵf& \mathrm{o}\mathrm{n}y=1+ϵf(x,t)\\ v=ϵ(f_t+u{f}_x)& \mathrm{o}\mathrm{n}y=1+ϵf(x,t)\\ v=ϵ(h_t+u{h}_x)& \mathrm{o}\mathrm{n}y=ϵh(x,t)\end{array}\\ & f(x,0)=0,u(x,y,0)=\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}},v(x,y,0)=0\end{array}\}.\end{array}$(11) ) shows that v and p are of order ε and hence also u. Hence we use the following scaling: (12) $\begin{array}{r}u\mapsto \sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}+ϵu,v\mapsto ϵv,p\mapsto ϵp,\end{array}$(12) where we avoid a new notation. We then obtain the system (13) $\begin{array}{r}\begin{array}{rl}{u}_x+v_y& =0\\ u_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}+ϵu)u_x+v(\sqrt{\frac{d}{g}}A+ϵu_y)& =-p_x\\ {\delta }^2(v_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}+ϵu)v_x+ϵv{v}_y)& =-p_y\\ u_y-{\delta }^2{v}_x& =0\end{array}\}\end{array}$(13) in $ϵh(x,t)<y<1+ϵf(x,t)$ and (14) $\begin{array}{r}\begin{array}{rl}& p=f\mathrm{o}\mathrm{n}y=1+ϵf(x,t)\\ & v=f_t+\left(\sqrt{\frac{d}{g}}A\right(1+ϵf(x,t))+\frac{B}{\sqrt{gd}}+ϵu)f_x\mathrm{o}\mathrm{n}y=1+ϵf(x,t)\\ & v=h_t+\left(\sqrt{\frac{d}{g}}Aϵh\right(x,t)+\frac{B}{\sqrt{gd}}+ϵu)h_x\mathrm{o}\mathrm{n}y=ϵh(x,t)\\ & f(x,0)=0,u(x,y,0)=0,v(x,y,0)=0\end{array}\}.\end{array}$(14) We now linearise this system by taking $ϵ\to 0$ and obtain (15) $\begin{array}{r}\begin{array}{rl}& \{\begin{array}{l}{u}_x+v_y=0\\ u_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}})u_x+\sqrt{\frac{d}{g}}Av=-p_x\\ {\delta }^2(v_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}})v_x)=-p_y\\ u_y-{\delta }^2{v}_x=0\end{array}\mathrm{i}\mathrm{n}0<y<1\\ & \begin{array}{ll}p=f& \mathrm{o}\mathrm{n}y=1\\ v=f_t+(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_x& \mathrm{o}\mathrm{n}y=1\\ v=h_t+\frac{B}{\sqrt{gd}}{h}_x& \mathrm{o}\mathrm{n}y=0\end{array}\\ & f(x,0)=0,u(x,y,0)=0,v(x,y,0)=0\end{array}\}.\end{array}$(15) By the first equation in (15(15) $\begin{array}{r}\begin{array}{rl}& \{\begin{array}{l}{u}_x+v_y=0\\ u_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}})u_x+\sqrt{\frac{d}{g}}Av=-p_x\\ {\delta }^2(v_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}})v_x)=-p_y\\ u_y-{\delta }^2{v}_x=0\end{array}\mathrm{i}\mathrm{n}0<y<1\\ & \begin{array}{ll}p=f& \mathrm{o}\mathrm{n}y=1\\ v=f_t+(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_x& \mathrm{o}\mathrm{n}y=1\\ v=h_t+\frac{B}{\sqrt{gd}}{h}_x& \mathrm{o}\mathrm{n}y=0\end{array}\\ & f(x,0)=0,u(x,y,0)=0,v(x,y,0)=0\end{array}\}.\end{array}$(15) ), there exists a stream function ψ such that (16) $\begin{array}{r}u={\psi }_y,v=-{\psi }_x\mathrm{i}\mathrm{n}0<y<1.\end{array}$(16) If we plug this in (15(15) $\begin{array}{r}\begin{array}{rl}& \{\begin{array}{l}{u}_x+v_y=0\\ u_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}})u_x+\sqrt{\frac{d}{g}}Av=-p_x\\ {\delta }^2(v_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}})v_x)=-p_y\\ u_y-{\delta }^2{v}_x=0\end{array}\mathrm{i}\mathrm{n}0<y<1\\ & \begin{array}{ll}p=f& \mathrm{o}\mathrm{n}y=1\\ v=f_t+(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_x& \mathrm{o}\mathrm{n}y=1\\ v=h_t+\frac{B}{\sqrt{gd}}{h}_x& \mathrm{o}\mathrm{n}y=0\end{array}\\ & f(x,0)=0,u(x,y,0)=0,v(x,y,0)=0\end{array}\}.\end{array}$(15) ), we get (17) $\begin{array}{r}\begin{array}{rl}& \{\begin{array}{l}{\psi }_{yy}+{\delta }^2{\psi }_{xx}=0\\ {\psi }_{yt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xy}-\sqrt{\frac{d}{g}}A{\psi }_x=-p_x\\ {\delta }^2({\psi }_{xt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xx})=p_y\end{array}\mathrm{i}\mathrm{n}0<y<1\\ & \begin{array}{ll}p=f& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-f_t-(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_x& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-h_t-\frac{B}{\sqrt{gd}}{h}_x& \mathrm{o}\mathrm{n}y=0\end{array}\\ & f(x,0)=0,{\psi }_y(x,y,0)=0,{\psi }_x(x,y,0)=0\end{array}\}.\end{array}$(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 (17(17) $\begin{array}{r}\begin{array}{rl}& \{\begin{array}{l}{\psi }_{yy}+{\delta }^2{\psi }_{xx}=0\\ {\psi }_{yt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xy}-\sqrt{\frac{d}{g}}A{\psi }_x=-p_x\\ {\delta }^2({\psi }_{xt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xx})=p_y\end{array}\mathrm{i}\mathrm{n}0<y<1\\ & \begin{array}{ll}p=f& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-f_t-(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_x& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-h_t-\frac{B}{\sqrt{gd}}{h}_x& \mathrm{o}\mathrm{n}y=0\end{array}\\ & f(x,0)=0,{\psi }_y(x,y,0)=0,{\psi }_x(x,y,0)=0\end{array}\}.\end{array}$(17) ) with respect to x and t, respectively, to obtain (18) $\begin{array}{r}{\psi }_{xyt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xxy}-\sqrt{\frac{d}{g}}A{\psi }_{xx}=-p_{xx}\end{array}$(18) and (19) $\begin{array}{r}{\psi }_{ytt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xyt}-\sqrt{\frac{d}{g}}A{\psi }_{xt}=-p_{xt}.\end{array}$(19) Furthermore, we can differentiate the fifth equation in (17(17) $\begin{array}{r}\begin{array}{rl}& \{\begin{array}{l}{\psi }_{yy}+{\delta }^2{\psi }_{xx}=0\\ {\psi }_{yt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xy}-\sqrt{\frac{d}{g}}A{\psi }_x=-p_x\\ {\delta }^2({\psi }_{xt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xx})=p_y\end{array}\mathrm{i}\mathrm{n}0<y<1\\ & \begin{array}{ll}p=f& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-f_t-(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_x& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-h_t-\frac{B}{\sqrt{gd}}{h}_x& \mathrm{o}\mathrm{n}y=0\end{array}\\ & f(x,0)=0,{\psi }_y(x,y,0)=0,{\psi }_x(x,y,0)=0\end{array}\}.\end{array}$(17) ) with respect to x to obtain (20) $\begin{array}{r}{\psi }_{xx}=-f_{xt}-(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_{xx}\mathrm{o}\mathrm{n}y=1.\end{array}$(20) Observe that the fourth equation in (17(17) $\begin{array}{r}\begin{array}{rl}& \{\begin{array}{l}{\psi }_{yy}+{\delta }^2{\psi }_{xx}=0\\ {\psi }_{yt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xy}-\sqrt{\frac{d}{g}}A{\psi }_x=-p_x\\ {\delta }^2({\psi }_{xt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xx})=p_y\end{array}\mathrm{i}\mathrm{n}0<y<1\\ & \begin{array}{ll}p=f& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-f_t-(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_x& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-h_t-\frac{B}{\sqrt{gd}}{h}_x& \mathrm{o}\mathrm{n}y=0\end{array}\\ & f(x,0)=0,{\psi }_y(x,y,0)=0,{\psi }_x(x,y,0)=0\end{array}\}.\end{array}$(17) ) implies that we can exchange the derivatives of f in (20(20) $\begin{array}{r}{\psi }_{xx}=-f_{xt}-(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_{xx}\mathrm{o}\mathrm{n}y=1.\end{array}$(20) ) by the corresponding derivatives of p. But then we can insert Equations (18(18) $\begin{array}{r}{\psi }_{xyt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xxy}-\sqrt{\frac{d}{g}}A{\psi }_{xx}=-p_{xx}\end{array}$(18) ) and (19(19) $\begin{array}{r}{\psi }_{ytt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xyt}-\sqrt{\frac{d}{g}}A{\psi }_{xt}=-p_{xt}.\end{array}$(19) ) in (20(20) $\begin{array}{r}{\psi }_{xx}=-f_{xt}-(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_{xx}\mathrm{o}\mathrm{n}y=1.\end{array}$(20) ) to obtain $\begin{array}{rl}{\psi }_{xx}& ={\psi }_{ytt}+2(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}}){\psi }_{xyt}+{(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})}^2{\psi }_{xxy}\\ & -\sqrt{\frac{d}{g}}A(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}}){\psi }_{xx}-\sqrt{\frac{d}{g}}A{\psi }_{xt}\end{array}$on y = 1. We introduce the constant (21) $\begin{array}{r}C=\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}},\end{array}$(21) which is the non-dimensional horizontal velocity of the flow at the surface. This allows us to rewrite the previous equation as $\begin{array}{r}(1+\sqrt{\frac{d}{g}}AC){\psi }_{xx}={\psi }_{ytt}+2C{\psi }_{xyt}+C^2{\psi }_{xxy}-\sqrt{\frac{d}{g}}A{\psi }_{xt}\end{array}$on y = 1. By observing the quadratic expression on the right-hand side and introducing the differential operator (22) $\begin{array}{r}S=C{\mathrm{\partial }}_x+{\mathrm{\partial }}_t,\end{array}$(22) we can rewrite this again as $\begin{array}{r}(1+\sqrt{\frac{d}{g}}AC){\psi }_{xx}=S^2{\psi }_y-\sqrt{\frac{d}{g}}A{\psi }_{xt}.\end{array}$So we obtain (23) $\begin{array}{r}\begin{array}{rl}{\psi }_{yy}+{\delta }^2{\psi }_{xx}=0,& \text{in }0<y<1\\ (1+\sqrt{\frac{d}{g}}AC){\psi }_{xx}=S^2{\psi }_y-\sqrt{\frac{d}{g}}A{\psi }_{xt},& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-h_t-\frac{B}{\sqrt{gd}}{h}_x,& \mathrm{o}\mathrm{n}y=0\\ {\psi }_x(x,y,0)=0,& {\psi }_y(x,y,0)=0,\end{array}\}\end{array}$(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 (17(17) $\begin{array}{r}\begin{array}{rl}& \{\begin{array}{l}{\psi }_{yy}+{\delta }^2{\psi }_{xx}=0\\ {\psi }_{yt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xy}-\sqrt{\frac{d}{g}}A{\psi }_x=-p_x\\ {\delta }^2({\psi }_{xt}+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}}){\psi }_{xx})=p_y\end{array}\mathrm{i}\mathrm{n}0<y<1\\ & \begin{array}{ll}p=f& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-f_t-(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_x& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-h_t-\frac{B}{\sqrt{gd}}{h}_x& \mathrm{o}\mathrm{n}y=0\end{array}\\ & f(x,0)=0,{\psi }_y(x,y,0)=0,{\psi }_x(x,y,0)=0\end{array}\}.\end{array}$(17) ) that $\begin{array}{r}{f}_x=p_x=-{\psi }_{yt}+C{\psi }_{xy}-\sqrt{\frac{d}{g}}A{\psi }_x\end{array}$on y = 1 and hence (24) $\begin{array}{r}{f}_t=-{\psi }_x-C{f}_x=C{\psi }_{yt}-C^2{\psi }_{xy}+(\sqrt{\frac{d}{g}}AC-1){\psi }_x\end{array}$(24) on y = 1. Together with the initial condition $\begin{array}{r}f(x,0)=0\end{array}$we can recover f by $\begin{array}{r}f(x,t)={\int }_0^t{f}_t(x,\tau )\mathrm{d}\tau ,\end{array}$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 $\begin{array}{r}\hat{\phi }(\xi ,y,t)=\frac{1}{\sqrt{2\pi }}{\int }_{\mathbb{R}}\phi (x,y,t){\mathrm{e}}^{-ix\xi }\mathrm{d}x\end{array}$and $\begin{array}{r}\tilde{\phi }(\xi ,y,\omega )=\frac{1}{2\pi }{\int }_{\mathbb{R}}{\int }_{\mathbb{R}}\phi (x,y,t){\mathrm{e}}^{-i(x\xi +t\omega )}\mathrm{d}x\mathrm{d}t.\end{array}$We also define Fourier multipliers: Let $m:\mathbb{R}\to \mathbb{C}$ be some function. We define $m\left(D\right)$ by $\begin{array}{r}\left(m\right(D\left)\phi \right)\left(x\right)=\frac{1}{\sqrt{2\pi }}{\int }_{\mathbb{R}}m\left(\xi \right)\hat{\phi }\left(\xi \right){\mathrm{e}}^{ix\xi }\mathrm{d}\xi ,\end{array}$or equivalently $\begin{array}{r}\left(\widehat{m\left(D\right)\phi }\right)\left(\xi \right)=m\left(\xi \right)\hat{\phi }\left(\xi \right),\end{array}$where $D=-i{\mathrm{\partial }}_x$, whenever this is defined. Note that $m\left(D\right)$ maps real-valued functions to real-valued functions whenever the condition (25) $\begin{array}{r}m(-\xi )=\overline{m\left(\xi \right)}\end{array}$(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 (23(23) $\begin{array}{r}\begin{array}{rl}{\psi }_{yy}+{\delta }^2{\psi }_{xx}=0,& \text{in }0<y<1\\ (1+\sqrt{\frac{d}{g}}AC){\psi }_{xx}=S^2{\psi }_y-\sqrt{\frac{d}{g}}A{\psi }_{xt},& \mathrm{o}\mathrm{n}y=1\\ {\psi }_x=-h_t-\frac{B}{\sqrt{gd}}{h}_x,& \mathrm{o}\mathrm{n}y=0\\ {\psi }_x(x,y,0)=0,& {\psi }_y(x,y,0)=0,\end{array}\}\end{array}$(23) ) to obtain (26) $\begin{array}{r}\begin{array}{rl}{\tilde{\psi }}_{yy}-{\delta }^2{\xi }^2\tilde{\psi }=0,& \mathrm{i}\mathrm{n}0<y<1\\ (1+\sqrt{\frac{d}{g}}AC){\xi }^2\tilde{\psi }=Q^2{\tilde{\psi }}_y-\sqrt{\frac{d}{g}}A\xi \omega \tilde{\psi },& \mathrm{o}\mathrm{n}y=1\\ \xi \tilde{\psi }=-\omega \tilde{h}-\frac{B}{\sqrt{gd}}\xi \tilde{h},& \mathrm{o}\mathrm{n}y=0\end{array}\},\end{array}$(26) where (27) $\begin{array}{r}Q(\xi ,\omega )=C\xi +\omega \end{array}$(27) with the relation (28) $\begin{array}{r}\widetilde{S\varphi }=iQ\text{∼}\varphi .\end{array}$(28) Since (26(26) $\begin{array}{r}\begin{array}{rl}{\tilde{\psi }}_{yy}-{\delta }^2{\xi }^2\tilde{\psi }=0,& \mathrm{i}\mathrm{n}0<y<1\\ (1+\sqrt{\frac{d}{g}}AC){\xi }^2\tilde{\psi }=Q^2{\tilde{\psi }}_y-\sqrt{\frac{d}{g}}A\xi \omega \tilde{\psi },& \mathrm{o}\mathrm{n}y=1\\ \xi \tilde{\psi }=-\omega \tilde{h}-\frac{B}{\sqrt{gd}}\xi \tilde{h},& \mathrm{o}\mathrm{n}y=0\end{array}\},\end{array}$(26) ) is just a linear ODE of second order in y for fixed ξ and ω, we can, somewhat explicitly, write down the solution: (29) $\begin{array}{r}\tilde{\psi }(\xi ,y,\omega )=D_1(\xi ,\omega ){\mathrm{e}}^{\delta \xi y}+D_2(\xi ,\omega ){\mathrm{e}}^{-\delta \xi y},\end{array}$(29) where $D_1,D_2$ are determined by the boundary conditions. We can then reconstruct $\tilde{f}$ by applying the Fourier transform on the first equation in (24(24) $\begin{array}{r}{f}_t=-{\psi }_x-C{f}_x=C{\psi }_{yt}-C^2{\psi }_{xy}+(\sqrt{\frac{d}{g}}AC-1){\psi }_x\end{array}$(24) ) to obtain $\begin{array}{r}\omega \tilde{f}(\xi ,\omega )=-\xi \tilde{\psi }(\xi ,1,\omega )-C\xi \tilde{f}(\xi ,\omega )\end{array}$and hence (30) $\begin{array}{r}\tilde{f}(\xi ,\omega )=\frac{-\xi }{C\xi +\omega }\tilde{\psi }(\xi ,1,\omega )=\frac{-\xi }{Q(\xi ,\omega )}\tilde{\psi }(\xi ,1,\omega ).\end{array}$(30) We ignore the potential singularity for the time being. We know that the functions $D_1$ and $D_2$ are uniquely determined. Once we have formulae for them, we get a formula for $\tilde{\psi }$ by inserting in (29(29) $\begin{array}{r}\tilde{\psi }(\xi ,y,\omega )=D_1(\xi ,\omega ){\mathrm{e}}^{\delta \xi y}+D_2(\xi ,\omega ){\mathrm{e}}^{-\delta \xi y},\end{array}$(29) ). Then we can insert this formula in (30(30) $\begin{array}{r}\tilde{f}(\xi ,\omega )=\frac{-\xi }{C\xi +\omega }\tilde{\psi }(\xi ,1,\omega )=\frac{-\xi }{Q(\xi ,\omega )}\tilde{\psi }(\xi ,1,\omega ).\end{array}$(30) ) to obtain an expression for $\tilde{f}$ depending on $\tilde{h}$.

We insert (29(29) $\begin{array}{r}\tilde{\psi }(\xi ,y,\omega )=D_1(\xi ,\omega ){\mathrm{e}}^{\delta \xi y}+D_2(\xi ,\omega ){\mathrm{e}}^{-\delta \xi y},\end{array}$(29) ) in (26(26) $\begin{array}{r}\begin{array}{rl}{\tilde{\psi }}_{yy}-{\delta }^2{\xi }^2\tilde{\psi }=0,& \mathrm{i}\mathrm{n}0<y<1\\ (1+\sqrt{\frac{d}{g}}AC){\xi }^2\tilde{\psi }=Q^2{\tilde{\psi }}_y-\sqrt{\frac{d}{g}}A\xi \omega \tilde{\psi },& \mathrm{o}\mathrm{n}y=1\\ \xi \tilde{\psi }=-\omega \tilde{h}-\frac{B}{\sqrt{gd}}\xi \tilde{h},& \mathrm{o}\mathrm{n}y=0\end{array}\},\end{array}$(26) ) to obtain $\begin{array}{r}\xi \left(D_1\right(\xi ,\omega )+D_2(\xi ,\omega \left)\right)=-(\omega +\frac{B}{\sqrt{gd}}\xi )\text{∼}h(\xi ,\omega )\end{array}$and $\begin{array}{rl}& (\xi +\sqrt{\frac{d}{g}}AQ(\xi ,\omega \left)\right)\left(D_1\right(\xi ,\omega ){\mathrm{e}}^{\delta \xi }+D_2(\xi ,\omega \left){\mathrm{e}}^{-\delta \xi }\right)\\ & =\delta Q(\xi ,\omega {)}^2\left(D_1\right(\xi ,\omega ){\mathrm{e}}^{\delta \xi }-D_2(\xi ,\omega \left){\mathrm{e}}^{-\delta \xi }\right).\end{array}$This linear equation in $D_1(\xi ,\omega )$ and $D_2(\xi ,\omega )$ is solved easily and by (30(30) $\begin{array}{r}\tilde{f}(\xi ,\omega )=\frac{-\xi }{C\xi +\omega }\tilde{\psi }(\xi ,1,\omega )=\frac{-\xi }{Q(\xi ,\omega )}\tilde{\psi }(\xi ,1,\omega ).\end{array}$(30) ) we have $\begin{array}{rl}\text{∼}f(\xi ,\omega )& =\frac{-\xi }{Q(\xi ,\omega )}\text{∼}\psi (\xi ,1,\omega ))\\ & =\frac{\delta Q(\xi ,\omega )(\omega +\frac{B}{\sqrt{gd}}\xi )}{\left(\delta Q\right(\xi ,\omega {)}^2\cosh \left(\delta \xi \right)-(\xi +\sqrt{\frac{d}{g}}AQ(\xi ,\omega \left)\right)\sinh \left(\delta \xi \right))}\text{∼}h(\xi ,\omega ),\end{array}$which we will write as (31) $\begin{array}{r}\text{∼}f(\xi ,\omega )=\frac{{\delta }^{2}Q(\xi ,\omega )(\omega +\frac{B}{\sqrt{gd}}\xi )\text{∼}h(\xi ,\omega )}{\left({\delta }^{2}Q\right(\xi ,\omega {)}^2\cosh \left(\delta \xi \right)-(\delta \xi +\sqrt{\frac{d}{g}}A\delta Q(\xi ,\omega \left)\right)\sinh \left(\delta \xi \right))}.\end{array}$(31) We will also write this as (32) $\begin{array}{r}(S^2+\frac{D-i\sqrt{\frac{d}{g}}AS}{\delta }\tanh (\delta D\left)\right)f=\frac{S({\mathrm{\partial }}_t+i\frac{B}{\sqrt{gd}}D)}{\cosh \left(\delta D\right)}h.\end{array}$(32) One would need to consider the roots of the denominator in (31(31) $\begin{array}{r}\text{∼}f(\xi ,\omega )=\frac{{\delta }^{2}Q(\xi ,\omega )(\omega +\frac{B}{\sqrt{gd}}\xi )\text{∼}h(\xi ,\omega )}{\left({\delta }^{2}Q\right(\xi ,\omega {)}^2\cosh \left(\delta \xi \right)-(\delta \xi +\sqrt{\frac{d}{g}}A\delta Q(\xi ,\omega \left)\right)\sinh \left(\delta \xi \right))}.\end{array}$(31) ) in order to justify (32(32) $\begin{array}{r}(S^2+\frac{D-i\sqrt{\frac{d}{g}}AS}{\delta }\tanh (\delta D\left)\right)f=\frac{S({\mathrm{\partial }}_t+i\frac{B}{\sqrt{gd}}D)}{\cosh \left(\delta D\right)}h.\end{array}$(32) ) rigorously. However for the case A = 0, it turns out that the singularity cancels.

In the case A = 0 formula (32(32) $\begin{array}{r}(S^2+\frac{D-i\sqrt{\frac{d}{g}}AS}{\delta }\tanh (\delta D\left)\right)f=\frac{S({\mathrm{\partial }}_t+i\frac{B}{\sqrt{gd}}D)}{\cosh \left(\delta D\right)}h.\end{array}$(32) ) has the simpler form (33) $\begin{array}{r}(S^2+\frac{D}{\delta }\tanh (\delta D\left)\right)f=\frac{S^2}{\cosh \left(\delta D\right)}h.\end{array}$(33) We will extract a formula for f whenever it solves the following type of equation. Let $\tau :\mathbb{R}\to \mathbb{R}$ be an even function such that $\tau \left(\xi \right)>0$ for all $\xi \ne 0$. Note that in particular (25(25) $\begin{array}{r}m(-\xi )=\overline{m\left(\xi \right)}\end{array}$(25) ) is satisfied. We try to derive a formula for the solution f of (34) $\begin{array}{r}\begin{array}{rl}(S^2+\tau (D\left)\right)f& =\theta \\ \underset{t\to -\mathrm{\infty }}{lim}f& =0\end{array}\},\end{array}$(34) where θ is a given forcing term. We will assume f and θ to be in $\mathcal{S}({\mathbb{R}}^2,\mathbb{R})$ and the limit in (34(34) $\begin{array}{r}\begin{array}{rl}(S^2+\tau (D\left)\right)f& =\theta \\ \underset{t\to -\mathrm{\infty }}{lim}f& =0\end{array}\},\end{array}$(34) ) is adequate for functions in the Schwartz class [18]. We can decompose the operator $S^2+\tau \left(D\right)$ into (35) $\begin{array}{r}{S}^2+\tau \left(D\right)=(S+i\sqrt{\tau \left(D\right)})(S-i\sqrt{\tau \left(D\right)}),\end{array}$(35) since S and $i\sqrt{\tau \left(D\right)}$ commute. If we now set $\begin{array}{r}{f}_1=(S+i\sqrt{\tau \left(D\right)})f\end{array}$then we see by (35(35) $\begin{array}{r}{S}^2+\tau \left(D\right)=(S+i\sqrt{\tau \left(D\right)})(S-i\sqrt{\tau \left(D\right)}),\end{array}$(35) ) and the fact that we can commute the two factors that $\begin{array}{r}(S-i\sqrt{\tau \left(D\right)})f_1=\theta .\end{array}$We can now filter $f_1$ by the group ${\mathrm{e}}^{it\sqrt{\tau \left(D\right)}}$: Let $f_2={\mathrm{e}}^{-it\sqrt{\tau \left(D\right)}}{f}_1$. A calculation shows that $\begin{array}{r}S{f}_2=(C{\mathrm{\partial }}_x+{\mathrm{\partial }}_t)f_2={\mathrm{e}}^{-it\sqrt{\tau \left(D\right)}}\theta .\end{array}$This is an inhomogeneous transport equation with solution $\begin{array}{r}{f}_2(x,t)={\int }_{-\mathrm{\infty }}^t{\mathrm{e}}^{-is\sqrt{\tau \left(D\right)}}\theta (x+C(s-t),s)\mathrm{d}s.\end{array}$This yields $\begin{array}{r}\left(\sqrt{\tau \left(D\right)}f\right)(x,t)=\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(D\right)}}\theta \right(x+C(s-t),s\left)\mathrm{d}s\right).\end{array}$If we write out the Fourier multipliers explicitly, we get (36) $\begin{array}{r}f(x,t)=\frac{1}{\sqrt{2\pi }}\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\int }_{\mathbb{R}}\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}}{\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }\hat{\theta }\right(\xi ,s\left)\mathrm{d}\xi \mathrm{d}s\right).\end{array}$(36) We will deal with the possible singularity in the for us relevant case: We have $\tau \left(\xi \right)=\frac{\xi }{\delta }\tanh \left(\delta \xi \right)$ and $\theta =\frac{S^2}{\cosh \left(\delta D\right)}h$. We see that τ is even and $\tau \left(\xi \right)>0$ for $\xi \ne 0$. We have furthermore that τ is smooth with derivatives of at most polynomial growth at infinity. We will use that (37) $\begin{array}{r}\underset{\xi \to 0}{lim}\frac{|\xi |}{\sqrt{\tau \left(\xi \right)}}=1.\end{array}$(37) Now observe that the function $\frac{1}{\sqrt{\tau \left(D\right)}}\theta$ is real-valued and hence we can subtract $\left(\frac{1}{\sqrt{\tau \left(D\right)}}\theta \right)(0,\cdot )$ in (36(36) $\begin{array}{r}f(x,t)=\frac{1}{\sqrt{2\pi }}\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\int }_{\mathbb{R}}\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}}{\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }\hat{\theta }\right(\xi ,s\left)\mathrm{d}\xi \mathrm{d}s\right).\end{array}$(36) ) to obtain (38) $\begin{array}{r}f(x,t)=\frac{1}{\sqrt{2\pi }}\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\int }_{\mathbb{R}}\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }-1}{\sqrt{\tau \left(\xi \right)}}\hat{\theta }\right(\xi ,s\left)\mathrm{d}\xi \mathrm{d}s\right).\end{array}$(38) Clearly, $\hat{\theta }(\xi ,s)$ remains bounded around 0. For the quotient, note that we can rewrite it as follows: $\begin{array}{r}\left|\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }-1}{\sqrt{\tau \left(\xi \right)}}\right|=\left|\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }-1}{\xi }\right|\left|\frac{\xi }{\sqrt{\tau \left(\xi \right)}}\right|\end{array}$The second factor converges to 1 by (37(37) $\begin{array}{r}\underset{\xi \to 0}{lim}\frac{|\xi |}{\sqrt{\tau \left(\xi \right)}}=1.\end{array}$(37) ). For the first factor, note that the function $\begin{array}{r}\xi \mapsto {\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }\end{array}$is left and right differentiable in 0 and hence the first factor remains bounded as $\xi \to 0$ for any s. We conclude that (38(38) $\begin{array}{r}f(x,t)=\frac{1}{\sqrt{2\pi }}\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\int }_{\mathbb{R}}\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }-1}{\sqrt{\tau \left(\xi \right)}}\hat{\theta }\right(\xi ,s\left)\mathrm{d}\xi \mathrm{d}s\right).\end{array}$(38) ) holds (at least pointwise).

Note that we cannot do the same for the case $A\ne 0$: The operator on the left hand side of (32(32) $\begin{array}{r}(S^2+\frac{D-i\sqrt{\frac{d}{g}}AS}{\delta }\tanh (\delta D\left)\right)f=\frac{S({\mathrm{\partial }}_t+i\frac{B}{\sqrt{gd}}D)}{\cosh \left(\delta D\right)}h.\end{array}$(32) ) is given by $\begin{array}{r}{S}^2+(1-C^2)\frac{D}{\delta }\tanh \left(\delta D\right)-\frac{iC}{\delta }{\mathrm{\partial }}_t\tanh \left(\delta D\right)\end{array}$and there does not seem to be an apparent way to obtain a root for $\begin{array}{r}(1-C^2)\frac{D}{\delta }\tanh \left(\delta D\right)-\frac{iC}{\delta }{\mathrm{\partial }}_t\tanh \left(\delta D\right).\end{array}$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 $\delta \to 0$. 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 $\delta \to \mathrm{\infty }$) 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 $\delta \to 0$, (33(33) $\begin{array}{r}(S^2+\frac{D}{\delta }\tanh (\delta D\left)\right)f=\frac{S^2}{\cosh \left(\delta D\right)}h.\end{array}$(33) ) becomes (39) $\begin{array}{r}(S^2-{\mathrm{\partial }}_x^2)f=S^{2}h.\end{array}$(39) The initial conditions are on the one hand $f(x,0)=0$ from (15(15) $\begin{array}{r}\begin{array}{rl}& \{\begin{array}{l}{u}_x+v_y=0\\ u_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}})u_x+\sqrt{\frac{d}{g}}Av=-p_x\\ {\delta }^2(v_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}})v_x)=-p_y\\ u_y-{\delta }^2{v}_x=0\end{array}\mathrm{i}\mathrm{n}0<y<1\\ & \begin{array}{ll}p=f& \mathrm{o}\mathrm{n}y=1\\ v=f_t+(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_x& \mathrm{o}\mathrm{n}y=1\\ v=h_t+\frac{B}{\sqrt{gd}}{h}_x& \mathrm{o}\mathrm{n}y=0\end{array}\\ & f(x,0)=0,u(x,y,0)=0,v(x,y,0)=0\end{array}\}.\end{array}$(15) ). On the other hand, this gives $f_x(x,0)=0$ and inserting this in the sixth relation in (15(15) $\begin{array}{r}\begin{array}{rl}& \{\begin{array}{l}{u}_x+v_y=0\\ u_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}})u_x+\sqrt{\frac{d}{g}}Av=-p_x\\ {\delta }^2(v_t+(\sqrt{\frac{d}{g}}Ay+\frac{B}{\sqrt{gd}})v_x)=-p_y\\ u_y-{\delta }^2{v}_x=0\end{array}\mathrm{i}\mathrm{n}0<y<1\\ & \begin{array}{ll}p=f& \mathrm{o}\mathrm{n}y=1\\ v=f_t+(\sqrt{\frac{d}{g}}A+\frac{B}{\sqrt{gd}})f_x& \mathrm{o}\mathrm{n}y=1\\ v=h_t+\frac{B}{\sqrt{gd}}{h}_x& \mathrm{o}\mathrm{n}y=0\end{array}\\ & f(x,0)=0,u(x,y,0)=0,v(x,y,0)=0\end{array}\}.\end{array}$(15) ) for t = 0 yields $f_t(x,0)=0$. By using Duhamel's principle (similar to [19, p. 80–81]), we get the formula (40) $\begin{array}{r}f(x,t)=\frac{1}{2}{\int }_0^t{\int }_{x-(t-s)(C+1)}^{x-(t-s)(C-1)}{S}^{2}h(r,s)\mathrm{d}r\mathrm{d}s.\end{array}$(40) We will assume that h can be separated in the following way: $h(x,t)=a\left(t\right)b\left(x\right)$ for $a\in C^2(\mathbb{R},[0,\mathrm{\infty }\left)\right)$, $b\in C^2(\mathbb{R},\mathbb{R})$ such that $a\left(t\right)=0$ for $t\le 0$ and $a\left(t\right)=1$ for $t>t_0$ (where $t_0$ represents the duration of the earthquake), and with $b\left(x\right)=0$ for $x\notin (-L,L)$ modelling a localised tsunami source. Inserting this into (40(40) $\begin{array}{r}f(x,t)=\frac{1}{2}{\int }_0^t{\int }_{x-(t-s)(C+1)}^{x-(t-s)(C-1)}{S}^{2}h(r,s)\mathrm{d}r\mathrm{d}s.\end{array}$(40) ) yields the formula $\begin{array}{rl}f(x,t)& =a\left(t\right)b\left(x\right)\\ & +\frac{1}{2}{\int }_0^{t}a\left(s\right)\left[b^{\prime }\right(x-(t-s)(C-1))-b^{\prime }(x-(t-s)(C+1)\left)\right]\mathrm{d}s,\end{array}$which is a lot simpler. In the limiting case $t_0\searrow 0$, 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) $\begin{array}{r}f(x,t)=\frac{C^2}{C^2-1}b\left(x\right)+\frac{b(x-t(1+C\left)\right)}{2(1+C)}+\frac{b(x+t(1-C\left)\right)}{2(1-C)},\end{array}$(41) for $x\in \mathbb{R}$ and t>0. We can draw some insightful conclusions:

• at each instant t>0 after initiation, the generated wave is localised as $f(x,t)=0$ for $|x|\ge L+t(1+C)$;

• the surface wave consists of one stationary part (which can be disregarded, since $C\ll 1$) 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 $\sqrt{gd}-B$ in the physical variables) and the wave travelling to the right moves with non-dimensionalized speed 1 + C (corresponding to $\sqrt{gd}+B$ 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 $\frac{1}{2(1-C)}$ and $\frac{1}{2(1+C)}$, respectively.

This behaviour is illustrated in Figures 2 and  3. 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 [11]. Since $B\ll \sqrt{gd}$, i.e. $C\ll 1$, the model for $B\ne 0$ 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.

Figure 2. Model shortly after the instantaneous upward and downward thrust.

Figure 3. Expected long-term behaviour.

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 (38(38) $\begin{array}{r}f(x,t)=\frac{1}{\sqrt{2\pi }}\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\int }_{\mathbb{R}}\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }-1}{\sqrt{\tau \left(\xi \right)}}\hat{\theta }\right(\xi ,s\left)\mathrm{d}\xi \mathrm{d}s\right).\end{array}$(38) ). Ultimately, we will see that this approach does not work, since the method only works for too large t, see (48(48) $\begin{array}{r}T\gg \frac{\lambda }{\sqrt{gd}{\delta }^2}\approx 2.5\times {10}^6\mathrm{s}\approx 700\mathrm{h},\end{array}$(48) ).

We cannot do the same for B = 0 or the general case, since we lack a formula which would need to correspond to (36(36) $\begin{array}{r}f(x,t)=\frac{1}{\sqrt{2\pi }}\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\int }_{\mathbb{R}}\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}}{\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }\hat{\theta }\right(\xi ,s\left)\mathrm{d}\xi \mathrm{d}s\right).\end{array}$(36) ) or (38(38) $\begin{array}{r}f(x,t)=\frac{1}{\sqrt{2\pi }}\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\int }_{\mathbb{R}}\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }-1}{\sqrt{\tau \left(\xi \right)}}\hat{\theta }\right(\xi ,s\left)\mathrm{d}\xi \mathrm{d}s\right).\end{array}$(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 (38(38) $\begin{array}{r}f(x,t)=\frac{1}{\sqrt{2\pi }}\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\int }_{\mathbb{R}}\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }-1}{\sqrt{\tau \left(\xi \right)}}\hat{\theta }\right(\xi ,s\left)\mathrm{d}\xi \mathrm{d}s\right).\end{array}$(38) ). For this we try to use the stationary-phase principle. Recall that we set $\tau \left(\xi \right)=\frac{\xi }{\delta }\tanh \left(\delta \xi \right)$ and $\theta =\frac{S^2}{\cosh \left(\delta D\right)}h$. We first rewrite (38(38) $\begin{array}{r}f(x,t)=\frac{1}{\sqrt{2\pi }}\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\int }_{\mathbb{R}}\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }-1}{\sqrt{\tau \left(\xi \right)}}\hat{\theta }\right(\xi ,s\left)\mathrm{d}\xi \mathrm{d}s\right).\end{array}$(38) ) as $\begin{array}{rl}f(x,t)& =\frac{1}{\sqrt{2\pi }}\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\int }_{\mathbb{R}}\frac{{\mathrm{e}}^{i(t-s)\sqrt{\tau \left(\xi \right)}}{\mathrm{e}}^{ix\xi }{\mathrm{e}}^{iC(s-t)\xi }-1}{\sqrt{\tau \left(\xi \right)}}\hat{\theta }\right(\xi ,s\left)\mathrm{d}\xi \mathrm{d}s\right)\\ & =\frac{1}{\sqrt{2\pi }}\mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\int }_{\mathbb{R}}\frac{{\mathrm{e}}^{-is(\sqrt{\tau \left(\xi \right)}-C\xi )}{\mathrm{e}}^{it(\sqrt{\tau \left(\xi \right)}+(\mathcal{X}-C\left)\xi \right))}-1}{\sqrt{\tau \left(\xi \right)}}\hat{\theta }\right(\xi ,s\left)\mathrm{d}\xi \mathrm{d}s\right),\end{array}$where $\mathcal{X}=\frac{x}{t}$. The derivative of the phase factor $\sqrt{\tau \left(\xi \right)}+(\mathcal{X}-C)\xi$ is given by $\begin{array}{r}\left[\sqrt{\tau }{]}^{\prime }\right(\xi )+\mathcal{X}-C\end{array}$and one can check that (42) $\begin{array}{r}{\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}=\frac{\sinh \left(\delta \xi \right)\cosh \left(\delta \xi \right)+\delta \xi }{2\delta {\cosh }^2\left(\delta \xi \right)}\sqrt{\frac{\delta \cosh \left(\delta \xi \right)}{\xi \sinh \left(\delta \xi \right)}}\end{array}$(42) and (43) $\begin{array}{r}{\mathrm{\partial }}_{\xi }^2\sqrt{\tau \left(\xi \right)}=-\frac{1}{{\tau }^{\frac{3}{2}}\left(\xi \right)}\frac{{(\sinh (\delta \xi )\cosh (\delta \xi )-\delta \xi )}^2+4{\delta }^2{\xi }^2{\sinh }^2\left(\delta \xi \right)}{4{\delta }^2{\cosh }^4\left(\delta \xi \right)}<0.\end{array}$(43) One computes the limits $\begin{array}{rl}& \underset{\xi \searrow 0}{lim}{\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}=1,\\ & \underset{\xi \nearrow 0}{lim}{\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}=-1,\\ & \underset{\xi \to \mathrm{\infty }}{lim}{\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}=0,\\ & \underset{\xi \to -\mathrm{\infty }}{lim}{\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}=0,\end{array}$and concludes that $\xi \mapsto {\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}$ is strictly decreasing on $(0,\mathrm{\infty })$ from the asymptotic value 1 towards the asymptotic value 0 and on $(-\mathrm{\infty },0)$ from the asymptotic value 0 to the asymptotic value $-1$. So we conclude that $\begin{array}{r}\left[\sqrt{\tau }{]}^{\prime }\right({\xi }_0)+\mathcal{X}-C=0\end{array}$is possible if $\mathcal{X}-C\in (-1,0)\cup (0,1)$ for exactly one ${\xi }_0$ for given $\mathcal{X}$. So as long as $\mathcal{X}$ is in the corresponding range, we can write ${\xi }_0={\xi }_0\left(\mathcal{X}\right)$. The stationary-phase principle gives $\begin{array}{rl}f(x,t)& \sim \frac{1}{\sqrt{t\tau \left({\xi }_0\right)|\left[\sqrt{\tau }{]}^{″}\right({\xi }_0)|}}\cdot \\ & \mathrm{I}\mathrm{m}\left({\int }_{-\mathrm{\infty }}^t{\mathrm{e}}^{-is(\sqrt{\tau \left({\xi }_0\right)}-C{\xi }_0)}\hat{\theta }\right({\xi }_0,s\left){\mathrm{e}}^{it(\sqrt{\tau \left({\xi }_0\right)}+(\mathcal{X}-C){\xi }_0+\sigma \frac{\pi }{4}))}\mathrm{d}s\right)\end{array}$for large t, where σ is the sign of $\left[\sqrt{\tau }{]}^{″}\right({\xi }_0)$, which is just $-1$ by (43(43) $\begin{array}{r}{\mathrm{\partial }}_{\xi }^2\sqrt{\tau \left(\xi \right)}=-\frac{1}{{\tau }^{\frac{3}{2}}\left(\xi \right)}\frac{{(\sinh (\delta \xi )\cosh (\delta \xi )-\delta \xi )}^2+4{\delta }^2{\xi }^2{\sinh }^2\left(\delta \xi \right)}{4{\delta }^2{\cosh }^4\left(\delta \xi \right)}<0.\end{array}$(43) ). We see that this is a Fourier transform (with respect to time) and hence get $\begin{array}{rl}f(x,t)& \sim \frac{\sqrt{2\pi }}{\sqrt{t\tau \left({\xi }_0\right)|\left[\sqrt{\tau }{]}^{″}\right({\xi }_0)|}}\cdot \\ & \mathrm{I}\mathrm{m}\left(\text{∼}\theta \right({\xi }_0,\sqrt{\tau \left({\xi }_0\right)}-C{\xi }_0\left){\mathrm{e}}^{it(\sqrt{\tau \left({\xi }_0\right)}+(\mathcal{X}-C){\xi }_0-\frac{\pi }{4}))}\right).\end{array}$We have $\begin{array}{rl}\text{∼}\theta ({\xi }_0,\sqrt{\tau \left({\xi }_0\right)}-C{\xi }_0)& =-\frac{Q({\xi }_0,\sqrt{\tau \left({\xi }_0\right)}-C{\xi }_0{)}^2}{\cosh \left(\delta {\xi }_0\right)}\text{∼}h({\xi }_0,\sqrt{\tau \left({\xi }_0\right)}-C{\xi }_0)\\ & =-\frac{\tau \left({\xi }_0\right)}{\cosh \left(\delta {\xi }_0\right)}\text{∼}h({\xi }_0,\sqrt{\tau \left({\xi }_0\right)}-C{\xi }_0)\end{array}$and hence get (44) $\begin{array}{rl}f(x,t)& \sim \frac{-\sqrt{2\pi \tau \left({\xi }_0\right)}}{\cosh \left(\delta {\xi }_0\right)\sqrt{t|\left[\sqrt{\tau }{]}^{″}\right({\xi }_0)|}}\cdot \\ & \mathrm{I}\mathrm{m}\left({\mathrm{e}}^{it(\sqrt{\tau \left({\xi }_0\right)}+(\mathcal{X}-C){\xi }_0-\frac{\pi }{4}))}\text{∼}h\right({\xi }_0,\sqrt{\tau \left({\xi }_0\right)}-C{\xi }_0\left)\right).\end{array}$(44) Finally, we try to extract the asymptotics with respect to δ: We fix a $\xi >0$. A look at (42(42) $\begin{array}{r}{\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}=\frac{\sinh \left(\delta \xi \right)\cosh \left(\delta \xi \right)+\delta \xi }{2\delta {\cosh }^2\left(\delta \xi \right)}\sqrt{\frac{\delta \cosh \left(\delta \xi \right)}{\xi \sinh \left(\delta \xi \right)}}\end{array}$(42) ) gives that $\delta \mapsto {\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}$ is $C^{\mathrm{\infty }}$ and clearly even in δ. Since we have $\underset{\delta \searrow 0}{lim}{\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}=1$, the Taylor expansion with respect to δ around 0 has the form $\begin{array}{r}{\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}=1+\alpha \left(\xi \right){\delta }^2+O\left({\delta }^3\right),\end{array}$for some coefficient $\alpha \left(\xi \right)$. One could find this coefficient through tedious calculations, but one can instead notice that $2\alpha \left(\xi \right)$ is the coefficient of ${\delta }^2$ of the expansion of $({\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}{)}^2$. We have $\begin{array}{rl}{\left({\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}\right)}^2& =\frac{\sinh \left(\delta \xi \right)}{4\delta \xi \cosh \left(\delta \xi \right)}+\frac{1}{2{\cosh }^2\left(\delta \xi \right)}+\frac{\delta \xi }{4\sinh \left(\delta \xi \right){\cosh }^3\left(\delta \xi \right)}\\ & =1-{\delta }^2{\xi }^2+O\left({\delta }^3\right)\end{array}$and hence we conclude ${\mathrm{\partial }}_{\xi }\sqrt{\tau \left(\xi \right)}=1-\frac{1}{2}{\delta }^2{\xi }^2+O\left({\delta }^3\right)$. We conclude that for $\mathcal{X}-C\in (-1,0)$ the stationary-phase point ${\xi }_0>0$ with $\left[\sqrt{\tau }{]}^{\prime }\right({\xi }_0)+\mathcal{X}-C=0$ satisfies $\delta {\xi }_0=O\left(1\right)$. Inserting this into (43(43) $\begin{array}{r}{\mathrm{\partial }}_{\xi }^2\sqrt{\tau \left(\xi \right)}=-\frac{1}{{\tau }^{\frac{3}{2}}\left(\xi \right)}\frac{{(\sinh (\delta \xi )\cosh (\delta \xi )-\delta \xi )}^2+4{\delta }^2{\xi }^2{\sinh }^2\left(\delta \xi \right)}{4{\delta }^2{\cosh }^4\left(\delta \xi \right)}<0.\end{array}$(43) ) yields (45) $\begin{array}{r}|{\mathrm{\partial }}_{\xi }^2\sqrt{\tau \left({\xi }_0\right)}|=O\left(\delta \right).\end{array}$(45)

4.2.2. Some typical physical parameters

We want to see whether stationary-phase analysis, using the formula $\begin{array}{rl}f(x,t)& \sim \frac{-\sqrt{2\pi \tau \left({\xi }_0\right)}}{\cosh \left(\delta {\xi }_0\right)\sqrt{t|\left[\sqrt{\tau }{]}^{″}\right({\xi }_0)|}}\cdot \\ & \mathrm{I}\mathrm{m}\left({\mathrm{e}}^{it(\sqrt{\tau \left({\xi }_0\right)}+(\mathcal{X}-C){\xi }_0-\frac{\pi }{4}))}\text{∼}h\right({\xi }_0,\sqrt{\tau \left({\xi }_0\right)}-C{\xi }_0\left)\right),\end{array}$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 [8]): $\begin{array}{r}a=1\mathrm{m},d=4\mathrm{k}\mathrm{m},\lambda =200\mathrm{k}\mathrm{m}.\end{array}$This leads to $ϵ=0.00025$ and $\delta =0.02$ which point to linear wave approximation. Applying the stationary-phase principle is justified if (46) $\begin{array}{r}t|{\sqrt{\tau }}^{″}\left(\xi \right)|\gg 1.\end{array}$(46) We have by (45(45) $\begin{array}{r}|{\mathrm{\partial }}_{\xi }^2\sqrt{\tau \left({\xi }_0\right)}|=O\left(\delta \right).\end{array}$(45) ), $|{\sqrt{\tau }}^{″}\left(\xi \right)|=O\left(\delta \right)$, which translates to (47) $\begin{array}{r}t\gg \frac{1}{{\delta }^2}.\end{array}$(47) In physical variables this means (48) $\begin{array}{r}T\gg \frac{\lambda }{\sqrt{gd}{\delta }^2}\approx 2.5\times {10}^6\mathrm{s}\approx 700\mathrm{h},\end{array}$(48) which takes way too long to be applicable. If one were to relax (47(47) $\begin{array}{r}t\gg \frac{1}{{\delta }^2}.\end{array}$(47) ) to a condition of the form $\begin{array}{r}t\ge \frac{1}{\delta },\end{array}$one would arrive at $\begin{array}{r}T\ge \frac{\lambda }{\sqrt{gd}\delta }\approx 5\times {10}^4\mathrm{s}\approx 14\mathrm{h}\end{array}$which would barely be applicable in the case where the epicentre of the seaquake is far off the shore. Nevertheless, relaxing (47(47) $\begin{array}{r}t\gg \frac{1}{{\delta }^2}.\end{array}$(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).