Waves generated by a moving bottom

Get Complete Project Material File(s) Now! »

Waves at the surface of a liquid can be generated by various mechanisms: wind blowing on the free surface, wavemaker, moving disturbance on the bottom or the surface, or even inside the liquid, fall of an object into the liquid, liquid inside a moving container, etc. In this chapter, we concentrate on the case where the waves are created by a given motion of the bottom. One example is the generation of tsunamis by a sudden seafloor deformation.
There are diﬀerent natural phenomena that can lead to a tsunami. For example, one can mention submarine slumps, slides, volcanic explosions, etc. In this chapter we use a submarine faulting generation mechanism as tsunami source. The resulting waves have some well-known features. For example, characteristic wavelengths are large and wave amplitudes are small compared with water depth.
Two factors are usually necessary for an accurate modelling of tsunamis: information on the magnitude and distribution of the displacements caused by the earthquake, and a model of surface gravity waves generation resulting from this motion of the seafloor. Most studies of tsunami generation assume that the initial free-surface deformation is equal to the vertical displacement of the ocean bottom. The details of wave motion are neglected during the time that the source operates. While this is often justified because the earthquake rupture occurs very rapidly, there are some specific cases where the time scale of the bottom deformation may become an important factor. This was emphasized for example by Trifunac and Todorovska [TT01], who considered the generation of tsunamis by a slowly spreading uplift of the seafloor and were able to explain some observations. During the 26 December 2004 Sumatra-Andaman event, there was in the northern extent of the source a relatively slow faulting motion that led to significant vertical bottom motion but left little record in the seismic data. It is interesting to point out that it is the inversion of tide-gauge data from Paradip, the northernmost of the Indian east-coast stations, that led Neetu et al. [NSS+05] to conclude that the source length was greater by roughly 30% than the initial estimate of Lay et al. [LKA+05]. Incidentally, the generation time is also longer for landslide tsunamis.
Our study is restricted to the water region where the incompressible Euler equations for potential flow can be linearized. The wave propagation away from the source can be investigated by shallow water models which may or may not take into account nonlinear eﬀects and frequency dispersion. Such models include the Korteweg-de Vries equation [KdV95] for unidirectional propagation, nonlinear shallow-water equations and Boussinesq-type models [Bou71b, Per66, BBM72].
Several authors have modeled the incompressible fluid layer as a special case of an elastic medium [Pod68, Kaj63, Gus72, AG73, Gus76]. In our opinion it may be convenient to model the liquid by an elastic material from a mathematical point of view, but it is questionable from a physical point of view. The crust was modeled as an elastic isotropic half-space. This assumption will also be adopted in the present study.
The problem of tsunami generation has been considered by a number of authors: see for example [Car71, vdD72, BvdDP73]. The models discussed in these papers lack flexibility in terms of modelling the source due to the earthquake. The present study provides some extensions. A good review on the subject is [Sab86].
Here we essentially follow the framework proposed by Hammack [Ham73] and others. The tsunami generation problem is reduced to a Cauchy-Poisson boundary value problem in a region of constant depth. The main extensions given in the present work consist in three-dimensional modelling and more realistic source models. This approach was followed recently in [TT01, THT02], where the mathematical model was the same as in [Ham73] but the source was diﬀerent.
Most analytical studies of linearized wave motion use integral transform methods. The complexity of the integral solutions forced many authors [Kaj63, Kel63] to use asymptotic methods such as the method of stationary phase to estimate the far-field behaviour of the solutions. In the present study we have also obtained asymptotic formulas for integral solutions. They are useful from a qualitative point of view, but in practice it is better to use numerical integration formulas [Fil28] that take into account the oscillatory nature of the integrals. All the numerical results presented in this section were obtained in this manner.
One should use asymptotic solutions with caution since they approximate exact solu-tions of the linearized problem. The relative importance of linear and nonlinear eﬀects can be measured by the Stokes (or Ursell) number [Urs53]: U := a/h = a , (kh)2 k2h3 where k is a wave number, a a typical wave amplitude and h the water depth. For U ≫ 1, the nonlinear eﬀects control wave propagation and only nonlinear models are applicable. Ursell [Urs53] proved that near the wave front U behaves like U ∼ t 3 .
Hence, regardless of how small nonlinear eﬀects are initially, they will become important. Recently, the methodology used in this thesis for tsunami generation was applied in the framework of the Boussinesq equations [Mit07] over uneven bottom. These equations have the advantage of being two-dimensional while the physical problem and, consequently, the potential flow formulation are 3D. On the other hand, Boussinesq equations represent a long wave approximation to the complete water-wave problem.

Source model

The inversion of seismic wave data allows the reconstruction of permanent deformations of the sea bottom following earthquakes. In spite of the complexity of the seismic source and of the internal structure of the earth, scientists have been relatively successful in using simple models for the source. One of these models is Okada’s model [Oka85]. Its description follows.
The fracture zones, along which the foci of earthquakes are to be found, have been described in various papers. For example, it has been suggested that Volterra’s theory of dislocations might be the proper tool for a quantitative description of these fracture zones [Ste58]. This suggestion was made for the following reason. If the mechanism involved in earthquakes and the fracture zones is indeed one of fracture, discontinuities in the displacement components across the fractured surface will exist. As dislocation theory may be described as that part of the theory of elasticity dealing with surfaces across which the displacement field is discontinuous, the suggestion makes sense.
As is often done in mathematical physics, it is necessary for simplicity’s sake to make some assumptions. Here we neglect the curvature of the earth, its gravity, temperature, magnetism, non-homogeneity, and consider a semi-infinite medium, which is homogeneous and isotropic. We further assume that the laws of classical linear elasticity theory hold.
Several studies showed that the eﬀect of earth curvature is negligible for shallow events at distances of less than 20◦ [BMSS69, BMSS70, SM71]. The sensitivity to earth topogra-phy, homogeneity, isotropy and half-space assumptions was studied and discussed recently [Mas03]. A commercially available code, ABACUS, which is based on a finite element model (FEM), was used. Six FEMs were constructed to test the sensitivity of deformation predictions to each assumption. The author came to the conclusion that the vertical layer-ing of lateral inhomogeneity can sometimes cause considerable eﬀects on the deformation fields.
The usual boundary conditions for dealing with earth problems require that the surface of the elastic medium (the earth) shall be free from forces. The resulting mixed boundary-value problem was solved a century ago [Vol07]. Later, Steketee proposed an alternative method to solve this problem using Green’s functions [Ste58].

READ Near surface geophysical imaging and characterisation

Volterra’s theory of dislocations

In order to introduce the concept of dislocation and for simplicity’s sake, this section is devoted to the case of an entire elastic space, as was done in the original paper by Volterra [Vol07].
Let O be the origin of a Cartesian coordinate system in an infinite elastic medium, xi the Cartesian coordinates (i = 1, 2, 3), and ei a unit vector in the positive xi−direction. A force F = F ek at O generates a displacement field uki (P, O) at point P , which is determined by the well-known Somigliana tensor

Dislocations in elastic half-space

When the case of an elastic half-space is considered, equation (1.5) remains valid, but we have to replace σijk in Tik by another tensor ωijk . This can be explained by the fact that the elementary solutions for a half-space are diﬀerent from Somigliana solution (1.1).
The ωijk can be obtained from the displacements corresponding to nuclei of strain in a half-space through relation (1.2). Steketee showed a method of obtaining the six ωijk fields by using a Green’s function and derived ω12k, which is relevant to a vertical strike-slip fault (see below). Maruyama derived the remaining five functions [Mar64].
It is interesting to mention here that historically these solutions were first derived in a straightforward manner by Mindlin [Min36, MC50], who gave explicit expressions of the displacement and stress fields for half-space nuclei of strain consisting of single forces with and without moment. It is only necessary to write the single force results since the other forms can be obtained by taking appropriate derivatives. The method consists in finding the displacement field in Westergaard’s form of the Galerkin vector [Wes35]. This vector is then determined by taking a linear combination of some biharmonic elementary solutions. The coeﬃcients are chosen to satisfy boundary and equilibrium conditions. These solutions were also derived by Press in a slightly diﬀerent manner [Pre65].

Finite rectangular source

Let us now consider a more practical problem. We define the elementary dislocations U1, U2 and U3, corresponding to the strike-slip, dip-slip and tensile components of an arbitrary dislocation. In Figure 1.1 each vector represents the direction of the elementary faults. The vector D is the so-called Burger’s vector, which shows how both sides of the fault are spread out: D = u+ − u−.
A general dislocation can be determined by three angles: the dip angle δ of the fault (0 ≤ δ ≤ π), the slip or rake angle θ (0 ≤ θ ≤ π), and the angle φ between the fault plane and Burger’s vector D. When dealing with a geophysical application, an additional angle, the azimuth or strike (see Figure 1.2 for strike angle definition), is introduced in order to provide an orientation of the fault. The general situation is schematically described in Figure 1.3.
For a finite rectangular fault with length L and width W occurring at depth d (Fig-ure 1.3), the deformation field can be evaluated analytically by a change of variables and by integrating over the rectangle. This was done by several authors [Oka85, Oka92, Chi63, SM74, IS79]. Here we give the results of their computations. The final results are repre-sented below in compact form, using Chinnery’s notation to represent the substitution f (ξ, η) = f (x, p) − f (x, p − W ) − f (x − L, p) + f (x − L, p − W ),

Curvilinear fault

In the previous subsection analytical formulas for the free-surface deformation in the special case of a rectangular fault were given. In fact, Volterra’s formula (1.5) allows to evaluate the displacement field that accompanies fault events with much more general geometry. The shape of the fault and Burger’s vector are suggested by seismologists and after numerical integration one can obtain the deformation of the seafloor for more general types of events as well.
Here we will consider the case of a fault whose geometry is described by an elliptical arc (see Figure 1.7).
The numerical integration was performed using a 9-point two-dimensional Gauss-type integration formula. The result is presented on Figure 1.8. The parameter values are given in Table 1.2.
The example considered in this subsection may not be physically relevant. However it shows how Okada’s solution can be extended. For a more precise modeling of the faulting event we need to have more information about the earthquake source and its related parameters.
After having reviewed the description of the source, we now switch to the deformation of the ocean surface following a submarine earthquake. The traditional approach for hy-drodynamic modelers is to use elastic models similar to the model we just described with the seismic parameters as input in order to evaluate the details of the seafloor deformation. Then this deformation is translated to the free surface of the ocean and serves as initial condition of the evolution problem described in the next section.

Table of contents :

Introduction
Propagation of tsunamis
Classical formulation
Dimensionless equations
Shallow-water equations
Boussinesq equations
Korteweg–de Vries equation
Energy of a tsunami
Tsunami run-up
1 Tsunami generation
1.1 Waves generated by a moving bottom
1.1.1 Source model
1.1.2 Volterra’s theory of dislocations
1.1.3 Dislocations in elastic half-space
1.1.3.1 Finite rectangular source
1.1.3.2 Curvilinear fault
1.1.4 Solution in fluid domain
1.1.5 Free-surface elevation
1.1.6 Velocity field
1.1.6.1 Pressure on the bottom
1.1.7 Asymptotic analysis of integral solutions
1.1.8 Numerical results
1.2 Comparison of tsunami generation models
1.2.1 Physical problem description
1.2.2 Linear theory
1.2.3 Active generation
1.2.4 Passive generation
1.2.4.1 Numerical method for the linear problem
1.2.5 Nonlinear shallow water equations
1.2.6 Mathematical model
1.2.7 Numerical method
1.2.8 Numerical method for the full equations
1.2.9 Comparisons and discussion
1.2.10 Conclusions
1.3 Tsunami generation by dynamic displacement of sea bed
1.3.1 Introduction
1.3.2 Mathematical models
1.3.2.1 Dynamic fault model
1.3.2.2 Fluid layer model
1.3.3 Numerical methods
1.3.3.1 Discretization of the viscoelastodynamic equations
1.3.3.2 Finite-volume scheme
1.3.4 Validation of the numerical method
1.3.5 Results of the simulation
1.3.6 Conclusions
2 Dissipative Boussinesq equations
2.1 Introduction
2.2 Derivation of the Boussinesq equations
2.2.1 Asymptotic expansion
2.3 Analysis of the linear dispersion relations
2.3.1 Linearized potential flow equations
2.3.2 Dissipative Boussinesq equations
2.3.3 Discussion
2.4 Alternative version of the Boussinesq equations
2.4.1 Derivation of the equations
2.5 Improvement of the linear dispersion relations
2.6 Regularization of Boussinesq equations
2.7 Bottom friction
2.8 Spectral Fourier method
2.8.1 Validation of the numerical method
2.9 Numerical results
2.9.1 Construction of the initial condition
2.9.2 Comparison between the dissipative models
2.10 Conclusions
3 Two phase flows
3.1 Introduction
3.2 Mathematical model
3.2.1 Sound speed in the mixture
3.2.2 Equation of state
3.3 Formal limit in barotropic case
3.4 Finite volume scheme on unstructured meshes
3.4.1 Sign matrix computation
3.4.2 Second order scheme
3.4.2.1 Historical remark
3.4.3 TVD and MUSCL schemes
3.4.3.1 Green-Gauss gradient reconstruction
3.4.3.2 Least-squares gradient reconstruction method
3.4.3.3 Slope limiter
3.4.4 Diffusive fluxes computation
3.4.5 Solution interpolation to mesh nodes
3.4.6 Time stepping methods
3.4.7 Boundary conditions treatment
3.5 Numerical results
3.5.1 Convergence test
3.5.2 Falling water column
3.5.3 Water drop test case
3.6 Conclusions
4 Viscous potential flows
4.1 Introduction
4.2 Anatomy of dissipation
4.3 Derivation
4.3.1 Dissipative KdV equation
4.4 Dispersion relation
4.4.1 Discussion
4.5 Attenuation of linear progressive waves
4.6 Numerical results
4.6.1 Approximate solitary wave solution
4.6.2 Discussion
4.7 Conclusion
Direction for future research
A VFFC scheme
A.1 Discretization in the finite volume framework
A.1.1 The one dimensional case
A.1.2 Extension to the multidimensional case
A.1.3 On the discretization of source terms
A.2 On the discretization of boundary conditions
Bibliography