Tsunami Modeling

We have been developing finite volume wave propagation methods for the shallow water equations particularly for the purpose of tsunami modeling. In order to accurately model tsunami propagation and inundation, numerical methods for the shallow water equations must have certain properties that more traditional methods may not possess. See the description at the bottom of the page for more details. This research is conducted with Randall LeVeque, and is the numerical part of a larger collaborative effort with Marsha Berger at NYU, Harry Yeh at Oregon State University and others. This work is funded in part by the NSF and the DOE.

TsunamiClaw---an extension of the freely available Clawpack and AMRCLAW software---has been used to model tsunami propagation and inundation for hazard mitigation and scientific research. Validation and maintenance of this code is on an ongoing basis.

TsunamiClaw will soon be incorporated into a more general software package GeoClaw. This is a new and ongoing project. The project page is at www.geoclaw.org.

Simulation of the 2004 Indian Ocean Tsunami

The 2004 Indian Ocean tsunami was computed using TsunamiClaw---an extension of the Clawpack software. The tsunami was generated dynamically by moving the sea floor surface at the start of the computation using a model of the spatial-temporal pattern of the Sumatra-Andaman fault rupture. This model was provided by the Seismolab at CalTech. For a description of the numerical method, see the section at the bottom of this page.

Click below for movies of the tsunami simulations:

(These are pretty large---cache once then play again)

Madras India, Harbor Inundation (click on right image for simulation)

The Jason-1 satellite altimetry of the tsunami

The Jason-1 and Topex/Poseidon satellites, under a joint project between NASA and the French Space Agency, orbit the earth collecting satellite altimetry data for oceanographic research. The satellites happened to be passing over the Indian Ocean on December 26 2004. Researchers at the space agencies were surprised that the satellites apparently detected the tsunami. Below on the left is a plot of the sea-surface height data (ssha) from cycle 108 in blue, which was the previous time the satellites were on the same pass (about 10 days earlier), and cycle 109 in red which occurred during the tsunami. A crude estimate of the tsunami therefore can be obtained by subtracting the two data sets. (Earlier cycles and later cycles match cycle 108 fairly consistently.) The figure on the right shows the subtracted data sets (in red) and my computed solution (in blue).

1993 Okushiri Island, Japan

In 1993 a deadly tsunami struck Okushiri Island Japan. Near the harbor town of Monai, extreme runups of 32 meters were observed. This tsunami was used as a benchmark problem for the 3rd International Workshop on Long Wave Run-up Models. The goal of the problem was not to compute the actual event, but to compare to laboratory data from a wave tank scale model.

Click below for movies:

More of our results and benchmark problems from this workshop are available at www.amath.washington.edu/~rjl/catalina04

The Numerical Method

We are solving the nonlinear shallow water equations, which are the commonly accepted approximation for tsunami propagation. Additionally we solve the equations in the physically relevant conservative form, which presents some difficulties for deep ocean propagation, but may be important in the near shore zone if there are steep gradients or shocks in the solution. I am developing a finite-volume method for this application, based on the wave-propagation algorithms of Randall LeVeque. These are high-resolution shock-capturing methods developed more generally for hyperbolic conservation laws. They have the property of 2nd order convergence to smooth solutions (high-resolution), yet maintain the ability to converge to discontinuities (shock-capturing), such as shock waves in gas dynamics, or propagating bores in shallow water flow. A key component of this class of methods is the solution of a Riemann problem at each grid cell interface at each time-step. I have developed some approximate Riemann solvers, specifically designed to overcome some difficulties associated with this application. Some of these difficulties are explained below.

The Shallow Water Equations

The shallow water equations are a nonlinear set of hyperbolic conservation laws governing a free-surface flow, such as a body of water under the influence of gravity. These integral conservation laws, or the resulting PDEs, are often used to model flows where the length scale of the flow is large compared to the depth of the flow. This makes them the relevant approximation for tsunami modeling, since, although the ocean is not typically thought of as shallow, it is shallow compared to the long wavelength of tsunamis.

Aside from the usual numerical difficulties encountered with hyperbolic conservation laws (shocks, non-uniqueness, etc.) further difficulties are encountered in the presence of variable bottom bathymetry, which contributes a source term to the momentum equation. Although source terms are not universally a problem, when modeling mostly flat pools of water over variable topography, the stationary steady-state is maintained by a non-trivial balance between the source term and the momentum flux. A robust numerical method must perfectly preserve this steady-state, as well as accurately resolve small perturbations to it. This difficulty is exacerbated in global scale tsunami modeling, where ocean depths may reach 3km or more, which means that tsunami waves are small perturbations to ocean depth. Use of a linearized model does not circumvent this problem, since when approaching the shore the topography quickly changes, and the waves become large compared to the water depth. In these regions the nonlinear properties of the waves are relevant, as is the true conservation form of the equations which models shocks or discontinuities correctly. It is in this near-shore zone that accuracy in the solution is of most importance, and this region presents a new difficulty of its own: evolving dry regions and movement of the shoreline. It is this movement that demonstrates the inundation of the tsunami, and capturing the shoreline correctly is the ultimate goal of the numerical model.

Adaptive Mesh Refinement

Tsunamis belong to a class of geophysical problems with vastly diverse spatial scales of interest, presenting a difficulty to numerical modeling. For instance, a global scale tsunami may result in wave run-up and inundation that varies greatly even along local stretches of coastline. The desired accuracy at local scales needed for inundation mapping or early warning systems requires grids that are not feasible to use at the global scale. Since the propagating energy of a global scale tsunami is often focused by bathymetry in unpredictable ways, regions needing grid refinement are not known a priori, preventing the use of fixed telescoping grids, or similar grid-types, for real-time tsunami prediction. We have used adaptive mesh refinement algorithms to numerically model global scale tsunamis. These algorithms allow regions with grids of varying refinement where the solution has steep gradients or other features of interest. The regions of refinement may move adaptively with the solution over time, allowing various resolutions in a single global scale computation. These adaptive mesh algorithms were originally developed by Marsha Berger for gas-dynamics and other hyperbolic equations exhibiting shocks. They have been modified extensively for this particular application.

· Home ·