CC BY
16
RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 8, ES6004, doi:10.2205/2006ES000216, 2006
Principles of numerical modeling applied to the tsunami problem
Yu. I. Shokin, L. B. Chubarov, Z. I. Fedotova, S. A. Beizel, and S. V. Eletsky
Institute of Computational Technologies, Siberian Branch of the RAS, Novosibirsk, Russia
Received 21 November 2006; accepted 20 December 2006; published 24 January 2007.
[1] The main principles of using computational tools in the solution of applied tsunami problems related to support the decisions of antirecessionary manages are formulated.
A comparative assessment is given of the efficiency of known computational models and programming systems developed on their basis for hydrodynamic tsunami modeling. Index
TERMS: 4564 Oceanography: Physical: Tsunamis and storm surges; 3285 Mathematical Geophysics: Wave propagation; 4255 Oceanography: General: Numerical modeling; 7209 Seismology: Earthquake dynamics;
KEYWORDS: tsunami generation, seismic source, wave propagation.
Citation: Shokin, Yu. I., L. B. Chubarov, Z. I. Fedotova, S. A. Beizel, and S. V. Eletsky (2006), Principles of numerical modeling applied to the tsunami problem, Russ. J. Earth. Sci., 8, ES6004, doi:10.2205/2006ES000216.
1. Introduction
[2] In our paper, we formulate the principles of using computational tools to solve applied problems of tsunami oriented to the support of decisions by managers in critical situations [Shokin and Chubarov, 1999].
[3] From hereon we shall use term “tsunami” to denote the entire complex of catastrophic impacts caused by the waves in the basins and adjacent territories of the coast. The nature and spectral characteristics of these waves are usually associated with long surface gravity tsunami waves.
[4] Similarly to other natural and anthropogenic catastrophes, tsunamis pass three stages in their development: stage before crisis, crisis, and stage after crisis.
[5] Each of these stages is characterized by its duration, spatial scale, localization, and specific impact on the population involved in the crisis. Each stage requires specific methods of control and formulation of the list of problems, the important place in the solution of which is occupied by the methods of computational modeling. Specific limitations are imposed on mathematical models, numerical algorithms, technologies of calculation, character and volume of the results, means of their presentation (including visualization) means, protocols, and addressing of their spreading.
[6] Leaving a number of provisions formulated above without consideration for a certain period of time we shall discuss the requirements, which the computational tools should satisfy. These tools are allowed by the scientific community for determining the vital characteristics of tsunami that are impossible or difficult to obtain using alternative sources (field
Copyright 2006 by the Russian Journal of Earth Sciences. ISSN: 1681-1208 (online)
studies, laboratory experiment, historical search, analytical methods, etc.).
[7] These requirements should determine the technologies of developing mathematical models, creating numerical algorithms, methods of processing the numerical results, field and historical data, means of visualization and user interfaces, their programming realization, testing of mathematical models, algorithms, and programs.
[8] It is assumed that the recommendations of the international community would be exclusively granted to those systems whose development, testing, documentation, and guidelines of use would satisfy the coordinated requirements.
[9] Such approach is not new. It operates successfully for a long time in the development of algorithmic and programming software in the fields related to the solution of problems influencing the health, safety, and life of people, reliability of creating and application of high risk constructions (medicine, avionics, navigation, etc.).
[10] Let us formulate the list of applied problems solved during each of the stages of tsunami.
• Stage before tsunami:
1. Zoning of the coast by the risk of hazard.
2. Forecast of tsunami probability in time and space.
3. Construction and improving of observational networks: seismic and geophysical.
• Crisis stage:
1. Processing of information obtained from observational networks (monitoring).
2. Operative calculation of arrival times and characteristics of dynamic impact.
• Stage after the crisis:
1. Archiving of data of tsunami observations.
2. Processing of information obtained from all sources.
3. Repeated solution of the problems of the previous stage with account for the newly accepted information with the objective of filling the databases with more detailed data.
4. Formulation of test problems caused by the passed crisis, which are intended for improving the calculation tools.
5. Calculations to perform the orders of persons involved into the control of the crisis situation to optimize the process of crisis liquidation and transition to restoring operations.
6. Preparation of specialized analytical reports related to estimating the specific experience of performing calculation works.
[11] The context of the problems listed above is determined, first, by the mandatory list formulated and accepted by the community of experts, and second, by additional list prepared by the crisis managers involved in specific work with account for the regional and local peculiarities of protected territories and specific crisis situation.
[12] A set of mathematical models, algorithms of calculation and programs is defined depending on the peculiarities mentioned above and formulation of the problems. Two extreme technologies of calculations, and the third one, which is a symbiosis of the first two are possible.
[13] 1. Calculations are carried out locally by the personnel using software tools, which were obtained in advance and thoroughly documented. The personnel were specially trained and received the required certificates. This method has obvious advantages and not less obvious disadvantages.
[14] Advantages. The work is carried out locally. Some parameters of the problems, models, and algorithms can be determined in a flexible manner with adjustment to the observed physical conditions. The major part of the customers is located locally, which decreases the load on telecommunication systems.
[15] Disadvantages. The necessity to perform extremely complicated work in the conditions of chaos and mess, which appears in the crisis situations. Impossibility to enlist the services of the specialists in mathematics, hydrodynamics, computational algorithms, visualization and interpretation of the results including the authors of the applied mathematical models, algorithms and programs. Impossibility of using powerful computational devices increasing the accuracy of calculations, their properties of detailed work, and speed of solving the problems. Unreliability of telecommunication channels needed to obtain large volumes of information (bathymetry, topography of the land surface, characteristics of underlying surface), and wide spreading of the results.
[16] 2. Solution of the problems including processing of the results at special centers certified by the community for performing the corresponding works. The advantages and disadvantages of this approach are exactly opposite to those listed for the first approach.
[17] 3. Finally, the third, combined approach, when different problems are solved in different places. Realization of this approach requires certain additional work to structure the problems and reflect this structure to the corresponding structure of specialized centers of calculation. Good perspectives are opened for optimizing the works, information flows, and databases. The general efficiency of the works increases.
[18] Now, let us discuss the computational tool, which is a result of the interaction of several objects:
[19] Mathematical model with a set of parameters (hydrodynamic: initial data, bathymetry, form of the boundaries, topography of the land, roughness of the surface, wind tension, coefficients of turbulent viscosity, etc., and mathematical values of “virtual” properties, and those specified in deduction of equations);
[20] Computational algorithm, which approximates the equations of the mathematical model with its parameters of spatial and temporal discretization, scheme and artificial viscosity and dispersion, with “additional” boundary conditions, etc.;
[21] Algorithm of preprocessing, which modifies a number of hydrodynamic parameters for their adjustment to the peculiarities of the computational algorithm (recalculation to alternative grids, modification of boundaries, interpolation, etc.);
[22] Algorithm of calculation controlling, which monitors irregular situations during the calculations (especially during long-term modeling) such as loss of stability and others, and introduces necessary changes (automatically or by means of the user interface) in the parameters of the computational algorithm (up to changing the model and algorithm) and provides the specified accumulation of the results for further processing and interpolation as well as ensures operative presentation of critically important information about the dynamics of the modeled phenomenon or computational process, etc.;
[23] Algorithms of post processing, which process the results of modeling, calculated the necessary functionals of the solution, provides visualization of high quality, and if necessary transmits automatically generated messages of necessary volume and character to the specified telecommunication channels. The same algorithms perform selection and sorting of the information for further analysis, its archiving and transmission for long-term storage to the information databases. Most likely, this list can be supplemented, and no doubt presented in more detail.
[24] The problems of application of the methods of computational modeling, appearing in the solution of applied tsunami problems are related to the peculiarities of the corresponding software tools.
2. Formulation of the Problems of Computational Modeling of Tsunami
[25] The formulation of the problem of computational modeling should include clear description of the objective of modeling (character of the problem to be solved), definition of the degree of accuracy of details, indication of the list of required functionals, form of presentation and volume of the results including those, which would be transmitted to the databases.
[26] Taking into account the aforesaid division of applied problems into preliminary, operative, and subsequent ones, the customers, formulators of the problem, and users change correspondingly. The types of the problems and functions of each of them should be determined in the corresponding regulating documents.
[27] Preliminary problems usually require preparing a list of protected basins and objects, tying of the corresponding input data, preparing digital charts (data arrays) satisfying the agreed standards, etc.
[28] The problems should be solved by means of performing a technological series of operations, which includes, as the first stage, test calculations using different grids (with different resolution of bathymetry), mathematical models (to evaluate the importance of nonlinearity, dispersion, dissipation, account for the spheroid form of the Earth, its rotations, etc.) computational algorithms, programming systems because the peculiarities of realization of algorithms sometimes have strong influence on the results. Only after comparing and thorough analysis of the obtained information including the available field and historical data it is possible to choose a model (actual equations, algorithm, and programming code), which can be used at the further stage to obtain the results, which would be used in practice.
[29] Local authorities could be the customers of preliminary problems (national, regional, or local), as well as extraordinary services of different agencies (including international or private), humanitarian organizations and individuals in the part related to their safety and property.
[30] Formulator of the problem should be qualified at a level needed for clear understanding of the goals of the planned computational experiments, their capabilities and limitations, and should have an access to the data needed for correct formulations of these problems.
[31] It is optimal that the functions of calculator would be performed by a specialist or organization, which has specialists in the field of computational modeling of wave hydrodynamics problems that participated in developing of at least one of the applied algorithms for the solution of a specific problem or in their programming realization. Frequently, this function is performed by a professional for adjacent fields (seismologists, geophysicists, employees of extraordinary services, etc.) but in any case it is extremely important to engage in the work (even as a consultant) a mathematician specializing in computational experiments.
[32] The user of the problems (frequently this is the same person as the formulator of the problems) should be qualified enough for adequate interpretation of the modeling results,
systematization of the obtained results, and recommendations for their use.
[33] The operative problems are usually solved in the conditions of limited time resource by the specialists of “limited” qualification, thus their formulations of the problems should be determined in advance as much as possible, while the solution and processing of the results should be automated to the greatest degree. Approbated models, algorithms, and programming codes should be used for the problems of such class. The procedure of approbation is not discussed here as well as the other problems related to the preparation and admission of the regulating documents.
[34] The characteristics of the subsequent problems are very similar to the characteristics of preliminary problems. Their specific features, to our opinion, are significant widening of the users by means of spreading of the results among the maximally possible number of applied specialists, investigators, and interested people.
3. The Choice of Mathematical Models. Some Thoughts Using the Example of Tsunami Wave Transformation Problem
[35] The choice of the adequate mathematical model depends on the specific properties of the solved problem, on the specific features of the basin, in which the problem is solved, and what is not less important, on the stage of the development of the phenomenon within one problem and one basin.
[36] Application of the most complete models is not always reasonable not only owing to a serious increase in the required computational resources and complications of the needed algorithms, but also due to the impossibility of adequate determination of all input parameters of such models. The correct choice can be made only on the basis of preliminary calculations close to the content of the applied problems using different models and different algorithms at different grids. Such work should be done during the preliminary stage using agreed test problems and “real” input data. Finally, zoning of the protected territories of the basin and coast can be made on the basis of mathematical models capable to provide adequate results for each type of the problems.
[37] Zoning by the available types of input data (determination of the grid size for providing the needed accuracy and degree of modification of the boundaries) should be done simultaneously with the above described work as well as zoning by other hydrodynamic parameters: roughness, characteristic wind stress, etc.
[38] The experience of the authors of this article based on the solution of research and applied problems of tsunami allows us to state that even the simplest mathematical models make possible obtaining adequate estimates in the initial stage of tsunami development, while the subsequent effects require a thorough work for their reproduction.
[39] An example of such kind is the investigation of tsunami transformation in “wash-tub” (see Figure 1) model region,
Figure 1. Bottom topography of the “wash-tub” calculation region.
which was introduced for the first time by [Gusyakov and Chubarov, 1985, 1987]. The convenient features of this basin is explained, in particular by the possibilities to study the peculiarities of different mathematical models and algorithms to reproduce the transformation of waves in the regions with different depths and wave interaction with different types of boundaries (open and rigid reflecting boundaries). Regardless its simplicity, the bottom topography in the “wash-tab” region agrees well with the distribution of depths along the Kuril-Kamchatka Trough. The bottom topography is specified by piecewise linear function, which varies from 10 m to 9000 m and depends only on coordinate y (it is assumed that coordinate system XOY is located so that OX axis is directed along rectilinear coast, and OY is directed normal to the coast in the direction of depth increase). The length of the basin in OX direction is 555,000 m and in OY direction is 320,000 m.
[40] Test problems solved in this model basin should facilitate answering the questions about sensitivity of real problems to the computational accuracy of the algorithms and to hydrodynamic accuracy of mathematical models: the two main problems of numerical modeling. Two problems were considered with this goal in mind:
[41] The first of them is propagation of a solitary wave homogeneous in OX direction, and the second, is the transformation of initial perturbation of the free surface finite in both directions. Below, all linear sizes are given in meters and time is in seconds.
[42] In the first problem, the form of the initial perturbation was specified by a well known relation:
n = a • cosh-2(Y) ,
Y / 3a ~(y - yo)
\j 4(H + a) H ,
y0 = 175,000, a = 1, u(t = 0) = 0. Conditions of full reflection were specified at the coastline (y = 0) and at lateral walls (x = 0 and x = 555,000); boundary y = 320,000 was
free. At time moment t = 0, the wave decomposed into two, one of which propagated in the direction of increasing coordinate, and the second propagated in the shallow water region.
[43] Nine calculated gauges were located along line x =
277,500 for detailed study of wave regimes. The following values of the y coordinate were chosen: yui = 0,
yu2 = 1250, yu3 = 2500, yu4 = 5000, yu5 = 10,000, yu6 = 47,500, yu7 = 112,500, yus = yo = 175,000, yug = 230, 000 (see Figure 2).
[44] Modeling was performed within classical linear and nonlinear (NL) models of shallow water [Stoker, 1959], weakly nonlinear model of Nwogu [Nwogu, 1993] and fully nonlinear dispersion (NLD) model [Wei and Kirby, 1995] (“one-layer” model of Liu-Lynett [Lynett and Liu, 2004]). Owing to the fact that the results obtained using the two latter models appeared almost identical, the graphs presented bellow represent the Nwogu model. This fact allows us to state that the effects of nonlinear dispersion terms are of low importance in the solved problems.
[45] Multi-parametric versions of differential schemes of MacCormac [Fedotova, 2006] and Adams [Wei and Liu, 1995; Lynett and Liu, 2004] were used as computational algorithms realizing the first two models. The Adams scheme was used for approximating the Nwogu model.
[46] The calculations using different hydrodynamic models were performed using the Adams scheme (with numbers of nodes in the direction OX — Nx =51 and in the OY — Ny = 513 direction). Records of gauges at the ninth (closest to the source) and fourth (most distant from the source) points calculated from linear and non-linear equations of shallow water and weakly non-linear dispersion Nwogu model are shown in Figure 3.
[47] The graphs in deep water part of the region (Figure 3a) show that all models give practically the same results. Nonlinear effects begin to manifest themselves when the wave enters the shallow water zone (Figures 3d, 3e), in which the linear model differs from the other models in the
■10000]--------1--------1---------1--------1--------1--------r
Figure 2. Location of gauges in problem 1.
Figure 3. Pressure gauge records calculated at the 9th (a) and 1st gauge points using linear (markers) and nonlinear models of shallow water and the Nwogu model.
description of the reflected wave. Pressure gauge records calculated in the shallow water region, beginning from the fourth demonstrate increasing difference of the linear model results also in the description of the leading incident wave. Almost complete coincidence of the results of nonlinear equa-
tions of shallow water in both NLD models at all points of gauge calculations points to the fact that dispersion effects are insignificant in the problem with such initial conditions.
[48] In the second problem, initial perturbation of the free surface was specified in the form of exponential “caps” of different length A using relation
n = a exp I —
x — x0 L
exp —
y — yo L
Figure 4. Location of gauges in problem 2.
where L = 16,000 for A = 50,000, L = 24,000 for A = 75,000 and L = 32,000 for A = 100,000; a = 2 (the length is determined on the basis of 10% section of the initial amplitude a from the foot of the wave). The second problem was solved within classical nonlinear model of shallow water [Stoker, 1959] and weakly nonlinear dispersion Nwogu model [Nwogu, 1993].
[49] At time moment t = 0, the center of perturbation was located at point with coordinates (x0,y0) = (465,000; 230,000). Condition of full reflection was specified at boundary y = 0; the other boundaries were assumed free. Twenty three gauges were located to record the results of calculation. A scheme of their location is shown in Figure 4.
[50] Spatial images of wave pattern (from above) for three time moments are shown in Figure 5. Black horizontal lines denote breaks of bottom topography lines.
2
2
Figure 6. Wave patterns for initial perturbations of different length: First column: A = 50,000, second column: A = 75,000, third column: A = 100,000; First line: t = 0, second line (NL): and third line (NLD): t = 900, fourth (NL): and fifth line (NLD): t = 1800.
[51] Waves propagate faster in the deepest (“upper”) part of the basin; break lines are characterized by the changes in the curvature of fronts (reflection). The wave reflected from a break line with the greatest angle (a ~ 14°) is clearly seen.
[52] Dispersion effects manifest themselves stronger in the case of “short” initial perturbation demonstrating the generation of not only one wave as in the case of nonlinear equations of shallow water but the appearance of wave packets propagating in different directions from the source. This effect manifests itself to a lesser degree for the perturbations of medium length. Finally, for the “longest” perturbation the dispersion trace appears practically invisible. This is confirmed by the records of gauges (see Figure 6). Wave packets, which follow the leading wave, are clearly seen in
Figure 5. Pressure gauge records at 11th (first column) and 16 (second column) points, A = 50, 000 (a), (b); A = 75, 000 (c), (d); A = 100, 000 (e), (f). Thin line denotes nonlinear equations of shallow water; heavy line denotes the Nwogu model.
Figures 6a,b. If the length of the initial perturbation increases, these wave packets gradually disappear.
[53] The character of free surface oscillations, which is seen in the gauge records in Figure 6, allows us to note the differences in the wave regimes at coastal points (11th gauge) and in the offshore direction (16th gauge). A simple signal is observed at the coast consisting of the leading wave and an elevation behind the wave (dispersion decomposes this fragment into a series of oscillations). At a distance from the coast, the leading wave is followed by an elevation (wave packet in dispersion approximation) and a wave reflected from the boundary of the underwater slope.
[54] A calculation for a longer period of time was performed to study the reflection process in the case of A = 50, 000 m. Pressure gauge records calculated at the coastal (9th) point, at a point close to the coast (11th), at a point close to the source (12th), and remote point (16th) are shown in Figure 7. These graphs demonstrate (9th point) that the existence of a lengthy shallow water interval in the pathway of the wave propagation leads to the appearance of a peculiar filter, which removes small scale dispersion oscillations and allows us to reproduce the coastal wave processes within the equations of the shallow water theory with the same efficiency. At the same time, the gauges located in the zones of greater depth (11th, 12th, and 16th gauges) demonstrate the necessity of the account for the dispersion in reproducing the initial phases of the wave regime already, whereas its final phases determined by the waves reflected from the break lines of the bottom located in deep water zones do not require application of NLD models. It is noteworthy that the amplitudes of the waves propagating over maximal depth (16th gauge) are almost five times smaller than the amplitudes of the waves propagating in the region of decreasing depth.
[55] The observation of the peculiarities of wave propagation along different pathways characteristic of the basin studied here is very informative. In this case, the material for such analysis is observation of the wave transformation when it propagates from the source (21st gauge) to the coast (9th gauge). This analysis characterizes the capabilities of the mathematical models and algorithms to reproduce the transformation of waves over the peculiarities of bottom topography, whose depth varies from 9000 m in the region adjacent to a depth of 10 m near the coast. A comparison of gauge records showed that the form of the leading wave does not change practically, its amplitude decreases strongly after it leaves the source but in the course of the further propagation it conserves its values, while the form of the wave packet that follows the leading wave becomes more complicated due to multiple interactions with the bottom topography and break lines of the bottom. The amplitudes of the waves in the wave packet increase as the depth decreases. We note that “hydrodynamic” dispersion caused by the presence of the corresponding terms in the equations of the mathematical model manifests itself only in a slight decrease of the amplitude of the leading wave and has no effect on its phase characteristics. Thus, nonlinear and nonlinear-dispersion models lead the similar results.
[56] The amplitude of the leading wave decreases strongly during the propagation over the deep water part of the basin.
0.6
Figure 7. Pressure gauge records at 9th, 11th, 12th, and 16th points. Thin solid lines denote nonlinear equations of shallow water; heavy dashed line denotes the Nwogu model.
Its form is strongly determined by the influence of dispersion, which is practically not seen during the subsequent stage of the process. This part of the results characterizes the capabilities of mathematical models and algorithms to describe the wave transformation in the deep regions located far from the boundaries.
[57] Wave regimes calculated at coastal points are very complicated. The displacements of the free surface at these points are calculated using the algorithm approximating the corresponding boundary conditions. They are determined from the assumption of the coast as a vertical wall. The analysis of calculations near the “reflecting” boundaries showed that “hydrodynamic” dispersion does not manifest itself, and the amplitude of the waves has a tendency to decrease, but not to such degree as over the “marine” pathway. The character of oscillations and their frequency properties do not practically change from one point to another.
[58] The calculation of one of the most important characteristics, the distribution of maximal amplitudes along the coast as seen in Figure 8, does not require enhanced mathematical models, at least for the specific set of “hydrodynamic” parameters.
[59] The provisions formulated above are confirmed by a series of spatial images shown in Figure 9 and Figure 10. The first of them demonstrates the results obtained using the Nwogu model, and the second using nonlinear model of
Figure 8. Maximal wave heights along the coast calculated using nonlinear model of shallow water (thin line with markers) and using the Nwogu model (heavy line).
shallow water. The figures are shown in traditional lexicographic sequence (from left to right and from top to bottom) beginning from time step 500 with an interval of 500 steps. Significant differences are found only in the first four figures, the next four differ significantly less, and beginning with the ninth figure the states of the surface in the basin are practically not distinguishable. Cell structures appearing during the interaction between waves propagating to the coast and waves reflected from the coast and from the break lines of bathymetry are worth attention. The account for the dispersion strengthens this effect especially in the zone of the flat shelf and adjacent regions of the basin.
4. Requirements to the Computations Algorithms in Tsunami Pproblems
[60] The problem of sensitivity of applied tsunami problems to hydrodynamic accuracy of mathematical models and to the computational accuracy of discrete algorithms that realize the models becomes the first in the numerical solution of such problems. In the first part of this study we used the
Figure 9. Dynamics of the free surface calculated using nonlinear-dispersion model.
Figure 10. Dynamics of the free surface calculated using nonlinear model of shallow water.
example of the problem of wave transformation in “wash-tub” model region, which was first considered by [Gusyakov and Chubarov, 1985], to analyze the peculiarities of different mathematical models reproducing the wave transformation in the regions with different depths and wave interactions with different types of boundaries (open and rigid reflecting boundaries).
[61] In the second half of the work, we study the problem of choosing the computational algorithms, which provide the required accuracy of the modeling results and formulate a number of test problems for estimating the adequacy of the developed algorithms and programs.
[62] Computational algorithms used in the problems of tsunami modeling can be divided, for example, into the algorithms designed for the solution of differential equations, algorithms of data processing (field, historical, etc.), and others. Here, we shall consider the algorithms of the first group based on finite difference methods. They are designed for calculation of wave amplitudes, velocities or transports, characteristics of the interaction with the coast, fields of maximal heights, and other functionals.
[63] The main requirements to the algorithms are directed to reach the needed accuracy, economical properties, and flexibility. Sometimes these requirements correlate well with each other but sometimes hardly solved contradictions arise. Such properties as stability margin, conservatism, invariant features, monotonicity, and lack of high frequency fluctuations facilitate gaining the satisfied accuracy in the need of using rough spatial discretization. Due to the approximation, computational algorithms inherit the parameters of the mathematical models (“hydrodynamic” and “mathematical”) and introduce their own “computational” parameters, among which the most important are discretization parameters (size of computational cells, basic functions), which determine to a great degree internal properties of the algorithm, as well as additional parameters, which are used to control the scheme dispersion and dissipation, and other properties, in particular those mentioned above. We should also note the situation, which appears when asymmetric algorithms are used (multi-step splitting methods, methods of the predictor-corrector type, etc.) when special alternation of individual stages of the algorithm is required to avoid
this asymmetry. In this case, symmetry is gained not at each time step.
[64] A transition to a discrete region of variation in spatial variables and time can be performed using different methods: by introducing a regular of irregular grid, which can be uniform or not uniform, either triangle, rectilinear, or curvilinear, fixed or varying in time, adjusted to the geometry of the calculation region or approximating it to a certain degree. The specific geometrical feature of calculation areas in tsunami problems, which reflects the specific features of the corresponding real basins, manifests itself in their multiply connected property and extremely complicated form of internal and external boundaries determined by the complexity of the coastline contours. In some situations, when the formulation of the problem requires the computation of the zone of inundation, such boundaries become mobile, which significantly complicates the problem. This led to the development of algorithms, which allow for the so-called continuous computation without distinguishing the water-land boundary (a type of the method of dummy regions).
[65] The efficiency of the computational algorithm depends strongly on the correctness of specifying the numerical boundary conditions. The conditions that correspond to the physics of the studied phenomenon are formulated at the moment when the mathematical model is being chosen for reproducing such processes as reflection from a wall, free pass, runup of waves on the shore including the moving one (landslide above water), or arrival of the wave into the computational area. It is desirable to approximate these “physical” boundary conditions so that the order of the accuracy of the algorithm as a whole is not changed. Usually, the specific property of the numerical method requires formulation of additional conditions, whose number and form are determined jointly by the algorithm and mathematical model. A theoretical investigation of the correctness of computational boundary conditions is possible only in the simplest cases, while in the solution of specific problems their experimental study is emphasized, which requires thoughtful tests and a large amount of computational experiments.
[66] Maintenance of the required accuracy of calculations can be gained using two methods, at least. The first of them is related to the use of more detailed input information. A quantitative measure in this case would be a spatial size, for example the size of the grid if regular discretization is used. The key problem here is the accuracy of input information, which consists of the sets of bathymetric data (digital charts), data on the roughness of the bottom and land topography, etc. A transition to a more detailed bathymetry introduces the effect of increasing accuracy, but here, one has to pay attention to the fact that more detailed charts of real basins are the result of applying certain interpolation algorithms to the initial data measured at the final number of points. Application of such procedures is inevitably accompanied by the accumulation of errors. Uniform decreasing of the grid in the entire computational area increases the resource capacity of the calculation and very frequently is not needed. An alternative to this approach is application of non-uniform grids of different type when condensation that provides increased accuracy of the results is performed, for example, in the zone of shallow depths or in the zone of their
sharp variations. Such solution of the problem is accompanied by certain difficulties, for example by the fact that the computational area is multiply connected or by the behavior of the solution at the boundaries between the fragments of the computational areas of different size of the grid.
[67] Application of the algorithms of increased order of approximation can be considered as the second method of increasing the accuracy. As a rule, the choice of computational methods of such class complicates significantly the algorithm, increases the number of points involved in the calculation, and complicates the formulation and realization of the boundary conditions. This aspect of providing the accuracy of calculations requires thorough grounds and detailed study to make a reasonable solution.
[68] Increasing of the accuracy of numerical modeling can be gained also by more accurate specification of the boundaries in the calculation area. This approach becomes most pressing in the case, when the objective of modeling is determination of the characteristics of wave regime at the boundary and near boundary points. Different methods exist to gain such increase in the accuracy, for example, by application of curvilinear grids adjusted to the contour of the boundary or by transition from traditional finite difference schemes to the methods of finite element type, finite volumes, and others, which use triangle or irregular quadrangular grids. Certain difficulties for the finite element methods in tsunami problems are related to the reproduction of a number of characteristic boundary conditions.
[69] The problem of providing the computational stability in tsunami problems has its own specific features. First of all, one should differentiate the local and global problems in time. As a rule, stability of calculations can be gained by standard methods in the solutions of local problems designed for modeling short-term processes occurring in limited basins (first of all by using theoretically stable methods that have reliable stability margins and mechanisms, which filter out non-physical high frequency oscillations). Longterm modeling, in the course of which the wave can propagate over transoceanic distances, requires additional software tools, permanent monitoring of the state of calculated fields with the objective of early distinguishing of potential sources of instability and their timely removal either automatically or as a result of interactive interruption of computational expert. The development of instability during long calculations is usually related to the so-called “nonlinear instability”. The time of its occurrence and spatial localization can be hardly performed in advance or this is practically impossible. In the test problems considered by the authors we made an attempt to correlate these effects with certain peculiarities of bottom topography. However no universal criterion has been developed yet.
[70] Stability of “short term” calculations depends strongly on the smoothness of bottom topography in the area considered in the study and smoothness of the initial data, which are input data. They are determined, for example, using other computational algorithms or as a result of processing of field measurements. The computational experience of the authors indicates that such dependence manifests itself to a certain degree when mathematical models are used, whose equations include high order derivatives (non-linear
dispersion models, account for elastic effects at the free surface in modeling the interaction of waves with floating platforms, airdromes, etc.). The problem of smoothness of the input information can be solved using special smoothing algorithms, which to a certain degree of sophistication include the mechanisms of artificial dissipation. If necessary, these algorithms can be used to remove high-frequency computational oscillations of the grid scale. Such work requires high computational qualification. If the work is performed by a dilettante, it frequently leads to the attempts to interpret purely computational effects as physical ones. An example of such type of errors would be given below. We note, however that the determination whether one or another effect in the results of numerical modeling is “physical” requires that each problem should be solved using sequentially decreasing grids. Only conservation of the found effect at any scale of discretization allows us to start interpreting it as an existing phenomenon. Naturally, we have to mention the correct choice of the time step and a possibility of its variation from one layer to another in the solution of nonlinear problems with attempts to maintain the “Courant number” as close to unity as possible. The solution of the problem of computational stability (prevention of the collapse of numerical solution) in the case of long-term calculation requires special thoroughness in the formulation of the modeling problem. For example, the mechanisms of smoothing and artificial dissipation, or others mentioned above could be an admissible means of providing stability if it is needed to reproduce the long-term dynamics of the most general characteristics of the phenomenon and determine general regularities in the variation of wave parameters in space and time. A different case is in the situations when it is required to determine absolute values of wave heights, velocities, and runup values. Such calculations should be preceded by thorough methodological computational experiments over real topography using the grids of different size, various computational algorithms, and different stabilization mechanisms.
[71] The considerations described above would be illustrated on the examples related to the solution of a number of test problems including those considered in the previous sections about wave transformations in the “wash-tub” basin.
5. Analysis of the Accuracy of Algorithms on the Example of the Problem of Tsunami Waves Transformation
[72] The objective of calculations, the results of which are described below, was the study of different approaches to provide the accuracy of algorithms. With this goal in mind the calculations were performed using the methods of different order of approximation and grids with different resolution. It was already shown in the first half of this paper that dispersion effects in the problem of the propagation of a solitary wave in the “wash-tub” basin appear insignificant. Thus, the further consideration is based mainly on the material obtained using classical equations of the shallow water theory.
[73] The calculations were carried out using the schemes
of MacCormac and Adams and a sequence of grids with sizes decreasing in the OY direction (the number of points was equal to Ny = 129 (Ax = 2h0), Ny = 257 (Ax = h0), Ny = 513 (Ax = h0/2), Ny = 1025 (Ax = h0/4)). The first scheme has the second order of approximation at internal points of the area. The condition of reflection at the wall (at the first node) was approximated by the first order. As to the Adams scheme, the authors used two versions. In the first of them at the second node and at (N — 1) node located near the boundary over the OY-axis, the first spatial derivatives in nonlinear terms of equations were approximated using the second order (central difference calculated over three points), while at the other internal nodes, the equations were approximated using the fourth order based at five points. In the second version the fourth order of approximation at near boundary nodes (using a five-point template displaced to the interior of the area) was conserved. The boundary condition at the first node in both cases was approximated using the second order.
[74] First of all, we shall describe the general characteristics of this quasi-one-dimensional process, which are shown in Figure 11 as gauge records calculated along the pathway of wave propagation to the coast (gauges 1-7) and in the seaward direction (gauges 8-10).
[75] A solitary wave with an amplitude of 1 m, which at the initial time moment was located over gauge 8, splits into two. One of them (“right”) propagates to increasing depths and after passing gauges 9 and 10 leaves the basin. The second (“left”) propagates to the coast. The “right” wave almost immediately interacts with a sharp break line of the topography, which leads to a distortion of the trailing front of the leading wave, which is manifested at gauges 9 and 10. In the final fragments of these records one can see a flat wave generated by the wave transformation over the shelf area (depth is 1000 m) reflected from a flat coastal slope and from a steep break line of topography at the end of the shelf area. Their interaction is clearly seen in the record of gauge 6, while the record at the next point (gauge 7) shows the separation of the wave mentioned above, which propagates in the seaward direction.
[76] Coastal gauges demonstrate the interaction of the incident wave with the coastal wall. During this interaction, the wave height increases approximately by a factor of 3 (gauge 1). The coastal gauges also demonstrate the formation of the wave reflected from the coast (gauges 2 and 4), propagation of the wave in the direction of increasing depths, and its corresponding transformation during multiple interactions with other components of the wave field and break lines of bottom topography (gauges 5 and 7).
[77] Deep water gauges do not record any difference in the characteristics of the wave field caused by different computational algorithms, thus, an increased order of the Adams scheme does not result in any advantages. No differences are observed in the results obtained using the Adams scheme within the nonlinear dispersion models and classical model of shallow water. This fact makes possible considering in the further description the latter ones as reference models. However, we note that the corresponding algorithm requires large computational resources and more sophisticated realization than the algorithm of the MacCormac scheme.
0 1000 2000 3000 4000 5000 0 1000 2000 3000 4000 5000
Figure 11. Transformation of solitary wave in the “wash-tub” basin. The results were obtained using the Adams scheme (h0/2). Scheme of location of gauges is shown on Figure 2.
[78] The effects of different approaches to the control of accuracy of calculations of waves at the coastal point are presented in Figure 12. It is clearly seen that an increase in the accuracy of the computational algorithm allows us to obtain acceptable results without decreasing the size of the grid (thin solid line in fragment (b)). This effect is observed both for the incident and reflected waves. The MacCormac scheme (fragment (a)) responses stronger to roughening of the grid. It is easy to see that the calculations performed using these grids can be used only for approximate estimates of the most general characteristics, such as, for example, the amplitude of the leading wave. Even in this case, more accurate schemes are preferable.
[79] Pressure gauge records shown in Figure 13 demonstrate that only slight difference in the description of the final fragment of gauge record observed at point 6 results from the serious deviations in the description of processes occurring at near boundary points (1-3). The calculations using a rough grid demonstrate not only a large error in the amplitude of the leading wave, but also a clear simplification of the form of free surface oscillations, which later manifests itself in the gauge records calculated far from the coast (points 4 and 5).
[80] The appearance of non-physical oscillations caused by insufficient resolution of the calculation grid is illustrated by the graphs shown in Figure 14 and Figure 15. It is manifested best of all if the Adams scheme is used together with a rough grid.
[81] The results of calculation of wave regime at the fifth point located in the coastal zone at a depth of 10 meters are shown in Figure 16. As was shown above, at this point the leading wave is practically separated from the wave reflected from the coast, and a time moment exists that divides these components of the wave field. However, it is easy to see (a) that application of rough grids (2h0 and h0) does not allow us to reproduce this phenomenon regardless increased accuracy of the Adams scheme. The graphs shown in fragment (b) evidence that the calculations with step h0 using the Adams scheme and h0/2 using the MacCormac scheme are close. A transition to the most exact grid (h0/4) is shown in fragment (c). The leading wave is described ideally by each of the algorithms, while the reflected wave is reproduced best by the most accurate algorithm at the best grid. The closest version to this one is the result obtained by a more exact algorithm over a less detailed grid. Finally, for comparison, the most “exact”, least “exact”, and “intermediate” gauge records are shown in fragment (d). Their comparison allows us to estimate the sense of increasing the computational resources needed to increase the accuracy in the description of quantitative and qualitative characteristics of the wave regime. Similar graphs for the first and fourth gauge records are shown in Figure 17. These graphs indicate that the MacCormac scheme, which is easy in realization, requires doubling of calculation points (compared to the Adams scheme) in order to provide a result acceptable by accuracy. We shall consider that the acceptable reference result is the curve calculated using the Adams scheme over a grid with step h0/2. Application of more exact scheme over a smaller size grid is necessary to reproduce some “fine” effects.
Figure 12. Dependence of the characteristics of the wave regime at point five on the accuracy of the algorithm near the boundary and on the size of the grid. MacCormac scheme of the second order of accuracy (a), Adams scheme with the second order of approximation at boundary nodes (b), Adams scheme with the fourth order of approximation at boundary nodes (c). Solid heavy curve: h0/4, dashed heavy curve: h0/2; solid thin curve: h0, dashed line: 2h0.
Figure 13. Pressure gauge records calculated using the MacCormac scheme with different steps (dashed line: h0; thin solid line: h0/2) of the second order of accuracy and Adams scheme of the fourth order with steps (solid line with marker: h0). Numerals in the graph field correspond to the numbers of gauges.
Figure 14. Point 8 (under the crest); the final part of gauge record is shown in the inset. Calculation using the scheme of Adams. Solid line is step of the grid h0/2, dashed line refers to 2h0.
6. Comparison Analysis of Two Programming Systems
[82] A comparison analysis of two programming systems (programming code packages) is an absolutely needed stage preceding the formulation of numerical modeling program. Here, we shall consider a specific example of such comparison and suggest a set of test problems for discussion, each of which is designed for determining and estimating specific characteristics of different models, algorithms, and program realizations.
[83] The main attention would be focused on the Nereus programming system developed by the authors and designed for numerical modeling of tsunami propagation from the source to the coast including wave interaction with island systems, coastal constructions, and flooding of the adjacent land [Eletsky et al., 2005].
[84] The capabilities of this system are compared with the capabilities of the program developed by specialists from Nizhniy Novgorod and Turkey on the basis of the TUNAMI code [Goto et al., 1997; Zaitsev et al., 2005] in its turn developed by Japan specialists in the 1980s. The recent modifications of this code claim to be a standard, thus they are widely spread in the community of the specialists involved in the solution of applied tsunami problems. Both systems are based on classical equations of the theory of shallow water. The algorithms realized by the systems are related to the class of explicit finite difference schemes. The Nereus and TUNAMI programs to one or another degree of flexibility allow us to perform calculations with or without account for nonlinear effects, rotation of the Earth, bottom and wind friction.
[85] The modeling results of Sumatra (2004) tsunami in the Indian Ocean became the material for comparison. The data obtained using one of the latest versions of TUNAMI system (we shall call it TUNAMI-M) were kindly given to us by E. N. Pelinovsky and A. I. Zaitsev, the active participants of the modification process and development of this code
[86] A number of tests were suggested by the authors to estimate the efficiency of programming systems. In all problems considered below, condition of zero transport (reflec-
tion) is specified along the coastline, which corresponds to the assumption about the existence of a vertical wall at the “land-sea” boundary. Thus, runup of waves is not considered in these problems.
[87] As was mentioned above, one of the main problems of numerical modeling of tsunami is reproducing the coastal wave regimes. In this relation, the problem about the quality of reproducing the boundary conditions at the coastlines, whose configuration complexity is close to the real ones, becomes especially important. This problem becomes the first, since the methods, which were used as the basis of many specialized (tsunami oriented) computational algorithms including the algorithms of the Nereus system, use uniform rectangular grids. Thus, the first test problems are directly related to the verification of the efficiency of numerical methods in the description of the interaction between waves and coastal structures of arbitrary form.
[88] Three of the problems considered here are fully model problems. The first problem (problem 1) describes the interaction of a solitary wave with the wall located normal to the direction of wave propagation. Schemes of the initial wave location relative to the coastal wall and location of gauges are shown in Figure 18.
[89] Here and above, we suggest using mainly gauge records to estimate the solution and perform comparative analysis. The modeling area is specified by a rectangle with sides
3267.000 m (in the direction of the wave propagation) and
2913.000 m. The corresponding computational grid has a size of 1090 by 972 points. The depth of the basin is constant equal to 1000 m. The vertical wall is located in the center of the modeling area. The front of the solitary wave with amplitude of 1 m is parallel to the wall. At the initial
Figure 15. Point 6: calculation using the scheme of
MacCormac (MC) (a); calculation using the scheme of Adams (A) with different resolution (b).
Figure 17. Points 1-4: calculation using the schemes of MacCormac and Adams with different resolution. Numerals Figure 16. Point 5: calculation using the schemes of in the graph field indicate the numbers of the points of
MacCormac and Adams with different resolution.
gauges.
Figure 18. Test problem 1: location of the initial perturbation of the free surface relative to the coast wall (a); scheme of location of gauges (b).
time moment, the wave crest is located over a point with coordinate 2451,000 m. Initial velocity field is specified so that the wave propagates without changing its form in the direction to the wall. Reflection conditions are specified at the lateral boundaries of the calculation area along the wave propagation, while at the back boundary free pass condition is specified. The physical time of the process is 24,000 seconds.
[90] Inclusion of a significant fragment consisting of land points, in which the calculations are not performed, is explained by the idea to conserve unique geometrical characteristics of the test problems.
[91] It is assumed in problem 2 that the coastal wall is located at an angle of 45° to the initial location of the wave front (see Figure 19), and the distance from the closest point of the wall to the front is 827,000 m.
[92] A classical problem of the interaction between a solitary wave with an amplitude of 1 m and a cylindrical island is considered as problem 3 (see Figure 20). The size of the calculation area and parameters of the grid are the same as in the previous problems. The cylindrical island has a radius of 110,000 m. Its center is located at a point with coordinates (816,750; 1456,500). Initial velocity field was specified zero. The other conditions were not changed. It is note-
Figure 20. Test problem 3: location of the initial perturbation of the free surface (a); scheme of location of gauges (b).
Figure 21. Test problem 4: initial perturbation of the free surface (a); bottom topography of the basin (b); scheme of location of gauges (c).
Marigrams
32-33-34-35-36
Figure 22. Test problem 5: bottom topography of the basin.
worthy that besides the analysis of the interaction with the coastal structures, the solution of the first three problems allows us to estimate the variations in the main characteristics of the solitary wave before, after, and during the interaction with the obstacles.
[93] The next series of the test problems is directed to distinguish the efficiency and principal possibility of using the test program systems for numerical modeling of tsunami. These are the Nereus and TUNAMI-M programs.
[94] Problems 4-6 are related to a certain degree to the modeling of the catastrophic Sumatra tsunami (26 December 2004).
[95] Problem 4 implies modeling of tsunami wave transformation in the Indian Ocean (see Figure 21). The initial perturbation of the free surface determined using the TOPICS program is located according to the generally accepted concepts about the source of the real event (Figure 21a). The calculation area (Figure 21b) extends from 70-110°E and from 10-25°S. The steps over spatial derivatives are chosen so that the computational grid contains 1090 and 972 nodes, respectively. The boundary conditions at external boundaries provide a free pass of the wave.
[96] Problem 5 is a simplification of problem 4: here, the real depths are changed to a constant value 1000 m (see Figure 22). The other parameters including the scheme of location of gauges (Figure 21c), the form and the location of initial perturbation remain unchanged. This problem is designed to distinguish the effects in the results of numerical modeling, which are caused exclusively by the interaction between the wave field and coastal boundaries. Such analysis is strongly simplified by introducing constant value of depth. It follows from the pattern in Figure 22 that even after simplification the problem of adequate reproducing of costal wave regime remains very complicated, which is caused by strong indentation of the coastal boundary and large number of small scale islands. This allows us to state
that after we obtained quantitatively and qualitatively close results using different programming systems, these systems can be recommended for the solution of applied problems if they passed the model tests (for example, problems 1-3 formulated above).
[97] The last test problem (problem 6) is the limiting simplification of problem 4. It is designed to analyze the influence of the form of the real initial displacement of the free surface on the components of the calculated wave field. The initial displacement of the free surface inherited from the two previous problems is located in a rectangular basin of permanent depth 1000 m. The geometrical parameters of the basin and its discrete image remain unchanged. As a rule, this test problem should be addressed to determine the possible causes of discrepancy between the results presented by different programming systems. The initial state of the wave field shown in Figure 23 is characterized by large gradients and poor smoothness of the carrier.
[98] The next series of figures presents the most general characteristics of wave regimes corresponding to the test problems considered above. No doubt, the results presented here are of purely illustrative character. A number of steps organizing testing and certification of programming systems designed for numerical modeling of tsunami should be made if serious work is planned. The materials should be prepared to use in the test programs including electronic presentation of the data in agreed formats. These materials would contain all the necessary input data and full results including the wave and velocity fields calculated at given time moments and gauge records at specified points.
[99] The series of images in Figure 24 demonstrates the transformation of a solitary wave during its interaction with a slant wall (problem 2). It is seen in these figures recorded with equal time intervals (the sequence is from left to right and from top to bottom) how reflected wave appears after the contact between the wave and the wall. The wave reflects according to the laws of geometrical optics. The angle
<
op
Figure 23. Test problem 6: initial perturbation of free surface and coastline contours of the nearest islands.
Figure 24. Results of numerical modeling of test problem 2 using the Nereus software system.
of reflection coincides with the incidence angle. The generated configuration is stable. It is conserved in general even after the reflection from the lateral wall, which results in the change of direction of motion of the entire wave construction. The formation of this V-shaped form and its motion without change of the form over a large distance with multiple reflections from lateral walls allow us to estimate the quality of realization of reflecting conditions both at the slant wall and at the lateral boundaries.
[100] The images of the free surface illustrated in Figure 25 demonstrate the well known stages of solitary wave evolution when at the initial time moment it decomposes into two waves (at zero initial velocities). One of them propagates through the open boundary. Here, we get a possibility to estimate the quality of the realization of the corresponding boundary condition. The other wave propagates in the direction to the cylindrical island. The amplitudes of both waves are half of the amplitude of the initial perturbation. The estimate of the wave height (amplitude) in the immediate vicinity of the island allows us to estimate the value of the scheme dissipation.
[101] As the wave approaches the island, it begins to interact with it, turns around the island and forms a system of concentrated reflected waves, which propagate in the opposite direction. The front of the wave after insignificant changes restores its characteristics and continues the motion in the initial direction, while a configuration is formed behind the front frequently called a “dove tail”. In addition to these qualitative results, important material for estimating the computational tools can be obtained from the analysis of the distribution of the maximal and minimal wave heights along the perimeter of the island. These results also facilitate the estimate of the quality of reflecting boundary conditions and determination whether the resolution of the computational grid is sufficient.
[102] The images shown in Figure 26 illustrate the most general concepts about the process determined by the conditions of the fourth test problem. These wave fields were calculated for the bottom topography closest to the real in the basin. Wave patterns clearly demonstrate the formation of the wave front directed to the south and southwest, its interaction with Sri Lanka Island and further propagation beyond
o o o <5
0 m
Figure 25. Results of numerical modeling of test problem 3 using the Nereus software system.
1 V 4 1g^j
11 [I > x
Figure 26. Results of numerical modeling of test problem 4 using the Nereus software system. A sequence of images of wave fields in the condition of Sumatra (2004) tsunami.
the calculation area in the direction to the African coasts. Part of the wave energy remains trapped in the straits and small islands in the vicinity of Sumatra Island. In general, the reproduced picture agrees qualitatively with the available information about the catastrophic event. However, we are not speaking about adequate correspondence to the real natural phenomenon. Most likely, the natural phenomenon
serves as a test setup to test the choice of reliable computational technologies.
[103] Comparative results of modeling in the solution of problems 4-6 using the Nereus (solid line) and TUNAMI-M (dashed line) are shown in Figures 27-31.
[104] As to the comparison of gauge records calculated using two different programming systems, it is reasonable to
20000 12000 Figure 27. Records of gauges 10 and 15 calculated for the 6th problem.
1 | ^ | r 17000 18000 19000 20000
Figure 28. Records of gauges 13 and 14 calculated for the 5th problem.
Figure 29. Records of gauges 20 and 25 calculated for the 4th problem.
Figure 30. Records of gauge 25 calculated for the 4th (a), 5th (b), and 6th (c) problems.
0.2
0.1 —
0 —
-0.1 —
0.2 —| _
I; 0.1 — 0.04 —
1 1 1 1 f ■ *
k, . i, A i i '■ Vi A ' r<\i v V ‘V.'M i1 \i '< * 0 — 1 life 'L J A ,11 i 0 —
(a) -0.1 — V ft t 1 '! * 1 « (b) -0.04 —
i 1 i 1 i 1 i 1 1 1 1 1 1
(C)
8000
10000 12000
14000
8000
12000
16000
20000 8000
12000
16000
n I 20000
Figure 31. Records of gauge 22 calculated for the 4th (a), 5th (b), and 6th (c) problems.
start the comparison from the simplest problem 6 (in the absence of variations in the depths and coastal boundaries). On the basis of the analysis of gauge records shown in Figure 27 we can speak not only about the qualitative but also about quantitative coincidence of the modeling results. We note that different conditions at the boundaries of the computational area are used in the calculations presented here. This is clearly seen in the record of gauge 15. Natural condition of free pass of the wave was imposed at the southern boundary in the calculations using the Nereus program, while in the TUNAMI-M programming system another boundary condition was specified.
[105] The comparison of the results of the solution of more complex problem 5 (see Figure 28), in which the real form of the coastal boundaries is conserved, while the depth remains constant, demonstrates that the deviation of the results becomes more notable, however, the leading waves are reproduced almost identically. In some cases (not included in the material described here) we can speak about complete qualitative coincidence of the results. The cause of the noted differences can be a risky joint use of linear and nonlinear models of shallow water in the algorithm of the TUNAMI-M program, which requires for the correct realization a thorough sewing the solutions with account for the difference in the directions of characteristics used for transition of solutions in hyperbolic equations. The Nereus programming system uses only nonlinear theory of shallow water.
[106] Finally, the results of modeling of problem 4 (see Figure 29), which is the closest to the real tsunami phenomenon, indicate that regardless that the leading waves are described almost identically, the results diverge stronger and stronger in the course of time.
[107] The analysis of the entire set of the materials of calculations (no doubt, there are points, in which solutions differ not so strongly) leads to a conclusion that the cause of these differences is most likely the different approaches of the authors of programming systems to testing of algorithms both at internal and boundary points. It is our opinion that distinguishing of such differences should not be limited only by stating of this fact. We should apply all forces to understand their cause and eliminate the error in algorithms and programs if there are any. Test problems 5 and 6 were suggested with this goal in mind. The analysis of their solutions
allowed us to determine the possible means of solving this crisis. We note that such work requires certain enthusiasm and personal will for the cooperation between the developers of the models, algorithms, and programs, as well as organization forcing on behalf of the entire community of the specialists.
[108] Concluding the description of the results of numerical solution of problems 4-6 we shall briefly describe the above mentioned analysis of the influence of real bottom topography, indentation of the coastal boundaries, and combination of these factors on the solution. Pressure gauges presented in Figure 30 and Figure 31 are the illustrations of the perspectives of such analysis using the materials of calculations. The presence of high frequency oscillations in the results obtained using the TUNAMI-M programming system is worth attention. The appearance of such oscillations in problem 6 with the bottom of constant depth and absence of internal boundaries points to the purely computational nature of such oscillations and impossibility of associating them with any physical sense.
[109] The further interaction of the “numerical” oscillations with irregularities of bottom topography and peculiarities of the coastal boundaries can become a cause of increasing distortion of the solution of more complex problems 4 and 5 and increasing discrepancy with the results of calculations using another programming system. In such cases, one should seriously treat the reveake problem and make efforts to overcome it.
7. Conclusion
[110] At present, many programming systems of different level of organization for numerical modeling of tsunami are available for the scientific community investigating the tsunami wave dynamics. The systems differ from each other in the variety of applied mathematical models and numerical algorithms as well as in the choice of the programming medium depending on the formulated problem, purpose of the programming product, and qualification of the possible user. Full review and comparative estimates over the entire variety of programming systems for tsunami modeling and
algorithms realized for these purposes are extremely difficult owing to a number of reasons. This study is one of the attempts to make a comparative estimate of the efficiency of the known computational models and programming systems of hydrodynamic tsunami modeling based on these models.
[111] Acknowledgments. This work was supported by the Council of grants of the President of the Russian Federation for supporting leading scientific schools of the Russian Federation (grant SS-9886.2006.9), Russian Foundation for Basic Research (grants 05-05-64460, 06-05-64869), Program of Multi-disciplinary integrated research of the Siberian Division of the RAS (project 2006-2.12), and Program of Interdisciplinary Integration research of the Siberian Division of the RAS (projects 2006-28, 2006-113).
References
Eletsky, S. V., Z. I. Fedotova, and L. B. Chubarov(2005), Computer model of tsunami waves, in Proc. of Tenth Baikal All-Russia Conf. “Information and Mathematical Technologies in Science, Technology, and Education”, Part 1, p. 138, ISEM SO RAN, Irkutsk.
Fedotova, Z. I. (2006), Application of MacCormac differ-
ential scheme for the problems of ling wave hydrodynamics, Special issue dedicated to 85 Anniversary of N. N. Yanenko, Computational Technologies, 11(Part 2), 53.
Goto, C., Y. Ogava, and F. Imamura (1997), Numerical method of tsunami simulation with the leap-frog scheme (English translation and preparation by N. Shuto), in IUGG/IOC Time
Project, IOC Manual and Guides, no. 35, p. 126, UNESCO, Paris.
Gusyakov, V. K., and L. B. Chubarov (1985), Tsunamis and earthquake mechanisms in the island arc regions, Sci. Tsunami Hazards, 3(1), 3.
Gusyakov, V. K., and L. B. Chubarov (1987), Numerical modeling of generation and propagation of tsunami in coastal zone, Izv. Phys. Solid Earth, 23(11), 53.
Lynett, P. J., and P. L.-F. Liu (2004), A two-layer approach to water wave modeling, Proc. R. Soc. London, Ser. A, 460, 2637.
Nwogu, O. (1993), Alternative form of Boussinesq equa-
tions for near shore wave propagation, J. Waterway, Port, Coastal and Ocean Engng., 119, 618, doi:10.1061/(ASCE)0733-950X(1993)119:6(618).
Shokin, Yu. I., and L. B. Chubarov (1999), Modern information technologies is an effective toll for solving crisis situations, ECO, (12), 15.
Stoker, J. (1959), Waves in Water, 618 pp., Inostrannaya
Literatura, Moscow.
Wei, G., and J. T. Kirby (1995), A time-dependent numerical code for extended Boussinesq equations, J. Waterway, Port, Coastal and Ocean Engng., 120, 251, doi:10.1061/(ASCE)0733-950X(1995)121:5(251).
Zaitsev, A. I., A. A. Kurkin, B. V. Levin, E. H. Pelinovsky, A. Yalchiner, Yu. I. Troitskaya, and S. A. Ermakov (2005), Modeling of the propagation of catastrophic tsunami (26 December 2004) in the Indian Ocean, Dokl. Ros. Akad. Nauk, 402(3), 388.
S. A. Beizel, L. B. Chubarov, S. V. Eletsky, Z. I. Fedotova, Yu. I. Shokin, Institute of Computational Technologies, Siberian Branch of the RAS, Novosibirsk, Russia
