CC BY
13
RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 8, ES5004, doi:10.2205/2006ES000213, 2006
A software complex for the numerical simulation of generation and propagation of tsunami waves in various marine areas from dynamic seismic sources
L. Yu. Kataeva, A. V. Romanov, R. Kh. Mazova, and I. V. Kozhevnikov Nizhny Novgorod State Technical University, Nizhny Novgorod, Russia Received 21 October 2006; accepted 30 October 2006; published 20 December 2006.
[1] A software complex is designed for the simulation of tsunami generation by a dynamic seismic source consisting of several vertically moving blocks of a preset configuration. The program interface allows the user to change the source kinematics of the blocks as a function of the seismic source under consideration. The software makes it possible to display the tsunami wave propagation in a water area (unconstrained by concrete conditions), perform its monitoring, interrupt and resume computations, depict the tsunami wave propagation front, display the distribution of maximum wave heights in a marine area, construct histograms of maximum wave heights along a coastline, and construct the velocity fields of particles in the traveling wave, wave front velocities, and mareograms at any given points of the area under consideration. Using the software capabilities, current data and their variation in time can be obtained at any point of the calculation domain, and all results (both final and intermediate) can be saved in the form of both numerical data and
plots or diagrams. INDEX TERMS: 3285 Mathematical Geophysics: Wave propagation; 4255 Oceanography:
General: Numerical modeling; 4594 Oceanography: Physical: Instruments and techniques; KEYWORDS: program interface, mathematical model, shallow water equations, software product.
Citation: Kataeva, L. Yu., A. V. Romanov, R. Kh. Mazova, and I. V. Kozhevnikov (2006), A software complex for the numerical simulation of generation and propagation of tsunami waves in various marine areas from dynamic seismic sources, Russ. J. Earth. Sci., 8, ES5004, doi:10.2205/2006ES000213.
Introduction
[2] At present, many numerical schemes exist together with the related interfaces designed in such a way that a user insufficiently skilled in programming can work with the program. Developed software complexes are not always conformable to the algorithm used by its author for implementing the concept of a given physical process. Our software complex is based on the original concept of a dynamic seismic source of various configurations moving in various times intervals and at various velocities, and its interface enables the most comprehensive implementation of this idea. As a mathematical model, we used the well-known nonlinear equations of shallow water [Garagash et al., 2003; Lobkovsky et al., 2005], and both the well-known numerical schemes [Bakhvalov, 1975; Marchuk et al., 1983] and efficient difference schemes based on the iterative interpolation method were used for numerical implementation of the model, which
Copyright 2006 by the Russian Journal of Earth Sciences. ISSN: 1681-1208 (online)
significantly reduced the computation time, without invoking methods of concurrent programming.
Main Equations
[3] To describe the generation and propagation of a wave, we use the nonlinear system of equations of shallow water
Ut + U • grad U + g • grad n = F nt +div( (H + n - B)U) = Bt .
(1)
Here grad(s) = Sxi + Syj, div(s) = (si)x + (s2)y, i and j are unit vectors directed along the x and y axes,
/
F
fv - g
\
-fu - g
u\/u2 + v2 Ch2 (H + n - B) v\/u2 + v2
\
Ch2{H + n - B) J
U = (“), f = 2ficos($) is the Coriolis parameter, fi is the velocity of the Earth, 6 is the latitude, g is the gravity
n+1 n g • At
uj = uj - 2TAX
L(nS+1 - nn+i)
n
- f • vij - g
uiW uj + vij2
vn+1 = vn. - g • At Uij Uij 2 • Ayp
Ch2 (Hij + n!j+1 - Bij+1) - (nlj+1 - nn+11j )
v'ij\/ un2 + vij2
[4] The inferred values are then corrected for the Earth’s curvature. The time step was chosen from the stability condition of the difference scheme in use.
[5] For an adequate implementation, we introduced a model grid with spatial intervals (in degrees) Ax and Ay and with a time step At. The spatial intervals were corrected for the Earth’s curvature: Axp = (Ax • n • REarth)/180, where REarth is the radius of the Earth, and the meridian intervals (in meters) between neighboring nodes were found from the relation Aypj = (Axp • n • cos (yn + j • Ay))/180; as mentioned above, the time step (in seconds) was calculated from the stability condition of the difference scheme.
[6] To describe free boundaries in the problem of tsunami wave propagation, we chose the well-known Sommerfeld condition, according to which the part of the wave field beyond the boundary is transferred along the outer normal with a constant velocity defined by the basin depth near the boundary and without changing the waveform. The normal component of the velocity at the free boundary is calculated from the known relation un = cq/(H + n), and the tangential velocity component vanishes, uT = 0.
Figure 1. Central control window.
acceleration (9.8 m s-2), Ch = (H+n-B)0A/sh is the Shesi coefficient, sh is the asperity factor, and B(x,y,t) is the law of motion of the basin bottom. If a point lies in the bottom motion region, the wave height is converted to the increment value with the corresponding sign. The found values of the wave height are then used for updating velocities:
Program Interface and Its Description
[7] The developed software complex can be easily used for modeling the generation and propagation of tsunami waves from a dynamic source in various water areas. Below we describe the possibilities of this complex.
[8] Figure 1 presents the central control window containing six modules:
[9] 1. Transformation of initial data;
[10] 2. Module loading the transformed file and initializing the required data arrays;
[11] 3. Module specifying calculation parameters and display parameters of results;
[12] 4. Module specifying dynamic sources and their characteristics;
[13] 5. Module displaying and saving results;
Ch2 (Hij + nij+1 - Bij+1) '
Figure 2. Example of display with a given grid.
p
Figure 3. Example of the velocity field display.
[14] 6. Computation module.
[15] Below we illustrate the operation of these modules in more detail.
[16] The module of transformation of initial data reads the initial data file, corrects the data, and converts them into the binary format. In addition, the module allows the user to choose a part of the given water area that will be used for calculations.
[17] Module 2 loads the data file into memory; the data are used for calculations and for the initialization of arrays required for the calculations and saving of their results.
[18] Module 3 provides the possibility of picture scaling with the use of the Scaling field. The Time and Out time are, respectively, the final time of computations and the time interval specifying the picture updating rate. The Test h field is used to specify the water depth for test computations. In other words, this program allows the user not only
Block 1 Block 2 Block 3
TimeBegin 15 20 0
TimeEnd 30 45 30
UpH eight 3 9 30
x1" 145 146.1 147.2
y1" 50.1 50.1 50.1
x2* 145 146.1 147.2
y2“ 50.5 50.5 50.5
x3“ 146 147.1 148.2
y3x 53.2 53.2 53.2
x4“ 146 147.1 148.2
y4“ 54 54 54
8locks count: [3 Close
Figure 4. Table specifying the dynamic source.
Figure 5. The display of the bottom topography.
to perform calculations in a given water area with an uneven bottom topography but also to compare their results with the case of an even bottom with the same coastlines, as well as to specify the grid (Figure 2). The VField tool enables the display of the current vector field of velocities with an adjustable step (Figure 3). The Size of diagram field is used to stretch or contract the picture: the size of the histogram window can be altered by changing the corresponding values in the Size of diagram field. The Isobath field specifies the depth of the calculation domain. Whenever this depth is reached, the wave is assumed to strike a vertical rigid wall from which it is reflected.
[19] Module specifying dynamic sources and their characteristics is controlled by the Show blocks table field, which allows the user to specify an arbitrary number of blocks. Each block is characterized by the initial and final times of its uplift, the uplift height, and the four coordinates of the block vertices, which can create a great diversity of block motion patterns (Figure 4).
[20] Module displaying and saving results enables the following functions.
[21] (i) The Show bottom tool enables the display of the bottom topography. The result can be saved in the bmp format (Figure 5).
[22] (ii) The Show Terr tool enables the display of the bottom topography and the dynamics of block motion. The result can be saved in the bmp format (Figure 6).
[23] In addition, the following functions of automatic saving of pictures in the process of computations are implemented.
Figure 8. The sample of Show data window tool.
Figure 6. The bottom topography and the dynamics of block motion.
[24] (iii) The Show max H and Show min H tools display the current maximum and minimum values of wave heights in the water area. These tools are active only if the process of computations is in progress or interrupted. The result can be saved in the bmp format.
[25] (iv) The Show Eta tool displays current values of wave heights in the water area. The tool is active only if the process of computations is in progress or interrupted. The result can be saved in the bmp format. If the auto option is inactive, results are not saved.
[26] (v) Several settings of file saving are available (Figure 7).
[27] (vi) Show data window can be active only during the
Figure 7. Settings of the file saving. (1) Save together with the vector field and grid; (2) save together with the vector field; (3) save together with the grid; (4) save without the grid and vector field.
computation process. If active, this tool displays a floating popup window showing the following values: h, ocean depth at a given point; u and v, velocities of wave particles at the given point along the X and Y axes, respectively; vis, current wave height at the given point; max vis and min vis, maximum and minimum ocean surface heights at the given point over the entire time of computations; C, wave front velocity (Figure 8).
[28] Thus, we have complete information on the water area throughout the computation process.
[29] (vii) The Show points tool is used to specify points for the analysis of the sea level oscillations caused by the passage of the tsunami wave (Figure 9). If the Show flag is set at unity, the plot at a given plot will be displayed. In all, 35 points can be selected for the analysis. Their number is controlled by the Number of points option. Any changes introduced in this window are ignored until the OK key is hit; otherwise, all unsaved changes will be lost upon the closure of the window. This window is very convenient because it enables the analysis of the ocean level evolution at each point of interest.
[30] (viii) The Show line form tool is used to specify the line along which the tsunami mareogram is constructed (Figure 10). The end points of the line and the desired time interval can be set using the tool. This changes the wave height along the chosen line. The results are entered into a table and are saved in a file.
[31] Computation module provides for the possibility of a new run, its saving, its interruption and renewal, and the test run, i.e., the computation in the same water area with a fixed depth value.
Figure 9. Specification of points for the analysis of the sea level oscillations.
Figure 10. Specification of the line along which the tsunami mareogram is constructed.
References
Bakhvalov, N. S. (1975), Numerical Methods (Analysis, Algebra, and Ordinary Differential Equations) (in Russian), 632 pp., Nauka, Moscow.
Garagash, I. A., L. I. Lobkovsky, O. R. Kozyrev, and R. Kh. Mazova (2003), Tsunami wave generation and
runup caused by a submarine landslide, Okeanologiya (in Russian), 4-3(2), 185.
Lobkovsky, L. I., R. Kh. Mazova, L. A. Garagash, L. Kataeva,
and N. S. Petrukhin (2005), Analysis of the 26 December 2004 earthquake and tsunami in the Indian Ocean on the basis of the subduction keyboard model, EGU General Assembly, EGU, Vienna, Austria, 24-30 April 2005.
Marchuk, A. G., L. B. Chubarov, and Yu. I. Shokin (1983), Numerical Simulation of Tsunami Waves (in Russian), 175 pp., Nauka, Novosibirsk.
L. Yu. Kataeva, I. V. Kozhevnikov, R. Kh. Mazova, A. V. Romanova, Nizhny Novgorod State Technical University, Nizhny Novgorod, 603600 Russia ([email protected])
