RUSSIAN JOURNAL OF EARTH SCIENCES VOL. 8, ES6004, doi:10.2205/2006ES000216, 2006

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 (reflection) 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.

2006ES000216-fig18
Figure 18
[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 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.

2006ES000216-fig19
Figure 19
[91]  It is assumed in problem 2 that the coastal wall is located at an angle of 45o 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.

2006ES000216-fig20
Figure 20
[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 noteworthy 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).

2006ES000216-fig21
Figure 21
[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-110o E and from 10-25o 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.

2006ES000216-fig22
Figure 22
[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).

2006ES000216-fig23
Figure 23
[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.

2006ES000216-fig24
Figure 24
[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 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.

2006ES000216-fig25
Figure 25
[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.

2006ES000216-fig26
Figure 26
[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 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.

2006ES000216-fig27
Figure 27
2006ES000216-fig28
Figure 28
2006ES000216-fig29
Figure 29
2006ES000216-fig30
Figure 30
2006ES000216-fig31
Figure 31
    
[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, 28, 29, 30, and 31.

[104]  As to the comparison of gauge records calculated using two different programming systems, it is reasonable to 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.


RJES

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.

Copyright 2006 by the Russian Journal of Earth Sciences

Powered by TeXWeb (Win32, v.2.0).