CC BY
15
RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 8, ES5003, doi:10.2205/2006ES000210, 2006
Numerical simulation of the 7 February 1963 tsunami in the Bay of Corinth, Greece
L. I. Lobkovsky1, R. Kh. Mazova2, I. A. Garagash3, and L. Yu. Kataeva2 Received 21 October 2006; accepted 30 October 2006; published 28 November 2006.
[1] The paper addresses the event that occurred in the Bay of Corinth, Greece, on 7 February 1963. The tsunami was produced by a submarine landslide in the mouth of the local Salmenikos River. The paper presents preliminary numerical results obtained in terms of an elastoplastic model for the surface water waves generated by the landslide and propagated onto the coast. The wave propagation is analyzed as a function of the landslide movement. The inferred results are compared with those obtained in terms of other models.
INDEX TERMS: 3070 Marine Geology and Geophysics: Submarine landslides; 3285 Mathematical Geophysics: Wave propagation; 4255 Oceanography: General: Numerical modeling; 4564 Oceanography: Physical: Tsunamis and storm surges; KEYWORDS: tsunami generation, shallow water equations, aseismic tsunami, sedimentary mass.
Citation: Lobkovsky, L. I., R. Kh. Mazova, I. A. Garagash, and L. Yu. Kataeva (2006), Numerical simulation of the 7 February 1963 tsunami in the Bay of Corinth, Greece, Russ. J. Earth. Sci., 8, ES5003, doi:10.2205/2006ES000210.
Introduction
[2] Recently the problem of tsunami generation due to seismic or aseismic submarine rockfalls or coastal landslides has received much attention. An example is the catastrophic tsunami produced by the 17 July, 1998 (M = 7.1) earthquake in Papua New Guinea (PNG) (e.g., see [Mazova et al., 2004]). Although the PNG tsunami generation mechanism is still debated, several observation and model data indicate that a submarine landslide was the tsunami cause [Mazova, 2003; Papadopoulos, 2000]. Moreover, detailed descriptions of other cases of landslide-induced tsunamis have been published over the last 50 years (for example, in the South Aegean Sea (Greece) in 1956; in Lituya Bay (Alaska) in 1958, in the city of Aegion (Greece) in 1963; in the Vayont valley (Italy) in 1963; in Norway and Nice (France) in 1979; in Skagway Harbor (Alaska) in 1994; in Izmit (Turkey) in 1999; and Fatu Hiva (French Polynesia) in 1999 (e.g. see [Fine et al, 1998; Mazova, 2003; Murti, 1981; Papadopoulos, 2000])).
[3] In particular, as was shown in [Papadopoulos, 2000], a landslide-induced tsunami is a rather frequent phenomenon
1P. P. Shirshov Institute of Oceanology, Russian Academy of Sciences, Moscow, Russia
2Nizhny Novgorod State Technical University, Nizhniy Novgorod, Russia
3Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russia
Copyright 2006 by the Russian Journal of Earth Sciences.
ISSN: 1681-1208 (online)
off the Greece coasts. For example, the tsunami wave that arose after the strong (M = 7.5) 7 July, 1956 earthquake in the South Aegean Sea was due to not only a seismic fault motion but also a seismically induced slide of submarine sedimentary masses. Another example is the local strong tsunami generated by the aseismic landslide of 7 February 1963, in the Aegion area of the western Bay of Corinth (central Greece).
[4] Statistical analysis and computations show that submarine and subaerial landslides can generate tsunami waves of a considerable height at the coast near their source [Mazova, 2003]. The length of such a wave is a few kilometers, and its height rapidly decreases due to frequency dispersion. On the contrary, the propagation range of tsunamis generated by a fault is very long, and their heights in open ocean are relatively small. Their wavelength is of the order a few hundreds of kilometers, while the height of wave decreases slowly because of the weakness of dispersion effects. A landslide-generated tsunami often results in catastrophic consequences [Mazova, 2003] and, therefore, the problem of landslide-generated waves along an open coastline or in a closed water area is of great practical interest for engineering in coastal zones.
[5] Significant progress has been made in the last decade in both numerical simulation of landslide-induced tsunamis and development of analytical models for an adequate description of the tsunami generation [Mazova, 2003]. The most often used are the rigid block model [Iwasaki, 1997] and viscous or viscoplastic models [Jiang and LeBlond, 1993, 1994]. However, as was noted by Fine et al, [1998], the rigid block model overestimates the responses of the water surface to submarine disturbances, whereas viscous/viscoplastic
Figure 1. Schematic map of the Bay of Corinth.
models underestimate it. Therefore, to adequately model a tsunami from a submarine landslide, a model is required that takes into account both the detailed structure of the landslide body and the mechanical properties of landslide-body constituents characterizing its behavior during the landslid-ing process.
[6] In this work, we perform numerical simulation of the generation and propagation of tsunami waves due to sliding of the submarine portion of a continental slope, with the landslide behavior being described in terms of an elasto-plastic model. As a case study, we chose the tsunami of 7 February 1963, in the Bay of Corinth (central Greece). We considered the shoreline segment shown by an arrow in the schematic map of the Bay of Corinth (Figure 1) [Mazova et al., 2004]. The landslide motion and the generation and propagation of the tsunami were computed in a 1-D approximation.
bay (Figure 1). Intense tectonic activity often associated with fractures in sedimentary masses leads to the formation of steep submarine slopes. As a result, generation of both seismic and aseismic (landslide-induced) tsunamis is an additional geophysically important feature of the Bay of Corinth.
[8] A typical landslide-induced tsunami in the Bay of Corinth was the wave of 7 February 1963, whose characteristics have been studied in detail [Papadopoulos, 2000]. The maximum wave height reached 5-6 m at the southern coast. Field observations [Galanopoulos et al., 1964] leave no doubts that the tsunami was caused by a sudden aseis-mic seaward landslide of the submarine slope overloaded by sedimentary masses (57,000 m3 in volume) that were accumulated in the mouth of the local Salmenikos River on the southern coast. The sliding masses were mostly submarine.
The Tsunami of 7 February 1963 in the Bay of Corinth (Greece)
[7] The Bay of Corinth is a scoop-like tectonic structure rapidly widening in the southerly direction due to the stress field, and this results in high seismicity of the region related to an E-W trending fault. The fault is located between the mainland of Greece in the north and the Peloponessus Peninsula in the south and is enclosed in a water basin supplying sediments the northern and southern coasts of the
Numerical Simulation and Results
[9] The model used for the simulation of the wave of 1963 is based on the assumption of layered sediments resting on a relatively rigid base (Figure 2). The parameters involved in the calculations were the density of a layer ps, shear modulus G, bulk modulus K, cohesion c, maximum friction angle ipa, and tensile strength a [Garagash et al., 2003]. It is supposed that the volume of the landslide body does not depend on the disturbance magnitude and is rather controlled by mechanical properties of its components, for example, the angle
Figure 2. Scheme used for calculating the submarine slope. The upper layer is the plastic landslide mass resting on an elastic substrate.
of internal friction. An additional dynamic effect can be a triggering mechanism exerting a pushing action on the landslide mass on a slope until the critical stress is exceeded.
[10] After the situation becomes unstable, the sliding process continues due to the accumulated potential energy. We introduced an initial dynamic action on the landslide process equivalent to a moderate earthquake and lasting for 6 s (Figure 2b) [Garagash et al., 2003].
[11] The stability of submarine slopes was analyzed using the FLAC code [Garagash et al., 2003] implementing an explicit 3-D finite difference scheme applied in geomechanical modeling. In contrast to the finite element method, the explicit difference scheme enables the modeling of nonlinear behavior of pore-saturated sedimentary masses under the condition of plastic flow above the yield strength [Garagash and Ermakov, 2001; Garagash et al., 2003]. The modeling was performed in two stages. At the first stage, we modeled the initial prestressed state of the sedimentary slope formed due to the gravity force and the seawater saturation under the sea pressure conditions. We considered layer-by-layer sliding of an upper part of the elastoplastic sedimentary layer on the slope surface that forms during the landslide process; a distinct interface is supposed to exist between water and the landslide body. The water density pw(x,z) is constant and the landslide density is a function of coor-
dinates: ps = ps (x,z). The sliding process of sedimentary masses is essentially controlled by the friction angle. In this work, the numerical simulation was performed for the maximum friction angle = 20°. Assuming that the process stops whenever the velocity of the leading front of landslide masses becomes equal to zero, we found that the landsliding process lasts for 45 s. The maximum velocity of the leading front of the landslide is about 6 m s-1 and is attained at about 40 s after the beginning of the process. Figures 3-5 show the evolution of particle velocities in the landslide body during the sliding of the sedimentary layer, the plastic strain of the landslide body during the sliding process, and the distribution of the shear strain for the angle = 20°and 32°. (The computations were performed for four time moments (10, 25, 40, and 55 s; the sliding process virtually stopped at the latter time moment.)
[12] At the second stage, we performed a numerical simulation of landslide-generated surface water waves on the basis of the nonlinear system of shallow water equations with the use of a specially modified explicit difference scheme including the check for stability conditions (see below).
[13] In the computations, the origin of coordinates was chosen on the shoreline; the seaward x axis coincided with the level of the undisturbed water surface; and the z axis was directed upward (Figure 6), where z = -h(x,t) is the
Figure 3. Evolution of particles displacements within the landslide body during its movement.
variable water depth and z = —h(x, 0) = hs(x) is the profile of the submarine slope that existed before the disturbance onset (t = 0). D(x,t) is the thickness of the landslide body: h(x,t) = hs(x,t) — D(x,t). The shape of the free water surface is denoted as z = n(x,t) with n(x, 0) = 0. Figure 5 shows the successive positions of the landslide surface at a time step of 10 s. It is well seen that the landslide moves downward on the slope, continuously changes its shape and, shifting to the right, reaches the final position at a distance of about 1300 m from the shoreline (x = 0). It is noteworthy that the upper part of the sliding body seems to move leftward as well, involving the dry coast. As a result, the
immobile shoreline point (x = 0) becomes mobile, moving leftward and reaching the position x = —200 m at the final moment of the landsliding process (t = 45 s). The same effect is observed in a purely viscous model of a landslide. However, a landslide in the viscoplastic model moves only downward on the slope [Jiang and LeBlond, 1993, 1994].
[14] The above modeling results obtained for the landslide movement were used for numerical simulation of the landslide-generated surface water wave. The wave generation was computated using the ordinary nonlinear system of shallow water equations, which can be written as
du du dn
"oT + uj)—+ = 0
at ox ox
Figure 4. Plastic deformation of the moving landslide body.
where x is the horizontal coordinate and u is the horizontal component of the water particle velocity in the wave. The coupling between the landslide masses and seawater during the surface water wave generation was described through the related continuity equation:
dn dD dhs d .
dt = ~at — ~dt — dx ((hs + n — D)u
in time and the second order in spatial coordinates:
,,n+1
= u"— du" (u
d ( n n\,/n I n r» n\
~7T~ (ti+1 — ti-1) +c(u;+1 + ui_1 — 2ui),
Br
n+1
nni — d
(2)
(n"+1 + hni+1 — Dn+1)
+1 ui+1
We should note here that, in contrast to the existing models [Fine et al., 1998; Iwasaki, 1997; Jiang and LeBlond, 1993, 1994], the second term on the right-hand side of (2) does not vanish because the landslide in the model under consideration moves on a surface formed by the sliding process rather than on a stationary slope surface, as is generally supposed [Fine et al., 1998; Iwasaki, 1997; Jiang and LeBlond, 1993, 1994].
[15] To solve numerically the system of equations (1) and (2), we used an explicit two-layer scheme of the first order
— (tn-1 + hni-1 — Dn-1)un-1
+
_(hn+1 — Dn + c(n?+1 + n?-:
n) — () — 2n.n)
Din
where
i = 1, Nx — 1; n = 1, Nt — 1; d =
At
2Ah
A 1Ah2 At 2At Ah2
A1
T
(3)
n
n
u
1
n
c
(a)
/ = 10 s
(p = 20
t = 25 s
/ = 40 s
t =55 s
Contour of Shear Strain Increment Gradient Calculation
HO. OOOOe+OOO to 5.0000e-002 1 .0QQQe-0Q1 to 1.5Q00e-0Q1
■ 2.00000-001 to 2.50000-00 1
■ 3.0000e-001 to 3,5000e-001 _ 4.00006-001 to 4.5000e-001 ~\ 5 .OOOOe-OO1 to 5.5000 e-00 1
8.0000e-001 to S.5Q00e-0Q1 7.00000-001 to 7.50000-001 7.5000e-001 to 7.8S76e-001 nterval = 5.0e-002
Contour of Shear Strain Increment
Gradient Calculation
1 .7489e-006 to 5.0000e-001 5.00000-001 to 7.50000-001 7.5000e-001 to 1.0000e+000 1 .00000+000 to 1 .25000+000 1 .2500e+000 to 1.5000e+000 1 .50006+000 to 1.75006+000 1 .75000+000 to 2.00000+000 2.00006+000 to 2.06536+000 nterval = 2.5e-001
Contour of Shear Strain Increment
Gradient Calculation
■ -9.7787e-001 to 0.0000e+000
■ 0.OOOOe+OOO to 5.0000e-001 H 5.00006-001 to 1.00006+000 H 1 .OOOOe+OOO to 1 ,5000e+000 B 1.5000e+000 to 2.00QQe+000 H 2.00000+000 to 2.50000+000
■ 2.500Oe+000 to 3.0000e+000 3.OOOOe+OOO to 3.5000e+000 3.50000+000 to 4.00000+000 4.000Qe+000 to 4.5QQQe+Q0Q 4.50000+000 to 5.00000+000 5.OOOOe+OOO to 5.S000e+000 5.50006+000 to 6-OOOOe+OOO 6.00000+000 to 6.50000+000 6.5000e+000 to 6.8849e+QQQ nterval = 5.00-001
Contour of Shear Strain Increment
Gradient Calculation
■ 3.5071e-006 to 1 .OOOOe+OOO
■ 1 .0000e+000 to 1 .50006+000 H 1 ,5000e+000 to 2.00000+000 I 2.QQQ0e+00Q to 2.500Qe+QQQ I 2.5000e+000 to 3.00000+000 I 3.OOOOe+OOO to 3 .S000e+000 H 3.5000e+000 to 4.00006+000 _ 4.OOOOe+OOO to 4.5000e+000 “ 4.50006+000 to 5.00006+000 _ 5.OOOOe+OOO to 5.50000+000 _ 5.5Q00e+QQQ to S.OOOOe+OQQ [J 6.OOOOe+OOO to 6.50000+000
6.5000e+000 to 7.OOOOe+OOO 7.OOOOe+OOO to 7.50000+000 7.5000e+000 to 7.9203e+000 nterval = 5.0e-001
(b)
t = 10 s
<p= 32°
Contour of Shear Strain Increment
Gradient Calculation
2.5842e-006 to 5.0000e-002 5.0000e-002 to 7.50006-002 7.5000e-002 to 1.00000-001 1 .0QQQe-Q01 to 1.25Q0e-QQ1 1 ,2500e-001 to 1.50000-001 1 .5000^-001 to 1.7500e-Q01 1 .7500e-001 to 2.00000-001 2.0000e-001 to 2.2500e-001 2.2500e-001 to 2.50006-001 2.5000e-001 to 2.75000-001 2.75006-001 to 2.8527e-001 nterval - 2.50-002
t = 25 s
Contour of Shear Strain Increment
Gradient Calculation
■ 7.0697e-006 to 2.0000e-001
■ 2.0000e-001 to 3.0000e-001 H 3.00006-001 to 4.00006-001 fj 4.0000e-001 to 5.0000e-001 M 5.0000e-001 to B.0000e-001
1 6.0000e-001 to 7.0000e-001
__ 7.QQQQe-QQ1 to 8.QQQQe-QQ1
_ 8.0000e-001 to 9.0000e-001 □ S.0000e-001 to 1 .0000e+000 .00006+000 to 1.10006+000 000e+000to 1 .1052e+000 -v’al = 1 .06-001
/ = 40 s
Contour of Shear Strain Increment
Gradient Calculation
HO.OOOOe+OOO to 1.0000e-001 2.00006-001 to 3.00006-001 m 4.00000-001 to 5.00000-001 _ e.Q0QQe-Q01 to 7.00Q0e-Q01 _ S.00000-001 to 9.00000-001 1.0000e+000 to 1.1000e+000 Vj 1 .20006+000 to 1.3000e+000 1,4000e+000 to 1 .S000e+000 1 .50006+000 to 1.57536+000 nterval = 1 ,0e-001
At = Ah,
(u2 + gh' ± 2u^gnj)
(5)
where Amin = min{A1, A2}, n is the time index of the point, and i is the index in the spatial variable.
Figure 5. Distribution of the shear strain during the landslide movement computed for four time moments.
and A1 and A2 are the coefficients of artificial viscosity subjected to the condition
(4)
Using the perturbation method, we obtained the following Figure 6. Model scheme of the submarine landslide move-stability condition for the given difference scheme: ment.
A
Figure 7. Evolution of the surface shape of the submarine slope; the blue and red lines show, respectively, the initial position of the slope and its position at the landslide immobilization time moment.
[16] Figure 7 illustrates the evolution of the surface water wave form above the moving landslide. As is well seen, the landslide generates a dipolar wave with a trough facing the coast and a crest moving seaward. This ap-
pears to be evident because the landslide movement on the seafloor is known to generate an elevation wave above the leading front of the landslide and a depression wave above
10
x
Figure 8. Evolution of the form of the landslide-generated surface water wave at various time moments.
-200 0 200 400 600 800 1000 1200 1400 1600 1800 2000
X
Figure 9. Velocities of particles in the surface water wave at various time moments.
its rear front [Mazova, 2003; Iwasaki, 1997]. The situation does not change if the landslide moves on an inclined surface: an elevation wave arises above the leading front of the landslide, and a depression wave arises above the landslide “tail”, thereby forming a run-down wave from the slope where landslide occurs. The run-down wave increases until the time moment t = 15 s. As is well seen, the length of this depression wave is three times as short as the wavelength of the elevation wave (the first wave crest). Beginning from t = 15 s, an additional crest (the second one in the generated wave) forms at the back front of the elevation wave; with the
further movement of the landslide, it reaches the free water surface (t = 30 s) and continues to build up. At the same time, the receding distance of the depression wave decreases and the trough disappears at t = 40 s, implying that the elevation wave approached the shoreline. At t = 45 s, the height of the crest moving toward the shoreline attains 6 m, whereas the height of the first, seaward crest is no more than 4 m (Figure 8). Figure 9 plots the velocity characteristics of the generated wave at various time moments. It is well seen that the maximum horizontal velocity of water particles in the wave receding from the dry shore does not exceed
Figure 10. Distribution of velocities on the surface water wave profile. An arrow shows the direction and relative value of the velocity at a given point of the wave front.
Figure 11. Evolution of the surface water wave after the termination of generation.
5ms
1 while this velocity is about 3ms 1 in the crest
region of the seaward elevation wave; in the newly formed second crest, the horizontal velocity is directed toward the shoreline and, at the final generation stage, when (the landslide stops), attains 6ms-1. More detailed distributions of the velocity on surface water wave profiles at time moments in the interval t = 15-45 s are presented in Figure 10. As is well seen at t = 15 s, the leading and trailing fronts of the first trough (the run-down wave) move toward each other. This pattern results in the formation of the second crest at the back front of the trough (at t = 25 s). As clearly seen from Figure 10, the velocities of both the back front of the smaller crest and the trough following it are directed toward the shoreline, while the first, larger crest moves seaward (the arrows in the figure). It is seen that, during the immobilization of the landslide (t = 45 s), the surface water wave is discomposed into two wave trains: the first moves toward the open sea, while the second train consistings of the two depression waves and the wave crest moves toward the dry shore (t = 55 s). The second wave train yields a vertical run-up of about 8 m at the isobath x = -200 m. Figure 11 shows the evolution of the shapes of the first and the second crests for time moments in the interval t = 45-65 s (after the landslide movement stops). The wave traveling into the sea retains its form and a height of about 4 m.
[17] Figure 12 presents 10 time moments of the generation of a tsunami wave due to underwater landslide movement, its propagation and run-up on the coast approximated by a plane slope. Computations showed that the heights of wave run-up of landslide origin are up to 6 m and the heights of the run-up at the opposite shoreline of the bay is of the
order of 4 m; these values are well consistent with available observational data.
Conclusion
[18] The data of observations in the Bay of Corinth indicate that, in the place of the landslide origin, water first receded from the shore and then returned, producing a run-up of 5-6 m. The numerical simulation of the tsunami wave runup performed in this study showed that estimates obtained after the generation and propagation of the wave to the isobath x = -200 m agree well with observational data, thereby supporting the validity of the numerical schemes chosen in this study. Moreover, as was shown by Jiang and LeBlond [1993, 1994], who considered viscous and viscoplastic models of a landslide and solved this problem in a long-wave approximation for both the landslide and the generated wave, the generated wave first recedes from the shoreline at the initial time moment of the landslide movement, after which an elevation wave starts forming. Thus, the solution patterns obtained by Jiang and LeBlond [1993, 1994] and in our work are qualitatively similar; however, the analysis of the landslide movement in terms of an elastoplastic model, incorporating real data on sediments at the landslide origin and taking into account the strength reduction in the ground during the development of plastic deformations, provides the most correct (at present time) results of numerical simulation of tsunami wave generation and propagation and run-up values fitting better the available observational data.
Figure 12. Generation, propagation, and run-up of a tsunami wave due to underwater landslide movement computed in the framework of the elastoplastic model.
References
Fine, I. V., A. B. Rabinovich, E. A. Kulikov, R. E. Tromson, and B. D. Bornhold (1998), Numerical modelling of landslidegenerated tsunamis with application to the Skagway Harbor tsunami of 3 November, 1994, in Proc. Int. Conf. on Tsunamis (Paris 1998), p. 211, CEA Press, Paris.
Galanopoulos, A., N. Delibasis, and P. Comninakis (1964), A sea wave from a slump set in motion without shock, Ann. Geol. Pays Hellen, 16, 93.
Garagash, I. A., and V. A. Ermakov (2001), Use of geological-geophysical models for the simulation of the crustal stress state: A case study of Sakhalin Island and the Northern Tien Shan, in Proc. III Scientific Conf. ”Seismicity Problems in the Far East (in Russian), p. 33, IVS FEB RAS Press , Khabarovsk.
Garagash, I. A., L. I. Lobkovsky, O. R. Kozyrev, and R. Kh. Mazova (2003), Generation and run-up of tsunami waves caused by a submarine landslide, Okeanologiya, 43(2), 185.
Iwasaki, S. I. (1997), The wave form and directivity of a tsunami generated by an earthquake and a landslide, Sci. Tsunami Hazards, 15(1), 23.
Jiang, L., and P. H. LeBlond (1993), Numerical modeling of an underwater Bingham plastic mudslide and the waves which it generates, J. Geophys. Res., 98(C6), 10303.
Jiang, L., and P. H. LeBlond (1994), Three-dimensional
modeling of tsunami generation due to a submarine mudslide, J. Phys. Oceanogr, 24(3), 559, doi:10.1175/1520-0485(1994)024<0559:TDM0TG>2.0.C0;2.
Mazova, R. Kh. (2003), Tsunamis generated by submarine landslides, Izvestiya RAEN PMM (in Russian), 4, 117.
Mazova, R. Kh., G. A. Papadopoulos, and L. Yu. Kataeva (2004), The analysis of aseismic tsunami 7 February 1963 in Corinth Gulf: One-dimensional theory, Izvestiya RAEN PMM (in Russian), 9, 63.
Murti, T. S. (1981), Seismic Marine Tsunami Waves (in
Russian), 342 pp., Gidrometeoizdat, Leningrad.
Papadopoulos, G. A., Ed. (2000), Historical Earthquakes and Tsunamis in the Corinth Rift, Central Greece, 128 pp., Inst. Geodynamics, National Observatory of Athens, 11810, Athens, Greece.
I. A. Garagash, Institute of Physics of the Earth, Russian Academy of Sciences, B. Gruzinskaya Str., Moscow, Russia
L. I. Lobkovsky, P. P. Shirshov Institute of Oceanology, Russian Academy of Sciences, 36 Nakhimovskiy pr., Moscow, Russia R. Kh. Mazova, L. Yu. Kataeva, Nizhny Novgorod State Technical University, Nizhny Novgorod, 603600 Russia (raissama-zova@mail.nnov.ru)
