# Analytical Solutions of Seismic-Generated Internal Tsunamis

Authors: Chi-Min Liu
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Internal tsunamis induced by an instantaneous seabed deformation are theoretically analyzed in this paper. Based on Hammack’s mathematical methods and findings, an in-depth and extensive study is performed to examine tsunami waveforms, especially at the initial stage. Waveforms of surface and internal tsunamis induced by three fundamental seabed uplifts, the rectangle-, cosine-, and sine-shape deformations, are solved to constitute an important base for analyzing waves generated by arbitrary seabed movements. A closed-form solution for the rectangle-shape deformation and analytic solutions for cosine and sine cases are obtained. The effects of spatial parameters on waveforms and the contributions of two frequency modes are investigated for the Boussinesq limit. The derived wave solutions not only improve the understanding of the formation of internal tsunamis but also provide an exact initial waveform for simulating wave propagation by various wave models.


---

## Body

## 1. Introduction

Tsunamis are one of the most catastrophic natural disasters for human beings. In this century, several destructive tsunamis, which include the 2004 Indian Ocean, 2010 Chile, and 2011 Japan Tohoku tsunamis, cause heavy casualties and destruction of infrastructures and buildings. In general, the development of tsunamis is distinguished into three phases [1]: generation by seabed deformation and propagation near the source, free propagation in the open ocean at large depths, and propagation and run-up in coastal regions. For tsunamis over a homogeneous single fluid, extensive studies have been implemented for tsunami formation [2–12], free propagation [4, 10, 13–16], and run-up in coastal regions [17–21].Among three phases in tsunami life, understanding the behaviors of tsunamis at the generation stage is crucial for accurately analyzing subsequent wave propagation as the error made in this phase cannot be corrected by numerical methods and will propagate in space and time [22]. In general, underwater earthquakes and submarine landslides are two major sources causing topographic variations to trigger tsunamis. The present paper will focus on the tsunami formation by underwater earthquakes by which 75 percent of historical tsunamis were generated [23]. Emphases are placed on the examination of wave features at the initial stage. In some past studies, the initial tsunami profile is assumed to be identical to the seafloor deformation due to the reasons that the sea depth is sufficiently shallow in comparison with tsunami wavelength, and the seabed deformation occurs instantaneously or within a very short period [9, 24].For a single-fluid system, Hammack [3] studied tsunami generation by theoretical modeling and laboratory measurements. The linear potential theory was adopted to examine tsunamis induced by both of the exponential- and cosine-type rectangle-shape seabed uplifts. He noted that the validity of linear theory requires that the total amplitude of the bed displacement must remain small compared with the uniform water depth. A time-size ratio is proposed to distinguish the movements into fast and creeping motions. In Hammack’s derivation, the term 1/coshkh (k is the wavenumber and h the water depth) appearing in wave solutions demonstrates that the difference between the initial waveform and seafloor deformation will increase as water depth becomes larger even the seafloor movement is instantaneous. It means that the conventional idea of directly assigning the seafloor displacement to be the initial waveform must be carefully considered and should be modified by the factor 1/coshkh especially for tsunamis generated in deeper water. Dutykh et al. [22] numerically examined the effect of 1/coshkh by simulating wave propagation by various numerical models. Their numerical results show that the difference between the initial waveform with the inclusion of the effects 1/coshkh and the seabed deformation is of the order of 20% for a specific scenario. Kervella et al. [5] applied three wave models, which include linear shallow-water equations, nondispersive nonlinear shallow-water equations, and Boussinesq equations, to simulate tsunami propagations for examining nonlinear effects and the influence of 1/coshkh. They found that, for most tsunami events, the linear theory is sufficient, and the dynamic effects of seabed motion should be taken into account.In addition to the rectangular-type seabed deformations, more topographic variations were also considered to study tsunami features. Todorovska and Trifunac [12] theoretically investigated the lateral spreading velocity of the deformed seabed for fast and slow submarine faulting. Le Gal et al. [6] generalized the generation mechanisms proposed by [3, 12] to study the effects of the rise time in the vertical direction and the rupture propagation velocity in the horizontal direction. Besides, the Gaussian distribution as a seabed topography with a constant lifting speed was considered by [9] to investigate the criteria for small source size and short source duration.Recently, Derakhti et al. [2] examined the tsunami generations by three seafloor deformations of different shapes and simulated the wave propagation using the fully nonlinear smooth particle hydrodynamic model. The effects of spatial and temporal parameters on waveforms in near fields are investigated. In addition to wave evolutions, velocity fields are also extensively examined. Liu [7] theoretically analyzed the waveforms of tsunamis induced by three fundamental seabed deformations including rectangle-, cosine-, and sine-shape seabed uplifts. For instantaneous seabed deformations, the waveforms at the initial stage are obtained and shown in a form without integrations. The analytical solutions not only constitute an important base for analyzing the tsunami generation by arbitrary topographic variations with the help of Fourier analysis but also provide a fast way and, more importantly, an exact initial condition for performing subsequent simulation by various wave models. Tsunamis induced by noninstantaneous seabed movements were also investigated in his study.In contrast with a large number of papers focusing on tsunamis’ generation and propagation at the free surface, internal tsunamis, which are generated by seafloor movements and then propagate between two fluids, attract comparatively less attention in past studies. Hammack [25] examined barotropic (surface) and baroclinic (internal) tsunamis induced by an instantaneous rectangle-shape seabed uplift in a two-fluid system. An exponential form of the seabed uplift is adopted for the generation. By applying the integral-transform method, the solutions for surface and interfacial waves were derived. In the Boussinesq limit which denotes the density ratio of upper and lower fluids approaches unity, the initial wave elevations at the deformation edge were presented. The barotropic-mode wave motion is found to be not affected by the small differences in potential density typical of ocean stratification, while the baroclinic wave amplitude is proportional to the vertical offset of the seafloor. Song [26] also applied the integral transforms to solve the same problem by examining three scenarios of different seafloor movements. Based on the derived integral solutions, near-field waveforms are plotted for specific spatial, temporal, and physical parameters. Though the above papers reveal some internal tsunami features in a two-fluid system, more detailed examinations are still required to improve our understanding of internal tsunamis. In addition to theoretical analyses, very few papers were presented for the observation of real internal tsunamis [27–29].This study aims to analyze internal and surface tsunamis generated by instantaneous seabed deformations in a two-fluid system. The organization of the present paper is as follows: in Section2 the linear potential theory is adopted to derive the solutions of tsunamis with the help of the Laplace–Fourier transform. Three different types of seabed deformations are considered and analyzed. The derived wave solutions contain the components of fast and slow modes. Section 3 examines the wave features in the Boussinesq limit. Waveforms at the initial stage are solved to be expressed in concise forms (for rectangle-shape uplift, a closed-form solution is obtained). The effects of spatial parameters on tsunami waveforms are analyzed. Components of two frequency modes are examined separately. Conclusions are made in Section 4.


## 2. Mathematical Modeling

Consider a one-dimensional horizontal (1DH) fluid system possessing two layers of fluids which are assumed to be irrotational and inviscid. The notations ofϕi, ρi, and hi represent the velocity potential, density, and undisturbed layer thickness in which the subscript i=1 and 2 denote the upper and lower fluids, respectively. The compressibility of fluid is neglected as its major effect is on the excitation of ocean acoustic waves [10]. The density ratio ρr=ρ1/ρ2 is assumed to be less than unity. The total depth is denoted by h=h1+h2. The free-surface and interfacial wave elevations are represented by η1 and η2, respectively. The horizontal coordinate is denoted by x and the vertical coordinate z starts at the undisturbed interface and points upward. The linear potential theory is adopted for analysis due to the assumptions of a very long wavelength of a tsunami and a very short duration of seafloor deformation. Based on the linear wave theory, the governing equations, boundary, and initial conditions are(1)∂2ϕ1∂x2+∂2ϕ1∂z2=0,∂2ϕ2∂x2+∂2ϕ2∂z2=0,∂ϕ1∂z=∂η1∂t,at,z=h1,∂ϕ1∂t+gη1=0,at,z=h1,∂ϕ1∂z=∂η2∂t,at,z=0,∂ϕ2∂z=∂η2∂t,at,z=0,ρ1∂ϕ1∂t+gη2=ρ2∂ϕ2∂t+gη2,at,z=0,∂ϕ2∂z=∂ζ∂t,at,z=−h2,η1=η2=ζ=0,at,t<0,where ζ is the seabed deformation and g the acceleration of gravity. By applying the Fourier–Laplace transform, defined as(2)f˜k,s=12π∫−∞∞eikxdx∫0∞fx,te−stdt,to the above equations and eliminating η1 and η2 from boundary conditions, we have(3)∂2ϕ˜1∂z2−k2ϕ˜1=0,(4)∂2ϕ˜2∂z2−k2ϕ˜2=0,(5)s2ϕ˜1+g∂ϕ˜1∂z=0,at,z=h1,(6)ρ1s2ϕ˜1+g∂ϕ˜1∂z=ρ2s2ϕ˜2+g∂ϕ˜2∂z,at,z=0,(7)∂ϕ2∂z=sζ˜,at,z=−h2.The elevations of surface and interfacial waves are determined by velocity potentials:(8)η˜1=−sgϕ˜1,at,z=h1,(9)η˜2=sρrϕ˜1−ϕ˜2g1−ρr,at,z=0.One can solveϕ˜1 and ϕ˜2 from equations (3) and (4) with the help of boundary conditions, equations (5)–(7) (details are referred to [25]). Substituting the solutions of ϕ˜1 and ϕ˜2 into equations (8) and (9) results in(10)η˜1k,s=ζ˜k,s⋅s4Θ⋅coshkh1⋅coshkh2,(11)η˜2k,s=ζ˜k,s⋅s4+s2gkT1Θ⋅coshkh2,where(12)Θ=s41+ρrT1T2+s2gkT1+T2+1−ρrg2k2T1T2,(13)T1=tanhkh1,T2=tanhkh2.Now the rectangle-shape seabed deformation is considered(14)ζx,t=ζ0⋅HB2−xHB2+xHt,where H. denotes the Heaviside step function, B the width of deformation, and ζ0 the constant lifting height. By substituting the transformed seabed deformation of equation (14)(15)ζ˜k,s=ζ02π⋅sinkB/2ks,into equations (10) and (11) and then taking the inverse Laplace transform, the surface and interfacial waveforms are obtained(16)η1x,t=ζ0Bπ∫0∞D1mcosωmt+D1ncosωntcoskxdk,(17)η2x,t=ζ0Bπ∫0∞D2mcosωmt+D2ncosωntcoskxdk,where frequencies ωm and ωn are(18)ωmωn=gkT1+T2±T1+T22−41−ρr1+ρrT1T2T1T221+ρrT1T2,which satisfy the dispersion relation(19)ω41+ρrT1T2−ω2gkT1+T2+1−ρrg2k2T1T2=0.Equations (16) and (17) reveal that the total waveforms are the combination of waves of all possible wavenumbers. The amplitude densities D1m, D1n, D2m, and D2n in the integrands of equations (16) and (17) are(20a)D1m=ωm2F1ωm2−ωn2,(20b)D1n=−ωn2F1ωm2−ωn2,(20c)D2m=ωm2F2−F3ωm2−ωn2,(20d)D2n=−ωn2F2+F3ωm2−ωn2,and(21a)F1=2sinkB/2kBcoshkh1coshkh21+ρrT1T2,(21b)F2=2sinkB/2kBcoshkh21+ρrT1T2,(21c)F3=2gT1sinkB/2Bcoshkh21+ρrT1T2.Let us come back to examine the wave solutions of equations (16) and (17). Each solution contains two components of different frequencies ωm and ωn implying two possible modes of wave motions. The former is commonly called the fast-mode (surface-mode) wave motion, and the latter indicates the slow-mode (internal-mode) wave motion [30]. For general ocean circumstance, ωm is much greater than ωn due to a slight difference of densities of both layers implying that the component of the fast mode will dominate the wave features. Detailed discussions are left to the next section.In addition to the rectangle-shape seabed uplift, another two fundamental seabed deformations are also considered:(22)ζcx,t=ζ0⋅cosnπxBHB2−xHB2+xHt,(23)ζsx,t=ζ0⋅sinnπxBHB2−xHB2+xHt,in which the subscripts c and s denote the cosine- and sine-shape seabed deformations, respectively. It is remarked that equations (14), (22), and (23) constitute a set of fundamental topographical variations which can be applied to analyze arbitrary seabed deformations with the help of Fourier analysis. The derivation and analysis for tsunamis induced by equations (22) and (23) are referred to Appendix for details.


## 3. The Boussinesq Limit

In this section, the analyses are performed to examine tsunamis’ features in real oceans where the Boussinesq limit (ρr⟶1) is approximately satisfied. In the Boussinesq limit, we focus on deriving the closed-form solutions of initial tsunami waveforms. The initial waveform means the wave elevation at the time when the seabed deformation is just completed (t=0+). Based on the closed-form solutions, corresponding wave features are investigated. As the accuracy of the simulation of wave propagation is significantly affected by initial waveforms, closed-form initial waveforms are expected to provide a fast and correct calculation of initial conditions for subsequent numerical simulation.The rectangle-shape seabed uplift described by equation (14) is adopted for analysis. First, a dimensionless parameter is introduced(24)∈=1−ρr<1.With the consideration of equation (24), the amplitude densities of equations (20a)–(20d) can be approximated by expanding them in terms of ∈0 and ∈1:(25a)D1m=2sinkB/2kBcoshkh1coshkh211+T1T2+∈T1T2T1+T22+1+T1T22T1+T221+T1T22+O∈2,(25b)D1n=−∈⋅2sinkB/2kBcoshkh1coshkh2T1T2T1+T22+O∈2,(25c)D2m=2sinkB/2kBcoshkh2T21−T12T1+T21+T1T2+∈⋅T1T21−T12T2−T1+T23+5T1T22+2T12T23T1+T231+T1T22+O∈2,(25d)D2n=2sinkB/2kBcoshkh2T1T1+T2+∈⋅T1T21−T12T1−T2+2T12T2T1+T23+O∈2.The leading order ofD1m,D2m, and D2n is O∈0 while D1n is of O∈, implying that, at least, the fast mode dominates surface wave behaviors. Besides, by the O∈0 analysis, the ratio of D1m to the summation of D2m and D2n is(26)D1mD2m+D2n=1coshkh1.Equation (26) indicates the initial amplitude of surface tsunami is less than that of the internal tsunami. In the following propagation, the waveforms contributed by D2m and D2n will gradually separate due to the difference of phase velocities of two frequency modes.To quantitatively analyze the contributions of two frequency modes in wave solutions, the amplitude densities accurate up toO∈0 for deeper and shallower waters are plotted in Figures 1 and 2, respectively. It is noted again that the amplitude density reflects the distribution of wave amplitudes of all wave numbers. In Figure 1, spatial parameter ratios, h/B=0.04 and h1/h=0.025, are adopted. These spatial parameters approximately meet the general in situ conditions (B=100km, h=4km, and h1=100m are typical scales for general underwater earthquakes). Solid and dash lines represent the surface and interfacial wave behaviors, respectively. The profiles of D1m and D2m are very close to each other indicating the leading wave behaviors might be very similar. It is also found that the slow-mode contributions for both surface and interfacial waves (D1n and D2n) are much smaller than those of fast mode ones (D1m and D2m). This elucidates that the fast mode dominates wave features of both surface and internal tsunamis.Figure 1
Amplitude densities forh/B=0.04 and h1/h=0.025.Figure 2
Amplitude densities forh/B=0.01 and h1/h=0.1.As for shallower waters,h/B=0.01 and h1/h=0.1 (corresponding to B=100km, h=1km, and h1=100m), results are plotted in Figure 2. Most behaviors of amplitude densities are quite similar to those shown in the deeper-water case (Figure 1). The major difference is the growth of D2n in the range of small kB. It means that the amplitude of leading interfacial waves of the fast mode will become smaller in shallower waters.In addition to amplitude densities, initial waveforms of surface and internal tsunamis are subsequently analyzed. As the initial waveform is calculated by integrating the amplitude density associated withcoskx at t=0+, as shown in equations (16) and (17), similar behaviors of initial waveforms and amplitude densities are reasonably expected. To this end, the following analysis majorly focuses on the derivation of closed-form solutions of initial waveforms because the closed-form solutions are of great help for providing exact initial conditions for subsequent numerical simulation.In the Boussinesq limit (∈⟶0), the initial waveform for surface waves can be simplified from equations (16), (25a), and (25b):(27)η1x,0+=2ζ0π∫0∞sinkB/2coskxkcoshkh1+h2dk,which can be further integrated to be a closed form [7]:(28)η1x,0+=2ζ0π⋅tan−1sinhπB/4hcoshπx/2h.Note that the depth ratio does not affect the initial waveform of the surface tsunami. It is also remarked that equation (28) is identical to equation (5) in [7] for a single-layer fluid. As for interfacial waves, the initial waveform is rewritten from equation (17):(29)η2x,0+=2ζ0π∫0∞sinkB/2coskxkcoshkh21+T1T2dk,which can be simplified to be(30)η2x,0+=2ζ0π∫0∞sinkB/2coskxk⋅coshkh1coshkhdk.Equation (30) can be further integrated into a closed form by applying the integration formula 4.122.1 in [31](31)η2x,0+=ζ0π⋅tan−1sinhπB+2πx/4hcosπh1/2h+tan−1sinhπB−2πx/4hcosπh1/2h.It is remarked that Hammack [25] presented the elevations of surface and interfacial waves only at the position located at the edge of deformation (x=B/2) in his equations (29) and (30), which can be obtained from the present equations (28) and (31) as special cases. From equation (31), it is shown that, as the depth ratio h1/h increases, the initial amplitude of internal tsunami, η20,0+, also increases. Now the magnitudes of η10,0+ and η20,0+ are compared:(32)η10,0+η20,0+=tan−1sinhπB/4htan−1sinhπB/4h/cosπh1/2h.It shows that the initial amplitude ratio of the surface wave to interfacial wave decreases as the depth ratioh1/h increases. Though equation (31) is always smaller than unity, if one considers spatial conditions in most in situ conditions of B≫h and the behavior of arctangent function with a very large argument, this amplitude ratio is almost identical to unity even though the depth ratio h1/h is very large. Similar behaviors have also been found previously when examining the amplitude densities.Since internal tsunamis are constituted by two frequency modes of motions, the initial waveform contributed by the slow mode with the amplitude densityD2n shown in equation (25b) is(33)η2nx,0+=2ζ0π∫0∞sinkB/2coskxk⋅T1T1+T2dk,which can be rewritten as(34)η2nx,0+=2ζ0π∫0∞sinkB/2coskxk⋅sinhkh1sinhkhdk.By using the identitysinhkh1=i−1sinikh1 and performing the product-to-sum formula of trigonometric functions, equation (34) can be integrated with the help of the integration formula 4.121.2 in [31]. It results in(35)η2nx,0+=ζ02πilncoshπB+2x+2ih1/4h⋅coshπB−2x+2ih1/4hcoshπB+2x−2ih1/4h⋅coshπB−2x−2ih1/4h,which can be further simplified to be a real closed-form solution:(36)η2nx,0+=ζ02πcos−1coshπB/2hcosπh1/h+coshπx/h2−sinhπB/2hsinπh1/h2coshπB/2hcosπh1/h+coshπx/h2+sinhπB/2hsinπh1/h2.For the limiting case of eitherh1/h⟶0 or h1/h⟶1, equation (35) will approach zero indicating the contribution of the slow mode to internal tsunami can be ignored. The contribution from the fast-mode amplitude density G2m can be obtained by subtracting equation (31) from equation (36), which is neglected herein.Behaviors of initial waveformsη1x,0+ and η2x,0+ are examined in Figure 3. The ratios of η2 to ζ0 and η1 to η2 are plotted in the left and right panels, respectively. Note that η1 and η2 are symmetric with respect to x=0 for the rectangle-shape uplift. Solid and dash lines represent the results of h1/h=0.1 (a thinner upper layer) and 0.9 (a thicker upper layer), respectively, while the blue- and red-colored lines denote the results of h/B=0.01 (shallower water) and 0.1 (deeper water), respectively. For internal tsunamis shown in the left panel, initial waveforms are closer to the seabed deformation (a rectangle-shape uplift) for either larger h1/h or smaller h/B. In other words, the dynamic effects of the seabed movements are stronger either for the case of either deeper water or thinner-upper layer. In the right panel, the ratio of η1 to η2 for the case of h1/h=0.1 almost remains constant in and out of the deformation range (x≤B/2). As for the case h1/h=0.9, the ratios are also close to unity near the deformation center, but the gap between the two waveforms becomes obvious around the edge of the deformation. Though the ratios for this case begin to quickly grow when x>B/2, it is noted that the elevations of both tsunamis are very small out of the deformation zone implying that the difference between elevations of two tsunamis is exaggerated in the ratio presentation.Figure 3
Waveforms of surface tsunamiη1x,0+ and interfacial tsunami η2x,0+.As there exist two frequency modes of motion for internal tsunamis, the initial waveform of internal tsunami is separated into the fast-mode and slow-mode components for analysis herein. In Figure4, the right panel shows the results of the slow-mode initial waveform η2nx,0+ and the left panel displays the results of the fast mode initial waveform η2mx,0+. The former is directly calculated by equation (36) while the latter is obtained by subtracting equation (36) from equation (31). Four cases with different ratios of spatial parameters are identical to those adopted in Figure 3. Above the deformation zone, the magnitudes of the fast mode are found to be much greater than those of the slow mode for all cases. For the slow-mode amplitude, it is found that the waveforms for h1/h=0.1 are comparatively smoother than those for h1/h=0.9 as the waveforms of the latter ones abruptly oscillate near the edge of the seabed deformation.Figure 4
Interfacial waveforms of two wave modes.


## 4. Conclusions

Internal and surface tsunamis generated by an instantaneous seabed deformation are theoretically investigated in this paper. A fundamental rectangle-shape uplift is adopted as the seafloor variation while cosine- and sine-shape uplifts are also analyzed in Appendix. The linear wave theory is applied to derive the analytic wave solutions for all times. Major emphases are placed on analyzing the initial waveforms and features. In the Boussinesq limit, the closed-form solution for the rectangle-shape uplift and analytic solutions for cosine- and sine-shape deformations are obtained. According to the derived solutions, major wave features and future works are summarized as follows:(1)
Effects of spatial parameters: For surface tsunamis, the depth ratioh1/h does not affect the initial waveform. The difference between the initial waveform and seabed deformation increases as h/B grows. As for internal tsunamis, the wave elevation increases when the depth ratio h1/h grows. The difference between the initial waveform and seabed deformation becomes large as either h1/h or h/B increases. Moreover, the initial amplitude of surface tsunami is less than but very close to that of the internal tsunami.(2)
Effects of two frequency modes: In a two-fluid system, there exist two frequency modes of wave motions. In the Boussinesq limit, wave features of surface tsunamis are almost dominated by the fast mode. As for internal tsunamis, both frequency modes contribute to waveforms while the fast mode has larger contributions than the slow mode. Ash1/h⟶0 or 1, the initial amplitude will be totally contributed by the fast-mode component.(3)
Tsunami propagation: In fact, the analysis of tsunami propagation is not substantially performed in the present paper. This study aims to provide a set of exact solutions at the initial stage as exact initial conditions for manipulating the following simulations. Based on the derived initial wave solutions, one may perform the theoretical analysis to examine wave propagation features. For leading waves in a single-fluid system, an approximate solution is theoretically derived and shown in an integral form of the Airy function for the rectangle-shape uplifts [14, 32]. A similar analysis might be implemented for the case of a two-fluid system. On the other hand, one can perform numerical simulations to calculate wave propagation by using present solutions as initial conditions.(4)
Future works: Though most of the physical characteristics of tsunamis, especially the surface tsunamis, have been extensively studied in past decades, efforts are still required to investigate more wave features of internal tsunamis. For example, how does the Coriolis force affect the initial waveforms and the propagation features? What are the tsunami behaviors in a density-stratified fluid? These topics will be studied in the near future.



---

*Source: 1113733-2021-10-06.xml*