Influence of timescales on the generation of seismic tsunamis

European Journal of Mechanics B/Fluids 65 (2017) 257–273

Contents lists available at ScienceDirect

European Journal of Mechanics B/Fluids journal homepage: www.elsevier.com/locate/ejmflu

Influence of timescales on the generation of seismic tsunamis Marine Le Gal a,∗ , Damien Violeau a,b , Michel Benoit c a

Saint-Venant Hydraulics Laboratory, University Paris-Est, Joint research unit EDF - Cerema - Ecole des Ponts, Chatou, France

b

EDF R&D, Chatou, France Institut de Recherche sur les Phénomènes Hors-Equilibre (IRPHE, UMR 7342), Aix-Marseille Univ., CNRS, Centrale Marseille, Marseilles, France

c

article

info

Article history: Received 6 April 2016 Available online 27 April 2017 Keywords: Tsunami seismic generation Linear model Resonance

abstract A new semi-analytical solution is derived in a linear potential framework for simulating tsunami-like water waves caused by a dynamic sea ground deformation. The movement generating the tsunami is represented by a simplified model of the sea-floor displacement caused by a schematic seismic event. This displacement generalises the earlier work of Hammack (1973) and Todorovska and Trifunac (2001), in the sense that both the rise time (tr ) in the vertical direction and the rupture propagation velocity (vp ) in the horizontal direction are simultaneously included in the model. The influence of both parameters is investigated by applying the linear solution to a large range of parameter values. An analysis of the maximum amplitude of generated waves allows the identification of a resonance √ for no rise time (tr ≈ 0) and a rupture propagation velocity close to the long wave celerity (vp ≈ gh), corroborating the resonance mechanism identified by Todorovska and Trifunac (2001). An empirical relationship is proposed to estimate the amplification factor as a function of the fault length. In a second step, an energetic ratio is defined representing the quantity of potential energy that is not taken into account properly when using the shallow water equations. The energetic ratio is estimated as the fraction of potential energy associated with wavelengths shorter than the long wave limit. This ratio is largest for rapid vertical motions, in particular in the vicinity of the resonant condition, suggesting that models based on the shallow water equations should not be used for these conditions. Numerical simulations using a nonlinear and dispersive wave propagation model are also performed, showing that the resonance effect still exists when nonlinear effects are included. Finally, the tsunami event of 1947 near New Zealand is analysed and shown to be close to the resonant condition, with seismic activity generating a free surface wave that could be much larger than the amplitude of the ground motion. © 2017 Elsevier Masson SAS. All rights reserved.

1. Introduction In the recent past, the catastrophic impacts of tsunamis have reminded us how important it is to study their dynamics. These natural events can cause human losses and large material damages, and have additional societal consequences as experienced, for instance, during the December 2004 Indonesian tsunami, and more recently during the March 2011 Japanese tsunami. Among all of the tsunamigenic generation mechanisms, seismic generation is far more frequent than other causes (landslides, volcanic eruptions, meteorological causes, meteorites). In this case, the wave is generated by a dynamic deformation of the

∗

Corresponding author. E-mail addresses: [email protected] (M. Le Gal), [email protected] (D. Violeau), [email protected] (M. Benoit). http://dx.doi.org/10.1016/j.euromechflu.2017.03.008 0997-7546/© 2017 Elsevier Masson SAS. All rights reserved.

sea floor. This generation mechanism has been the subject of numerous studies, and several models have been developed. Braddock et al. [1] proposed a linear model by describing the ground deformation with series of orthogonal functions. Later, Sabatier [2] built 2D and 3D linear models for the generation of a wave on a slope. Lee et al. [3] worked on a combined theoretical, numerical and experimental study. They used generalised Boussinesq equations and Korteweg–de Vries equations to generate nonlinear waves. They also identified a critical Froude number that Todorovska and Trifunac [4] highlighted with a linear propagation model. Levin and Nosov [5] completed a comprehensive review of existing work with studies of different effects. Usually, to model the generation of a tsunami, a data inversion is first performed to obtain the characteristics of the fault (length, width, orientation and slip), see Ji et al. [6]. The second step is to calculate the vertical deformation of the sea floor using the method of Okada [7]. Then, the perturbation of the free surface is assumed to be instantaneous and to have the same

258

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

amplitude as that of the sea floor. Thus the initial condition for the tsunami wave is given by this free surface deformation, which then propagates in the considered spatial domain [8]. This traditional approach neglects a certain number of effects such as the horizontal displacement, the compressibility of the water, nonlinear effects, the rotation of the Earth and the dynamics of the bottom deformation. The present work only considers the impact of the timescale parameters, based on the work of Hammack [9] and Todorovska and Trifunac [4]. Hammack [9] carried out a theoretical, numerical and experimental study of the impact of the rising time tr of the ground deformation (the time for the vertical deformation). Since Hammack, few studies have been led about the impact of this rise time, among them, Jamin et al. [10] set up an experimental and theoretical study and pointed out a temporal high-pass filter and spatial low-pass filter that are not taken into account in tsunami warning system. The same year, the impact of the rise time and seamounts on wave trapping was analysed by Stefanakis et al. [11] using an analytical solution. While the rise time is associated to the temporal vertical motion, a temporal horizontal parameter can also be defined: the rupture propagation velocity vp . Todorovska and Trifunac [4] evaluated its impact. They identified an amplification of the free surface deformation when this parameter coincides with the long wave celerity. Dutykh [12], as well as Dutykh et al. [13] and Dutykh et al. [14], used theoretical and numerical approaches to study the dynamic generation of tsunamis. These authors developed a linear solution to study the differences from a passive or active generation, and used the finite fault method to work on a real case tsunami. They also investigated the effects of the horizontal displacements. The purpose of the present paper is to identify the possible impact of the timescales on the generated wave. We investigate the influence of vertical and horizontal motion simultaneously, unlike Hammack [9] and Todorovska and Trifunac [4]. A dimensional analysis is first introduced in Section 2, to identify the dominant non-dimensional parameters to study. From that, a linear solution is developed, in Section 3, starting from Euler’s equations in a simple case. In order to validate the calculation of this theoretical solution, comparisons are made with numerical simulations. The developed solution is then used to study the amplitude of waves generated by bottom motion in Section 4, and to assess and discuss the validity of using the shallow water assumption in Section 5. These features are analysed for a range of values of the rise time tr and rupture velocity vp . Analyses of both linear solution (Section 4) and nonlinear numerical simulations (Section 6) show that Todorovska and Trifunac’s resonance only occurs when the rise time is small and that it persists when nonlinear effects are included. Finally, in Section 7, the present analyses are applied to discuss a tsunami event that occurred in Offshore Poverty Bay close to New Zealand in 1947. The main conclusions of the present study are summarised in Section 8. 2. Dimensional analysis In order to identify the dominant parameters of schematic tsunami generation, a dimensional analysis is first applied. The considered fluid domain is presented in Fig. 1. The analysis is limited to one horizontal direction, with an unbounded domain in the wave propagation direction x. Coriolis effects associated with the rotation of the Earth are not considered. The wavelength of a tsunami is much larger than the capillary length scale, thus surface tension effects are not taken into account. Moreover, the flow is assumed to be irrotational, and the fluid is assumed to be incompressible and inviscid. It is bounded above by the free surface at z = η(x, t ), and below by the sea bottom at z = −h + ζ (x, t ), where h is the initial water depth and ζ the deformation of the

Fig. 1. Definition of the fluid domain and coordinate system (x, z ), with h the initial depth, ζ (x, t ) the deformation of the sea floor and η(x, t ) the deformation of the free surface. Table 1 Variables and units of the problem.

η

Variables Units

h [m]

[m]

ζ0

g [m s−2 ]

[m]

tr [s]

L [m]

vp

[m s−1 ]

sea bottom. Initially, the fluid is at rest, with flat free surface and bottom elevations: η(x, t = 0) = 0 and ζ (x, t = 0) = 0. To simplify the problem, the ground motion is reduced to an uplift of a rectangular block that propagates with a given velocity vp along the x axis for x ∈ [0; L]. The vertical deformation reaches the maximum amplitude ζ0 in a finite time, the rising time tr , along the z axis. Fig. 2 gives a schematic representation of this movement. In Section 3.2, we will analyse the case of a sinusoidal rise motion. The movement of the sea bottom is defined by the parameters: ζ0 , tr , L and vp . The variables of the problem in the gravity field g are these four parameters with the free surface elevation η and the water depth h. The Vaschy–Buckingham, or Π theorem (see [15,16]) is used to identify the non-dimensional parameters. There are two independent units (time [s], and space [m]) and seven variables summarised in Table 1. Thus five non-dimensional parameters can be defined:

η∗ =

η

, 

ζ0∗ =

h

tr∗ = tr

g h

,

ζ0

(1)

h

L∗ =

L

(2)

h

vp vp∗ = √ .

(3)

gh

The free surface elevation depends on x and  t, thus nong x ∗ ∗ dimensional space and time: x = h and t = t h should be included in the parameters of the system (see sub-Section 3.2). To be consistent with previous studies, the non-dimensional t∗

√

temporal parameters are τ ∗ = Lr∗ = tLr gh, as chosen by Hammack [9] and vp∗ , as in [4]. τ ∗ represents the ratio between the vertical timescale tr and the time that the wave takes to propagate over the distance of deformation L, and vp∗ is the ratio between the horizontal timescale vp and the long wave celerity c = will thus seek to express the free surface elevation as:

  η∗ = Φ ζ ∗ , L∗ , τ ∗ , vp∗ , t ∗ , x∗ .

√

gh. We (4)

Hammack [9] worked only on the effects of τ , and Todorovska and Trifunac [4] only on the effects of vp∗ . In this study, the effects of both τ ∗ and vp∗ are simultaneously investigated. Note that Dutykh [12] also studied the impact of the Ursell number, defined by S = η∗ L∗2 with our notations. Usually, the horizontal dimension of the seismic source exceeds the water depth: L∗ ≫ 1. Based on the analysis of past events, L∗ can vary between 16 as for the 1946 Aleutian event [17] and 300 as for the 2004 Sumatra–Andaman event [18]. If the wavelength of the generated wave is considered equivalent to L as suggested by ∗

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

Levin and Nosov [5] and Wu [19], L∗ varies between 10 and 5000 in deep ocean [19]. Then, a long wave theory could be used to study the tsunami, as done by Wu [19] among others, but it will not be developed in this work, and moreover, its limit will be discussed in Section 5 with the analysis of consequences implied by the Shallow Water approximation.

259

This transform is applied to the system (6a), (6d), (7), giving:

φ¯ zz − k2 φ¯ = 0, φ¯ z = −

s

2

g

(14a)

φ¯

z = 0,

(14b)

3. Linear theory

φ¯ z = sζ¯ sφ¯ = −g η¯

3.1. Methodology

This system can be resolved analytically in (k, s) space to obtain:

We review here the method presented by Hammack [9] and also described by Todorovska and Trifunac [4] and Dutykh [12], among others. We consider the motion of the fluid domain defined in the previous section. Since the fluid is assumed to be incompressible and the flow irrotational, there exists a velocity potential φ(x, z , t ) that satisfies the Laplace equation and the boundary conditions as follows:

▽2 φ = 0, φz = ηt + φx ηx

(5a) z = η(x, t ),

(5b)

φt + | ▽ φ|2 + g η = 0

z = η(x, t ),

(5c)

φz = ζt + φx ζx

z = −h + ζ (x, t ),

(5d)

1

2

where subscripts t, x indicate partial derivatives. Eqs. (5b) and (5c) are the kinematic and dynamic free surface boundary conditions, respectively, and Eq. (5d) is the kinematic bottom boundary condition. By assuming small perturbations of the free surface and the bottom, the boundary conditions are linearised and the previous system of equations becomes:

▽2 φ = 0, φz = ηt φt + g η = 0 φz = ζt

(6a) z = 0,

(6b)

z = 0,

(6c)

z = −h.

(6d)

Note that the three linearised boundary conditions now apply at the undisturbed horizontal surfaces, i.e. at z = 0 for (6b) and (6c) and z = −h for (6d). Combining (6b) and (6c), the following free surface boundary condition is obtained:

φtt + g φz = 0 at z = 0.

F (f ) = f˜ (k) =

f (x)e−ikx dx,



1 2π

(8)

f˜ (k)eikx dk,

(9)

ℜ

and the Laplace transform of f (t ) for t > 0 is:

L(f ) = f(s) =



+∞

f (t )e−st dt ,

L−1 (f) = f (t ) =

γ +i∞



i2π

γ −i∞

F L(f ) = f¯ (k, s) =

e

f(s)est ds,

−ikx

(11)

=

2π

e ℜ

ikx

f ( x, t ) e

−st

dx

dt ,

1 i2π



γ +i∞ γ −i∞

f¯ (k, s)est dkds.

(16)

The final formula for the free surface deformation as a function of the bottom deformation is obtained in the physical space by taking the inverse Fourier–Laplace transform of (16):

η(x, t ) =

eikx



1 2π

ℜ

1



cosh kh i2π

γ +i∞ 2 ¯

s ζ (k, s)est

γ −i ∞

(s2 + ω2 )

dsdk.

(17)

3.2. Solution for a schematic uplift To go further we need to choose a schematic bed deformation

ζ (x, t ). Of special interest in the present study is the impact of tr and vp , as described in Section 2. The area of deformation is a segment of length L, horizontally deforming at velocity vp and taking a time tr of uplift, leading to:

ζ (x, t ) = ζ0 H (L − x)H (x)T (x, t ) (18)     x 1 x H + tr − t T (x, t ) = H t − vp vp 2       x x 1 − cos ωr t − +H t − − tr , vp vp (19) where ωr = t . Unlike Hammack [9] or Todorovska and Trir funac [4], both parameters are included in the present definition of T (x, t ). This definition coincides with the one of Hammack [9] when vp → ∞ (i.e. the horizontal movement of the bottom deformation is assumed to occur instantaneously). On the other hand, if tr vanishes (i.e. instantaneous vertical displacement), this definition corresponds to the 1D solution of Todorovska and Trifunac [4]. The Laplace–Fourier transform (12) is applied to ζ , and we find:

ζ0 2

(1 + e

−str

(12)

(13)

−L(ik+ vsp )

ω2 1−e ) 2 r 2 s(s + ωr ) ik +

s

.

(20)

vp

Substituting ζ¯ into Eq. (16), we obtain:

η( ¯ k, s) =

ζ0

s

2 (s2 + ω2 ) cosh kh −L(ik+ vsp )

0

F L−1 (f¯ ) = f (x, t )



, (s2 + ω2 ) cosh kh √ where ω = gk tanh kh is the linear dispersion relation.

+∞



ℜ

1

s2 ζ¯ (k, s)

η( ¯ k, s) =

ζ¯ (k, s) =

where γ is a real constant that insures the existence of the integral. Following Hammack [9], the Laplace–Fourier transform is built for a function f (x, t ) combining these two transforms:



  s2 −sg ζ¯ (k, s) cosh kz − sinh kz , (15) [s2 + gk tanh kh] cosh kh gk

thus,

(10)

0

1

(14d)

π

ℜ

F −1 (f˜ ) = f (x) =

(14c)

z = 0.

(7)

The 1D space Fourier transform of a function f (x) is defined as:



¯ k, z , s) = φ(

z = −h,

×

1−e

ik + vs p

.

(1 + e−str )

ωr2 s2 + ωr2 (21)

From here, only the Fourier transform of η, η˜ , can be found analytically. The integration over k can be performed analytically only in particular cases or by invoking some approximation due to the

260

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

Fig. 2. Definition of the sea floor uplift ζ (x, t ): 2(a) spatial profile of the sea floor and 2(b), the black dashed lines represent the motion defined by Todorovska and Trifunac [4], and the solid red lines are the motion used in this study. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

complexity of the expression with the associated dispersion relation, see Mei et al. [20]. To determine η˜ , we used the inverse Laplace transform and its properties (details of the calculations can be found in Appendix A). The non-dimensional variables are noted with ∗ . We define: k∗ = kh in addition to the definitions (3), leading to:

η˜ ∗ (k∗ , t ∗ ) =

ζ0∗

vp∗

2

cosh k∗

ωr∗2 − ω ∗2

ω ∗2 r

f ∗ (k∗ , t ∗ )



 +H (t − tr )f (k , t − tr )         ∗ ∗ L L   ∗ −ik∗ L∗ ∗ ∗ ∗  , (22)  − H t − e f k , t − × ∗ ∗  v v p p         L∗ L∗  −ik∗ L∗ ∗ ∗ ∗ ∗ ∗ ∗ f −H t − tr − ∗ e k , t − tr − ∗ vp vp ∗

∗

∗

∗

∗

∗

1



ik∗ vp∗ cos ω∗ t ∗ + ω∗ sin ω∗ t ∗ ω∗2 − k∗2 vp∗2  1 ∗ ∗ ∗ − ik∗ vp∗ e−ik vp t − ∗2 ωr − k∗2 vp∗2   ∗ ∗ ∗ × ik∗ vp∗ cos ωr∗ t ∗ + ωr∗ sin ωr∗ t ∗ − ik∗ vp∗ e−ik vp t .

4. Free surface deformation analysis (23)

And the non-dimensional free surface elevation is:

η (x , t ) = ∗

∗

∗

1



2π

eik

∗ x∗

η˜ ∗ dk∗ .

(24)

R

Again, one can verify that if tr tends to 0, this solution tends to the 1D solution of Todorovska and Trifunac [4], and if vp tends to infinity, the solution of Hammack [9] is recovered. As previously stated, this integral is difficult to perform, thus a numerical inverse Fourier transform is used to obtain the resulting free surface η(x, t ). In particular in this study, the integral is calculated numerically using Simpson’s method. 3.3. Treatment of singularities We note that the previous solution (22)–(23) of the Fourier transform of the free surface has some singularities that should be addressed. For this purpose, we take the limit of η˜ ∗ at the critical values of k∗ , defined by:



1. k∗1 = ±

ωr∗ , vp∗

2. k∗2a = 0, 3. if |vp∗ | < 1, k∗2b = p−1 ( ∗

−1

4. k3 = q

(ωr ), ∗

1

vp∗2

),

of the singularities k∗1 , k∗2 , k∗3 , respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

x and q(x) = gx tanh x. The limits are where p(x) = tanh x computed in Appendix B. Fig. 3 shows an example of these critical wavenumbers in the spectrum. One can see that |η˜ ∗ | remains continuous everywhere.

where: f ∗ (k∗ , t ∗ ) =

Fig. 3. An example of the free surface Fourier transform η˜ ∗ with critical k∗ for ∗ vp∗ = 10 and τ ∗ = 0.5 at t ∗ = T ∗ = vL ∗ + tr∗ . The black line is the modulus of the p ∗ Fourier transform |η˜ |. The vertical blue, green and red lines represent the locations

The linear solution derived in the previous section is now used to study the wave trains generated by a dynamical bottom motion of the form similar to expression (19). First (Section 4.1), the impact of the temporal parameters is investigated by looking the free surface deformation at the end of the ground motion, then particularly by looking the maximum amplitude in x. In a second step (Section 4.2), the propagation of the wave is studied in order to see if the impact of vp∗ and τ ∗ remains during this stage. Section 4.3 addresses the comparisons of the previous developed solution and numerical results from the code Misthyc. At this point of the study, these comparisons permit to validate the solution obtained with the calculations from the theory. 4.1. At the end of the bottom deformation By convention, we define the end of the generation period as ∗ the time when the ground stops moving, or t ∗ = T ∗ = vL ∗ + tr∗ . p

In the following analysis, the geometric parameter representing the horizontal extent of the bottom deformation is chosen as L∗ = 50. Since the model is linear, the value of ζ0∗ is not important (provided it is small), and therefore all subsequent results will show the ratio η∗ /ζ0∗ . The horizontal velocity vp∗ is varied between 0.5 and 50, and the rising time τ ∗ is varied between 0 and 5. As an example, if the dimensions of the 1992 Nicaragua event are considered [21], L = 250 km in an ocean of depth h = 5 km, the dimensionless numbers correspond to rupture velocities vp ∈ [110; 11 000] m s−1 and rise times tr ∈ [0; 5644] s.

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

(a) vp∗ = 0.5; τ ∗ = 0.

(b) vp∗ = 0.5; τ ∗ = 1.

(c) vp∗ = 0.5; τ ∗ = 2.

(d) vp∗ = 1; τ ∗ = 0.

(e) vp∗ = 1; τ ∗ = 1.

(f) vp∗ = 1; τ ∗ = 2.

(g) vp∗ = 2; τ ∗ = 0.

(h) vp∗ = 2; τ ∗ = 1.

(i) vp∗ = 2; τ ∗ = 2.

(j) vp∗ = 10; τ ∗ = 0.

(k) vp∗ = 10; τ ∗ = 1.

(l) vp∗ = 10; τ ∗ = 2.

(m) vp∗ = 50; τ ∗ = 0.

(n) vp∗ = 50; τ ∗ = 1.

(o) vp∗ = 50; τ ∗ = 2.

261

Fig. 4. Free surface profiles at t ∗ = T ∗ for L∗ = 50 and different values of vp∗ and τ ∗ (vp∗ increases from top to bottom and τ ∗ increases from left to right). The red line is the linear solution (22)–(24), and the dashed black line is the numerical simulation results from the linear version of the Misthyc code (see Section 4.3). The vertical scale of graph (d) for vp ∗ = 1 and τ ∗ = 0 differs from the others for clarity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 shows the free surface profiles for: vp∗ = 0.5, 1, 2, 10 and 50, and τ ∗ = 0, 1 and 2. The maximum free surface elevation ∗ ηmax /ζ0∗ varies with vp∗ and τ ∗ . If vp∗ = 50 (i.e. very large) and τ ∗ = 0, the free surface deformation is almost identical to the sea ∗ floor deformation (i.e. ηmax /ζ0 ≃ 1) since the deformation is nearly instantaneous in both horizontal and vertical directions. When τ ∗ increases, the wave begins to propagate before the end of the ground motion as noted by Jamin et al. [10], the wave amplitude is smaller, and the wave propagates in both directions (±x). When vp∗ decreases, an asymmetry appears: the wave propagating in the same direction as the ground deformation (+x) is larger. ∗ Profiles of ηmax /ζ0∗ as a function of τ ∗ , for several rupture ∗ velocities are plotted in Fig. 5(a). When τ ∗ increases, ηmax /ζ0∗ ∗ ∗ ∗ decreases. If τ > 2, the influence of vp is negligible. For vp =

10, 20, 50, the free surface profiles are almost identical, except when τ ∗ ∼ 1. Thus for vp∗ > 10, the motion in the x direction can be ∗ considered as nearly instantaneous. On Fig. 5(b), ηmax /ζ0∗ is plotted as a function of vp∗ for several values of rise time τ ∗ . In a similar way as Fig. 5(a), one can conclude that: (i) for vp∗ > 10, the curves reach a constant asymptote, (ii) around vp∗ = 1 and τ ∗ < 0.7, the resonance is clearly identifiable, (iii) for τ ∗ > 2, the curves are very flat, thus the influence of vp∗ seems negligible for these values of τ ∗ . ∗ Fig. 6 shows the dependency of ηmax /ζ0∗ on the temporal parameters in vp∗ − τ ∗ space. For small τ ∗ , the parameter vp∗ has a strong influence on the maximum amplitude as previously shown. As shown by Todorovska and Trifunac [4], a resonance appears for vp∗ in the range 0.7 to 2 (i.e. close to vp∗ = 1) generating ∗ a wave that reaches an amplitude ηmax /ζ0∗ larger than 1. The

262

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

∗ Fig. 5. Evolution of ηmax /ζ0∗ as a function of τ ∗ for different values of vp∗ (5(a), vertical axis in log scale) and as a function of vp∗ for different values of τ ∗ (5(b), horizontal axis in log scale) for L∗ = 50 at t ∗ = T ∗ .

surface at t ∗ = 100, 500 and 1000 are shown in Figs. 8, 9, 10, respectively, for vp∗ = 0.5, 1, 2, 10, 50 and τ ∗ = 0, 1, 2. Note that t ∗ may be interpreted as the ratio of the propagation distance to the initial water depth. At t ∗ = 100 (Fig. 8), two waves propagate in both directions for all parameters values. For vp∗ = 50, the two waves are almost symmetrical. When τ ∗ decreases, the amplitude decreases. For vp∗ = 1 and 2, the asymmetry of the propagation appears clearly. The wave propagating in the positive direction is larger. By superimposing the curves of vp∗ = 10 and vp∗ = 50 (not shown here), the results from vp∗ = 10 appear weakly non symmetrical, which confirms that the free surface deformation is asymmetrical when vp∗ decreases. These conclusions are also valid for t ∗ = 500 (Fig. 9) and t ∗ = 1000 (Fig. 10). Moreover, for τ ∗ = 0, the free surface deformations show the development of frequency dispersion increasing in time as shown by Hammack [9] and Stefanakis et al. [11]. Except for τ ∗ = 1 and vp∗ = 1 at t ∗ = 1000, this dispersion does not exist for larger τ ∗ . In the general case, each wave has the shape of a single hump of water, and the maximum amplitude does not change with time. For vp∗ = 1 and τ ∗ = 0 (resonance condition), the maximum amplitude decreases in time, but even at t ∗ = 1000, this maximum is still greater than 1: the impact of the resonance at the generation is still present. Fig. 11 shows the evolution of the free surface deformation η∗ /ζ0∗ as a function of x and t for vp∗ = 1 and τ ∗ = 0. In this graph, the frequency dispersion and wave asymmetry appear clearly. In the −x direction, a smaller but wider wave propagates almost constantly. In the +x direction, the amplified wave propagates with strong dispersion, and its maximum amplitude decreases during propagation. Fig. 12 shows the variations of the maximum amplitude ∗ ηmax /ζ0∗ as a function of vp∗ and τ ∗ for t ∗ = 100, 500 and 1000. The shape of the dependence does not change with time, except for small vp∗ and large τ ∗ because their T ∗ is over 100, or around vp∗ = 1 and small τ ∗ , where the amplitudes are greater than 1 (see the spatial profiles in Figs. 8, 9, 10). The resonance is preserved but these maxima decrease in time, which confirms the previous conclusion. 4.3. Comparison with numerical results

∗ Fig. 6. The maximum free surface amplitude, ηmax /ζ0∗ as function of vp∗ and τ ∗ for L∗ = 50 at t ∗ = T ∗ (colour in log scale). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

maximum free surface amplitude can reach six times the bottom deformation amplitude at vp∗ = 1 as shown in Fig. 5. Thus, when the rupture velocity is close to the wave velocity, the amplitude of the generated wave is much larger than the amplitude of the sea floor deformation at the beginning of the propagation phase. This resonance only exists for small values of τ ∗ . ∗ The deformation length is then varied, and η∗ /ζ0∗ and ηmax /ζ0∗ are compared to the case L∗ = 50 for vp∗ = 1 and τ ∗ = 0. On Fig. 7(a), wave profiles are plotted for different L∗ : while the wave in the negative direction keeps a constant amplitude, the wave in the positive direction increases with L∗ . This is consistent with the fact that a larger volume of water is displaced. The maximum amplitude of the deformation is also sensitive to L∗ (Fig. 7(b)). The maxima can be reasonably well fit with a power law of the form: ∗ ηmax /ζ0∗ = 0.414(L∗ )0.669 for L∗ ∈ [0; 5000]. (25) ∗ ∗ ∗ According to the fitting law (25), ηmax /ζ0 increases with L , as said

earlier.

4.2. Propagation stage The propagation of the generated wave is analysed as a function of vp∗ and τ ∗ , with fixed L∗ = 50. The spatial profiles of the free

In order to corroborate the proposed analytical solution, these results are compared to simulation results from a numerical model called Misthyc [22], which solves the Euler–Zakharov equations [23] for a homogeneous incompressible and inviscid fluid in an irrotational flow:

˜ ▽ η + w( ηt = − ▽ φ. ˜ 1 + (▽η)2 ), 1

1

2

2

˜ 2+ w φ˜ t = −g η − (▽φ) ˜ 2 (1 + (▽η)2 ), where φ˜ is the potential at the free surface and w ˜ is the vertical velocity at the free surface. These equations correspond to (5a)–(5d), rewriting the system using free surface variables. These equations are solved by using a 4th order Runge–Kutta (RK4) scheme in time, a spectral approach in the vertical direction, and high-order finite difference schemes (4th order) in the x direction. The spectral approach uses a base of Chebyshev polynomials of the first kind, and for the present simulations, the maximum polynomial order is 7. The choice of this model is justified by its capacity to represent well the propagation of nonlinear and dispersive waves (see [24]). Figs. 4 and 8–10 show comparisons between the present linear solution and the results from the linear version of Misthyc. On the overall, the curves globally match very well indicating a very close agreement between the theoretical solution and the results of the numerical simulations. At t ∗ = T ∗ , the relative errors (not shown here) are found to be less than 1% for every vp∗ as soon as τ ∗ > 0.

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

263

Fig. 7. Impact of the deformation length L∗ for vp∗ = 1 and τ ∗ = 0 at the end of the generation. 7(a) Spatial profiles of the free surface for different values of L∗ . 7(b) Red ∗ crosses represent the ηmax /ζ0∗ of the linear solution and the black line is the trend line 0.414(L∗ )0.669 .

(a) vp∗ = 0.5; τ ∗ = 0.

(b) vp∗ = 0.5; τ ∗ = 1.

(c) vp∗ = 0.5; τ ∗ = 2.

(d) vp∗ = 1; τ ∗ = 0.

(e) vp∗ = 1; τ ∗ = 1.

(f) vp∗ = 1; τ ∗ = 2.

(g) vp∗ = 2; τ ∗ = 0.

(h) vp∗ = 2; τ ∗ = 1.

(i) vp∗ = 2; τ ∗ = 2.

(j) vp∗ = 10; τ ∗ = 0.

(k) vp∗ = 10; τ ∗ = 1.

(l) vp∗ = 10; τ ∗ = 2.

(m) vp∗ = 50; τ ∗ = 0.

(n) vp∗ = 50; τ ∗ = 1.

(o) vp∗ = 50; τ ∗ = 2.

Fig. 8. Free surface profiles at t ∗ = 100 for L∗ = 50 and different values of vp∗ and τ ∗ . Same legend as Fig. 4.

264

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

(a) vp∗ = 0.5; τ ∗ = 0.

(b) vp∗ = 0.5; τ ∗ = 1.

(c) vp∗ = 0.5; τ ∗ = 2.

(d) vp∗ = 1; τ ∗ = 0.

(e) vp∗ = 1; τ ∗ = 1.

(f) vp∗ = 1; τ ∗ = 2.

(g) vp∗ = 2; τ ∗ = 0.

(h) vp∗ = 2; τ ∗ = 1.

(i) vp∗ = 2; τ ∗ = 2.

(j) vp∗ = 10; τ ∗ = 0.

(k) vp∗ = 10; τ ∗ = 1.

(l) vp∗ = 10; τ ∗ = 2.

(m) vp∗ = 50; τ ∗ = 0.

(n) vp∗ = 50; τ ∗ = 1.

(o) vp∗ = 50; τ ∗ = 2.

Fig. 9. Free surface profiles at t ∗ = 500 for L∗ = 50 and different values of vp∗ and τ ∗ . Same legend as Fig. 4.

The cases with τ ∗ = 0 (left columns on Fig. 4) correspond to an instantaneous vertical uplift speed, which can be represented only approximately in the numerical model. In practice a small but finite value of τ ∗ was used in Misthyc (τ ∗ = 0.001), leading to slightly larger errors in these cases. The results are thus consistent, which provides an independent validation of the solution developed in the previous section. 5. Discussion of the validity of the shallow water equations In this section, we discuss the validity of using the Shallow Water Equations (SWE) to model dynamic tsunami generation and propagation. The SWE are derived from the Euler equations by a depth-averaging procedure under the assumption of hydrostatic pressure. The latter is valid only for long waves, meaning that frequency dispersion is neglected. In the frame of a

seismic tsunami, it is commonly proposed that its characteristic wavelength is equivalent to the horizontal dimension of the source (see Section 2 and Wu [19]), thus it is generally much larger than the water depth, justifying the use of the SWE in such conditions. The question about dispersion effects during the tsunami process has already been asked and explored. Glimsdal et al. [25] proposed a ‘‘dispersion time’’ depending mainly on the initial water depth, the source width and the distance from the source region to the shore. This ‘‘dispersion time’’ gives a glimpse of the magnitude of the dispersive effects. In a similar aim, in this study, a parameter ε is defined as function of vp∗ and τ ∗ using the potential energy, as explained below, to measure the part of the energy lying in the dispersive range of k∗ poorly represented with the SWE. The long wave approximation, and thus the SWE, is generally considered to be valid for k∗ < 0.2 in the linear theory, which corresponds to wavelengths greater than approximately 30h. This

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

265

(a) vp∗ = 0.5; τ ∗ = 0.

(b) vp∗ = 0.5; τ ∗ = 1.

(c) vp∗ = 0.5; τ ∗ = 2.

(d) vp∗ = 1; τ ∗ = 0.

(e) vp∗ = 1; τ ∗ = 1.

(f) vp∗ = 1; τ ∗ = 2.

(g) vp∗ = 2; τ ∗ = 0.

(h) vp∗ = 2; τ ∗ = 1.

(i) vp∗ = 2; τ ∗ = 2.

(j) vp∗ = 10; τ ∗ = 0.

(k) vp∗ = 10; τ ∗ = 1.

(l) vp∗ = 10; τ ∗ = 2.

(m) vp∗ = 50; τ ∗ = 0.

(n) vp∗ = 50; τ ∗ = 1.

(o) vp∗ = 50; τ ∗ = 2.

Fig. 10. Free surface profiles at t ∗ = 1000 for L∗ = 50 and different values of vp∗ and τ ∗ . Same legend as Fig. 4.

threshold corresponds to the limit for √which the relative error gh and the linear celerity between the linear long wave celerity √ (g /k) tanh(kh) is less than 1%. On Fig. 13, the range of validity of the SWE approximation based on this threshold is illustrated with the potential energy field for vp∗ = 50 and τ ∗ = 0 at t ∗ = T ∗ . The solid blue line is the modulus of the free surface Fourier transform η˜ ∗ of the linear solution. Outside of the grey striped zone that corresponds to the zone of validity of the Linear SWE (LSWE), the long wave approximation is no longer valid, and the energy is not properly modelled by the SWE. To estimate the impact of the SWE approximation on the present tsunami waves, we study the total potential energy. The amount of energy that lies outside the range of validity of the SWE is simply called ‘‘residual’’ energy hereafter. The ratio ε between the residual energy El and the total energy Et is estimated. The total energy is calculated numerically by integrating η˜ ∗2 (k∗ ) over

infinity, and the ‘‘residual’’ energy is calculated by integrating η˜ ∗2 (k∗ ) over ] − ∞; −0.2] ∪ [0.2; ∞[ (i.e. the part of the spectrum not properly handled by the SWE):

ε(t ) = ∗

El (t ∗ ) Et (t ∗ )

˜ |k∗ |≥0.2 (η



=  +∞ −∞

) dk∗

∗ 2

(η˜ ∗ )2 dk∗

.

(26)

In Fig. 14, the ratio ε(t ∗ ) is presented as a function of vp∗ and τ ∗ at t ∗ = T ∗ . Here again, the resonance around vp∗ = 1 and small τ ∗ is visible. The energy ratio ε for vp∗ = 1 and τ ∗ = 0 reaches 0.61, which means that more than half of the potential energy is not properly accounted for at the end of the deformation if a shallow water model is used. For larger τ ∗ and vp∗ , the ratio is typically less than 0.1, indicating that shallow water models are more appropriate. Globally, the energy ratio ε decreases when τ ∗ increases. Profiles of ε as function of time for vp∗ = 0.5, 1, 2, 10 and 50 and ∗ τ = 0, 1, 2 are shown in Fig. 15. For the sake of clarity, the profiles

266

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

Fig. 11. Evolution in space and time of the free surface deformation η∗ /ζ0∗ for

vp∗ = 1 and τ ∗ = 0. The colour represents the wave amplitude.

for the different τ ∗ are plotted on different graphs, and each colour line represents a rupture velocity vp∗ . In all cases, three phases can be identified in the evolution of ε : an initial steady phase is followed by a decrease and then by a second steady state. The first phase, with t ∗ < L∗ , is nearly independent of τ ∗ and slightly sensitive to vp∗ : for vp∗ = 0.5, 1, 2, the initial value of ε is above 0.5, and decreases when vp∗ increases. During this phase, ε remains high (>0.1) thus SWE models are not recommended for this range of parameters. The transition phase, around t ∗ = L∗ , shows a rapid decrease of ε . Its value oscillates before it reaches an asymptotic value in the third phase (t ∗ > L∗ ). Except for vp∗ = 1, this limit value is the same for all vp∗ but varies with τ ∗ . When τ ∗ increases, this limit decreases: for τ ∗ = 0, ε is close to 0.05, for τ ∗ = 1, ε is around 10−4 ; and for τ ∗ = 2, ε reaches 5 10−5 . During this phase, and for these values of vp∗ , SWE models are appropriate. For vp∗ = 1 and τ ∗ = 1, 2, ε follows the same pattern but with a higher asymptotic value in the third phase than for the other rupture velocities. The evolution of ε is different when τ ∗ = 0. In this case, ε remains above 0.5 and seems to increase slightly. The resonance is still observed. Thus, when vp∗ = 1 and τ ∗ is very small, SWE models should not be used. 6. Numerical simulation of nonlinear propagation In this section, we investigate how the resonance observed in the linear model behaves when nonlinear effects are taken into account. The previous problem is now solved with the nonlinear version of Misthyc (see sub-Section 4.3). Initially the deformation is nearly linear with ζ0∗ = 0.001 (small deformation of the sea floor compared to the water depth). Fig. 16 shows the resulting ∗ distribution of ηmax /ζ0∗ for the same range of parameters (vp∗ and τ ∗ ) considered in Section 4. Then nonlinear effects are introduced by increasing ζ0∗ from 0.001 to 0.1. Indeed, if ζ0 increases ηmax

(a) t ∗ = 100.

(b) t ∗ = 500.

Fig. 13. The blue line represents the modulus of η˜ ∗ for vp∗ = 50 and τ ∗ = 0 at t ∗ = T ∗ . The grey striped zone shows the domain of validity of the long wave approximation. Waves with wavenumbers outside of this zone are not taken into account properly in the SWE model.

Fig. 14. Energy ratio of the residual energy El over the total energy Et , ε = E l , at t t ∗ = T ∗ as a function of vp∗ and τ ∗ . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) E

increases too because the volume of displaced water is bigger. ∗ Thus the nonlinearity, that increases with ηmax , also increases ∗ with ζ0 . Here, the solution is computed only for vp∗ = 1 and τ ∗ = 0. Fig. 17 shows the deformation of the free surface for ∗ different values of ζ0∗ at t ∗ = T ∗ . We can see that ηmax varies −3 −1 from 5.10 to 6.10 , increasing nonlinear effects. The numerical results for ζ0∗ = 0.001, 0.005, 0.01 are close to the linear solution. The general wave form is preserved, and the maximum amplitude increases slightly with ζ0∗ . For ζ0∗ = 0.05 and 0.1, the propagation is faster. The frequency dispersion is more pronounced, and the first generated wave seems to have split already in two for ζ0∗ = 0.1, leading to a lower maximum amplitude than observed in the other cases. In reality, this kind of deformation (ζ0∗ = 0.05 and 0.1) is very large for seismic tsunami generation: the ground deformation is usually on the order of 10 m for a water depth of 4000 m, leading to ζ0∗ on the order of 10−3 .

(c) t ∗ = 1000.

∗ Fig. 12. Maximum free surface amplitude ηmax /ζ0∗ as a function of vp∗ and τ ∗ for L∗ = 50 at different times t ∗ . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

Fig. 15. Time profiles of the energy ratio ε =

El

267

for τ ∗ = 0, 1, 2 (Fig. 15(a)–(c), respectively) for a range of rupture velocities vp∗ . The horizontal black line represents

√ Et ε = 0.5. The vertical black line is at t = (L/ gh)∗ . ∗

∗ Fig. 16. The maximum free surface amplitude (ηmax /ζ0∗ ) calculated with the nonlinear numerical model as function of vp∗ and τ ∗ for L∗ = 50 at t ∗ = T ∗ (colour in log scale). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

∗ Fig. 18. Temporal profiles of ηmax /ζ0∗ for vp∗ = 1, τ ∗ = 0 and L∗ = 50. The solid black line is the theoretical linear solution. The coloured dashed lines are nonlinear numerical results from Misthyc for ζ0∗ = 0.001, 0.005, 0.01 and 0.1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 17. Deformation of the free surface at t ∗ = T ∗ and for vp∗ = 1 and τ ∗ = 0: comparison of the linear solution (black line) and nonlinear numerical results. The coloured dashed lines represent the dimensionless nonlinear numerical results for different increasing non-linearity levels: ζ0∗ = 0.001, 0.005, 0.01, 0.1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 19. Free surface profiles from nonlinear numerical results Misthyc, for ζ0∗ = 0.1 for vp∗ = 0.8, 0.9, 1, 1.1, 1.26 at the end of the ground motion. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Temporal profiles of the maximum free surface amplitude are drawn in Fig. 18 for a range of ζ0∗ (from 0.001 to 0.1). Until the ∗ end of the ground deformation (t ∗ = T ∗ = 50 here), ηmax /ζ0∗ increases similarly, reaching an almost identical maximum for all ∗ the ground deformation amplitudes ζ0∗ (ηmax ∼ 6ζ0∗ ). Thus, the resonance remains when nonlinear effects are included. However, ∗ the decrease of ηmax is different. The slope is smaller for larger ζ0∗ . ∗ ∗ In addition, for ζ0 > 0.05, ηmax remains nearly constant in time. With nonlinearity, it can be expected that the resonance is shifted towards higher values of vp∗ . Indeed, the linear theory showed that the resonance occurs when vp is equal to the long wave celerity. In weakly √ nonlinear conditions, the latter can be approximated by c = gh(1 + ηmax /h)(Boussinesq, [26]). The peak of resonance ∗ . As the value of η should occur for vp∗ = 1 + ηmax max is unknown, preliminary tests are done for ζ0∗ = 0.1 and vp∗ ∈ [0.8; 1.26] at the end of the ground motion. The free surface profiles are plotted

Fig. 19. The amplitude of the motion with vp∗ = 1.1 is larger than the one with vp∗ = 1, as expected. This subject deserves a deeper analysis in a future work. For the sake of simplicity, we keep vp∗ = 1 thereafter. The free surface deformation after the end of the ground deformation: t ∗ = 100 and t ∗ = 500, is shown in Fig. 20. For small ζ0∗ , the shape of the waves generated by the nonlinear numerical model is similar to the linear solution, even though the propagation speed is slightly higher for ζ0∗ = 0.005 and 0.01. However, for ζ0∗ = 0.05 and 0.1, a different pattern is observed. Solitary waves ∗ appear, which explains the constant value of ηmax in time after some duration (Fig. 18) as solitary wave solutions do not change amplitude during propagation. A more in-depth analysis of the transition between the two regimes illustrated in Figs. 18 and 20, and the condition of appearance of solitary wave(s) will be a topic of future research. Here, a first comparison is done with numerical results from the algorithm of Dutykh and Clamond [27]. The shapes

268

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

(a) t ∗ = 100.

(b) t ∗ = 500.

Fig. 20. Comparison of the deformation of the free surface at t ∗ = 100 and 500 (20(a) and (b) respectively) for vp∗ = 1 and τ ∗ = 0 of the linear solution and the nonlinear numerical results. The thick black line is the theoretical linear solution and the coloured lines represent the nonlinear numerical results of Misthyc for different initial deformations: ζ0∗ = 0.001, 0.005, 0.01, 0.1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(a) ζ0∗ = 0.05.

(b) ζ0∗ = 0.1.

Fig. 21. Comparison between the shape of the generated wave computed with the nonlinear numerical code Misthyc (black full curve) and the shape of the solitary wave from the algorithm of Dutykh and Clamond [27] (dashed green curves) at t ∗ = 993. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

∗ of solitary wave created for a same ηmax are superimposed in Fig. 21 on the free surface profile of the generated wave with a deformation of ζ0∗ = 0.05, 0.1 at t ∗ = 993. Qualitatively, the shapes match very well, which confirms that these leading waves are of solitary type.

Table 2 Characteristics of the March 25 1947 event, which occurred in the Offshore Poverty Bay (NZ). Date Epicentre location

03/25/1947 38.92SS, 178.24E

Fault characteristics

7. Application to a real tsunami event

Length Width Depth of the fault Dip angle Rise time (fault’s rupture) Water depth near the source Rise time tr Rupture velocity vp Deformation of the ground ζ0

55 km 50 km 5–12 km 8° 10–20 s ∼1500 m (between 100 and 3000 m) ∼90 s (between 60 and 120 s) ∼175 m s−1 (between 150 and 300 m s−1 ) ∼1 m

‘‘Tsunami earthquake’’ is a category of events proposed by Kanamori [28] that generate larger tsunami than suggested by the seismic magnitude. They are generally caused by slow rupture velocities of 1 km/s or less. Here, one such event that occurred March 25 1947 near New Zealand is chosen as an application of the linear model (see [29,30]). Fig. 22 shows the domain of application and the location of the event. The fault rupture is close to the Hikurangi subduction margin and the tsunami was observed along the coast between Tokomaru Bay and Mahia Peninsula. The fault characteristics are summarised in Table 2. With a water depth near the source of 1500 m and a sea floor deformation ζ0 of 1 m, ζ0 ∼ 6 10−4 ≪ 1, thus the generation can be considered as linear. h If we assume that vp ∈ [150; 300] m s−1 with a mean value of 175 m s−1 , then vp∗ ∈ [0.87; 9.57] with a mean value of vp∗ = 1.44.

vp∗ is close to 1 and τ ∗ < 1, thus this event is near the resonance

3000 m. The largest vp∗ is reached for vp = 300 m s−1 and h = 100 m. Likewise, if tr ∈ [60; 120] s with a mean value of tr = 90 s, we have τ ∗ ∈ [0.034; 0.36] with a mean value of τ ∗ = 0.19. The smallest τ ∗ is reached for tr = 60 s and a depth of 100 m. The largest τ ∗ is reached for tr = 120 s and h = 3000 m. These numerical values were provided by William Power (personal communication). These parameter values are shown in the plot of the maximum amplitude (Fig. 23) for a fault length of L∗ = 35. The mean value of

In this article, a semi-analytical expression of the free surface elevation induced by a simple seismic-like motion of the sea floor has been developed using linear theory as a function of τ ∗ (rise time in the vertical direction) and vp∗ (rupture velocity in the horizontal direction), which are the dimensionless temporal parameters of a schematic fault retained after doing a dimensional analysis. The derived solution allowed studying two aspects: the ∗ maximum amplitude ηmax reached by the first wave (as a function of the amplitude of the bottom upthrust ζ0∗ ) and the validity of the long-wave assumption of the shallow water equations (SWE).

The smallest vp∗ is reached for vp = 150 m s−1 and a depth of

condition previously defined. Therefore, one can expect that the tsunami will be larger than predicted when using models that do not take into account kinematic deformations. This may explain the large water waves observed during this event, as suggested by Bell et al. [29]. 8. Conclusions

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

269

Fig. 22. Locations of the March and May 1947 Gisborne earthquakes, estimated tsunami run-up, main features of the plate boundary in New Zealand (inset), and location of seismographs of the 1947 New Zealand seismograph network (inset) are shown. Source: The figure is reprinted from the paper of Downes and Stirling [30].

∗ An analysis of the non-dimensional maximum amplitude ηmax /ζ0∗ showed that:

1. when τ ∗ is large (i.e. slow vertical motion), vp∗ does not impact the maximum amplitude. On the contrary, when τ ∗ is small (i.e. rapid vertical motion), the maximum amplitude is strongly influenced by vp∗ . 2. A resonance phenomenon exists near vp∗ = 1 and τ ∗ = 0, this result coincides with the one from Todorovska and Trifunac [4]. The free surface amplitude is amplified compared to the deformation of the ground, by a factor depending on L∗ . The maximum amplitude could be approximated by a power law of L∗ that was fit to the simulated data (see Eq. (25)). 3. The impact of the resonance remains during the propagation phase even if ηmax slightly decreases. The resonance phenomenon remains when nonlinear effects are included as verified by using a fully nonlinear and dispersive numerical model. Furthermore, if ζ0∗ is large enough, solitary waves appear after a certain time/distance and then remain stable in time. To evaluate the impact of timescales on the validity of the SWE, an energy ratio ε was defined evaluating the fraction of potential energy that is not properly modelled by the SWE (as this energy is associated with wavelengths shorter than the non-dispersive limit L ≈ 30h) over the total potential energy. This analysis shows that: 1. there is still a resonance around vp∗ = 1 and τ ∗ = 0. However, for these parameter values and L∗ = 50, the energy lying outside the shallow-water range is more than half of the total energy. 2. During the propagation phase, ε reaches an asymptotic state that depends on τ ∗ : the larger τ ∗ , the smaller ε (except when vp∗ = 1 and τ ∗ = 0, ε remains large). This analysis identifies the range of parameter values for which models based on the SWE should or should not be used.

Fig. 23. Maximum amplitude of the free surface for L∗ = 35 at t ∗ = T ∗ . The black segments represent the ranges of values of the temporal parameters τ ∗ and vp∗ for the New Zealand 1947 event (see Section 7). The middle point is the mean value of both parameters. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

To illustrate possible applications of this work, the tsunami event of 1947 close to New Zealand was analysed. Based on estimates of the physical parameters of this event, it was identified as being close to the resonance condition, and thus the free surface waves generated by the earthquake were likely to have been much larger than the amplitude of the ground motion. Acknowledgements This work was funded by the French project ‘‘Tsunamis in the Atlantic and the English ChanNel: Definition of the Effects through numerical Modelling’’ (TANDEM, ref.: ANR-11-RSNR-0023-01). The authors would like to thank Marissa Yates and Cécile Raoult

270

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

from the Saint-Venant Hydraulics Laboratory (France) for their help with the program Misthyc, as well as Anne Lemoine from the Bureau de Recherches Géologiques et Minières (BRGM, France), and Rebecca Bell and William Power from GNS Science, Lower Hutt (New Zealand) for their useful discussions and advice. Appendix A. Fourier transform of the free surface

Thus: if B = A ∗ e−ikvp t (1 − H (t − vL )) p

 ωr B = 2 ωr − ω2

  1   ω2 − k2 v 2 [ikvp cos ωt p 

In this appendix, the calculations from Eqs. (21) to (22) are shown. An inverse Laplace transform is applied. The latter satisfies:

+ ω sin ωt − ikvp e−ikvp t ]

L {f ∗ g } = L {f } L {g }

(A.1)

−

L−1 {f } ∗ L−1 {g } = L−1 {fg }

(A.2)

= cos ωtH (t )

(A.4)

= δ (t ) + δ (t − tr )

(A.5)

= sin (ωr t ) H (t )

(A.6)

(A.7)

In order to calculate C = B ∗ (δ(t ) + δ(t − tr )), we use: f ∗ δ(t − a) = f (t − a)H (t − a): 

where: t



f (τ )g (t − τ )dτ

f ∗g =

(A.3)

0

and

L

−1



s2 + ω 2

 −1

L

L −1

L−1



s

1+e

 −str

ωr





s2 + ωr2

   1 − e− vLp (ikvp +s)  ikvp + s







= e−ikvp t 1 − H t −



L



vp

,

      ωr 1  C = 2 ikvp cos ωt + ω sin ωt − ikvp e−ikvp t 2 2 2 2 ωr − ω   ω − k vp  

thus:

η( ˜ k, t ) =

ζ0

ωr

2 cosh kh

cos ωtH (t ) ∗ (δ (t ) + δ (t − tr ))

   L . ∗ sin (ωr t ) H (t ) ∗ e−ikvp t 1 − H t − vp

(A.8)

A = cos ωtH (t ) ∗ sin (ωr t ) H (t )

ωr (cos ωt − cos ωr t ) ωr2 − ω2    L A ∗ e−ikvp t 1 − H t − vp   cos ωt ∗ e−ikvp t     − cos ωt ∗ e−ikvp t H t − L   ωr vp    = 2    ωr − ω2 − cos ωr t ∗ e−ikvp t      L  −ikvp t + cos ωr t ∗ e H t− vp =

(A.9)

(A.10)

also: cos ωt ∗ e−ikvp t   1 = 2 ikvp cos ωt + ω sin ωt − ikvp e−ikvp t ω − k2 vp2

(A.11)

and:



cos ωt ∗ e−ikvp t H



e−iLk



t−

L

=

H t − vL p



ω2 − k2 vp2    L L ikvp cos ω t − + ω sin ω t − vp vp  . (A.12)

× −ikvp e−ikvp t

vp 



  1 ikvp cos ωr t + ωr sin ωr t − ikvp e−ikvp t ωr2 − k2 vp2    1 − −ikvp cos ωr t   ωr2 − k2 vp2      −ωr sin ωr t − ikvp e−ikvp (t −tr )  + H (t − tr )     1  +   ω2 − k2 v 2 ikvp cos ω(t − tr ) p  −ikvp (t −tr ) +ω sin ω(t − tr ) − ikvp e       L e−iLk ikvp cos ωr t −   1    vp   2  L  ωr − k2 vp2  +ωr sin ωr t − − ikvp e−ikvp t     v L  p       +H t − L  vp    e−iLk ikvp cos ω t −    v 1 p −      ω 2 − k2 v 2   L p +ω sin ω t − − ikvp e−ikvp t vp   L + H t − tr − vp       L −iLk e − ik v cos ω t − p r   1     vp  ω 2 − k2 v 2    L − ik v ( t − t )  r p r p −ω sin ω t − − ik v e r r p   v   p        L   . (A.14) − iLk ×  e ikvp cos ω t − tr −    v   p       1   −  L   ω 2 − k2 v 2   +ω sin ω t − t −   r   p vp   −

We then define A as:



1

[ikvp cos ωr t + ωr sin ωr t − ikvp e−ikvp t ] ωr2 − k2 vp2      H t − vL  L p −iLk − 2 e ik v cos ω t − p ω − k2 vp2 vp    L + ω sin ω t − − ikvp e−ikvp t vp      L −iLk   e ikvp cos ωr t − v   H t − vL    p  p   L (A.13) + 2 .    ωr − k2 vp2  +ωr sin ωr t −   vp −ikvp e−ikvp t

−ikvp e−ikvp (t −tr )

Combining (A.8) and (A.14) gives (22) and (24) in the dimensionless space.

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

Appendix B. Singularities

271

For case (ii), the same process is used to obtain:

 In this appendix, the singularities of η˜ are treated. We can identify three critical cases where η˜ is discontinuous: (i) k2 vp2 → ωr2 , (ii) k2 vp2 → ω2 , (iii) ω → ωr2 . As stated in paragraph 3.3, to maintain a continuous function, the limit of η˜ is taken when |k − kcritical |h < 10−6 . For case (i), D is defined as: D = ikvp cos ωr t + ωr sin ωr t − ikvp e−ikvp t . We have: D

ω − 2 r

v

k2 p2

→

1 2ωr

(ωr te−iωr t + sin ωr t ).

Thus, when k2 vp2 → ωr2 :

η˜ (k, t ) →

ζmax

vp

ωr2 2 cosh kh ωr2 − ω2

       1  ikvp cos ωt + ω sin ωt − ikvp e−ikvp t × 2 2 2  ω − k vp   

−

1 2ωr

(ωr te−iωr t + sin ωr t ) 

1

ikvp cos ω(t − tr )





 ω2 − k2 vp2       − ikvp (t −tr ) + H (t − tr )  +ω sin ω(t − tr ) − ikvp e     1  −iωr (t −tr ) (ωr (t − tr )e + sin ωr (t − tr )) − 2ωr   L +H t − vp        L −iLk e ik v cos ω t − p   1       vp  − 2  L  ω − k2 vp2   −ikvp t +ω sin ω t − − ik v e   p ×  vp  −i L ω          e vp r  L L −iωr t − vLp − ωr t − e + sin ωr t − 2ωr vp vp   L + H t − tr − vp       L e−iLk ikvp cos ω t − tr −  vp          1   −  L    ω2 − k2 v 2  +ω sin ω t − tr −     p   vp      − ik v ( t − t ) ×   . (B.1) −ikvp e p r       L   L −iωr t −tr − vp  −i L ωr ωr t − tr −  e  e vp    v p −         2 ωr  L + sin ωr t − tr − vp

    1 2 ζmax vp ωr  η˜ (k, t ) →  (ωte−iωt + sin ωt ) 2 cosh kh ωr2 − ω2  2ω    −

1



ωr2 − k2 vp2 

ikvp cos ωr t + ωr sin ωr t − ikvp e−ikvp t 1



(ω(t − tr )e−iω(t −tr ) + sin ω(t − tr ))



 2ω       1 + H ( t − tr )  +  ik v cos ω t + ω sin ω t p r r r  ωr2 − k2 vp2   +ikvp e−ikvp (t −tr )  −i L ω       L e vp −iω t − vLp ω t − − e   2ω  vp      L  + sin ω t −      vp  L       +H t −   L  vp  e−iLk ikvp cos ωr t −    vp     1   +   L     2 2 2  ωr − k vp  +ωr sin ωr t −  vp −ikvp e−ikvp t   L + H t − tr − vp        L −iω t −tr − vLp L ω t − t − e −i vp ω r    vp  e      −    L 2ω    + sin ω t − tr −   vp         . L (B.2) ×   e−iLk −ikvp cos ωr t −     v p      1   +   L   −ω sin ω t −  ωr2 − k2 vp2  r r   vp −ikvp e−ikvp (t −tr ) Moreover, if k → 0 and ω → 0:

η˜ (k, t ) →

ζmax

vp

ωr2 2 cosh kh ω − ω2  2 r

   1 × t − ω2 − k2 v 2 ikvp cos ωr t r p   + ωr sin ωr t − ikvp e−ikvp t    1 t − tr + ik v cos ω t p r   ωr2 − k2 vp2 + H (t − tr )    −ikvp (t −tr ) +ωr sin ωr t + ikvp e   L +H t − vp     L 1 − t− + 2   vp ωr −  k2 vp2       L  × e−iLk ikvp cos ωr t − ×   v p       L +ωr sin ωr t − − ikvp e−ikvp t vp

272

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273

 + H t − tr −

L



vp     L 1 − t − tr − + 2 2 2   ωr  − k vp   vp      L − iLk . (B.3) × e − ikvp cos ωr t − ×    v p        L − ikvp e−ikvp (t −tr ) −ωr sin ωr t − vp For case (iii), we define: 1 1 ωr2 ω − ω2 ω2 − k2 vp2 ωr2 − k2 vp2   ikvp e−ikvp t (ω2 − ωr2 ) +ikvp ωr2 cos ωt − ω2 cos ωr t     2 2 ωt )  ×  +k 2vp (cos ωr t − cos .  +ωω sin ωt − ωr ω2 sin ωr t  r +k2 vp2 (ωr sin ωr t − ω sin ωt )

F =

2 r

When ω → ωr , we obtain: F →

ωr2 (ωr2 − k2 vp2 )2   −ikvp e−ikvp t   2 2 2   +ikv cos ω t − (k vp − ωr )t sin ω t  p r r   2ωr       × . sin ωr t t 2   +ωr − cos ωr t   2ωr 2       t sin ω t r +k2 vp2 + cos ωr t 2ωr 2

Thus:

η˜ (k, t ) →

ζmax

vp

ωr2 2 cosh kh (ω − k2 vp2 )2 2 r

         (k2 vp2 − ωr2 )t  −ikvp t + ikvp cos ωr t − ×  − ikvp e sin ωr t  2ωr      + ωr2



sin ωr t 2ωr

−

t 2



cos ωr t + k2 vp2



sin ωr t 2ωr

−ikvp e−ikvp (t −tr )



+

t 2

cos ωr t





2 2 2   +ikv cos ω (t − t ) − (k vp − ωr )(t − tr )  p r r   2ωr          × sin ω ( t − t ) r r          sin ω ( t − t ) r r   2 + H ( t − tr )  +ωr  2ωr      ( t − tr )   − cos ω ( t − t ) r r   2        sin ω ( t − t ) r r 2 2   +k vp   2 ω r    ( t − tr ) cos ωr (t − tr ) +



2

  L −H t − vp 

−ikvp e 

  −ikvp t − vLp



        L    + ik v cos ω t − p r   v p          L 2 2 2   ( k v − ω ) t −   p r vp L    − sin ωr t −   2 ω v  r p ×          L L     t − sin ω t − r   vp vp L  +ω2   − cos ω t − r r   2 ω 2 v r p            L L     t− v L  2 2  sin ωr t − vp p  + cos ωr t − +k v p 2ωr 2 vp   L − H t − tr − vp     −ikv t −t − L −ikvp e p r vp          L    cos ω t − t − r r      vp             (k2 vp2 − ωr2 ) t − tr − vLp     +ikvp       −       2 ω r        L     × sin ωr t − tr −   vp        L sin ω t − t − r r   v p 2   +ωr  . (B.4) ×  2 ω r              t − tr − vL  L p   −  cos ωr t − tr −   2 vp            L sin ωr t − tr − v   p 2 2   + k v   p 2 ωr          L      t − tr − v L p    +  cos ωr t − tr − 2 vp

References [1] R.D. Braddock, P.V.D. Driessche, G.W. Peady, Tsunami generation, J. Fluid Mech. 59 (1973) 817–828. [2] P.C. Sabatier, On water waves produced by ground motions, J. Fluid Mech. 126 (1983) 27–58. [3] S.-J. Lee, G.T. Yates, T.Y. Wu, Experiments and analyses of upstream-advancing solitary waves generated by moving disturbances, J. Fluid Mech. 199 (1989) 569–593. [4] M.I. Todorovska, M.D. Trifunac, Generation of tsunamis by a slowly spreading uplift of the sea floor, Soil Dyn. Earthq. Eng. 21 (2001) 151–167. [5] B. Levin, M. Nosov, Physics of Tsunamis, Springer Science, 2009. [6] C. Ji, D.J. Wald, D.V. Helmberger, Source description of the 1999 Hector Mine, California, earthquake, part I: Wavelet domain inversion theory and resolution analysis, Bull. Seismol. Soc. Amer. 92 (4) (2002) 1192–1207. [7] Y. Okada, Internal deformation due to shear and tensile faults in a half-space, Bull. Seismol. Soc. Amer. 82 (2) (1992) 1018–1040. [8] M. Nosov, Tsunami waves of seismic origin: The modern state of knowledge, Izv. Atmos. Ocean. Phys. 50 (5) (2014) 474–484. [9] J. Hammack, A note on tsunamis: their generation and propagation in an ocean of uniform depth, J. Fluid Mech. 60 (1973) 769–799. [10] T. Jamin, L. Gordillo, G. Ruiz-Chavarría, M. Berhanu, E. Falcon, Experiments on generation of surface waves by an underwater moving bottom, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 471 (2178) (2015). [11] T.S. Stefanakis, F. Dias, C.E. Synolakis, Tsunami generation above a sill, Pure Appl. Geophys. 172 (3) (2015) 985–1002. [12] D. Dutykh, Mathematical modelling of tsunami waves (Ph.d. thesis), École normale supérieure de Cachan - ENS Cachan, (France), 2007.

M. Le Gal et al. / European Journal of Mechanics B/Fluids 65 (2017) 257–273 [13] D. Dutykh, F. Dias, Y. Kervella, Linear theory of wave generation by a moving bottom, C. R. Math. 343 (7) (2006) 499–504. [14] D. Dutykh, D. Mitsotakis, X. Gardeil, F. Dias, On the use of the finite fault solution for tsunami generation problems, Theor. Comput. Fluid Dyn. 27 (2013) 177–199. [15] A. Vaschy, Sur les lois de similitude en physique, Ann. Télégr. 19 (1892) 25–28. (in French). [16] E. Buckingham, On physically similar systems; Illustrations of the use of dimensional equations, Phys. Rev. 4 (1914) 345–376. [17] J. Johnson, K. Satake, Estimation of seismic moment and slip distribution of the April 1, 1946, Aleutian tsunami earthquake, J. Geophys. Res. 102 (11) (1997) 11, 765–11, 774. [18] Y. Fujii, K. Satake, Tsunami source of the 2004 Sumatra–Andaman earthquake inferred from tide gauge and satellite data, Bull. Seismol. Soc. Amer. 97 (1A) (2007) S192–S207. [19] T.Y. Wu, Long waves in ocean and coastal waters, J. Eng. Mech. 107 (EM3) (1981) 501–522. [20] C.C. Mei, M. Stiassnie, D.K.-P. Yue, Theory and Applications of Ocean Surface Waves, in: Advanced Series on Ocean Engineering, vol. 23, World Scientific Publishing, 2005. [21] K. Satake, Mechanism of the 1992 Nicaragua tsunami earthquake, Geophys. Res. Lett. 21 (23) (1994) 2519–2522.

273

[22] M.L. Yates, M. Benoit, Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves, Internat. J. Numer. Methods Fluids 77 (10) (2015) 616–640. [23] V.E. Zakharov, Stability of periodic waves of finite amplitude on the surface of a deep fluid, J. Appl. Mech. Tech. Phys. 9 (1968) 190–194. [24] C. Raoult, M. Benoit, M.L. Yates, Validation of a fully nonlinear and dispersive wave model with laboratory non-breaking experiments, Coast. Eng. 114 (2016) 194–207. [25] S. Glimsdal, G.K. Pedersen, C.B. Harbitz, F. Løvholt, Dispersion of tsunamis: does it really matter? Nat. Hazards Earth Syst. Sci. 13 (6) (2013) 1507–1526. [26] J. Boussinesq, Théorie de l’intumescence liquide appelée onde solitaire ou de translation se propageant dans un canal rectangulaire, C. R. Acad. Sci., Paris 72 (1871) 755–759. (in French). [27] D. Dutykh, D. Clamond, Efficient computation of steady solitary gravity waves, Wave Motion 51 (1) (2014) 86–99. [28] H. Kanamori, Mechanism of tsunami earthquakes, Phys. Earth Planet. Inter. 6 (1972) 346–359. [29] R. Bell, C. Holden, W. Power, X. Wang, G. Downes, Hikurangi margin tsunami earthquake generated by slow seismic rupture over a subducted seamount, Earth Planet. Sci. Lett. 397 (2014) 1–9. [30] G. Downes, M.W. Stirling, Groundwork for development of a probabilistic tsunami hazard model for New Zealand, in: International Tsunami Symposium 2001, 2001, pp. 293–301.