Extreme Inundation Flows
During The Hokkaido-Nansei-Oki Tsunami

Vasily V. Titov
PMEL, NOAA, Seattle,Washington

Costas Emmanuel Synolakis
University of Southern California, Los Angeles

Abstract. The tsunami generated by the July 12, 1993 Hokkaido-Nansei-Oki Mw = 7.8 earthquake produced in Japan the worst local tsunami-related death toll in fifty years [Hokkaido Tsunami Survey Group, 1993; Shuto and Matsutomi, 1995], with estimated 10-18m/sec overland flow velocities and 30m-high wave runup. These extreme values are the largest recorded in Japan this century [Hokkaido Tsunami Survey Group, 1993; Shuto and Matsutomi, 1995] and are among the highest ever documented for non-landslide generated tsunamis. We model this event to confirm the estimated overland flow velocities, and we find that current state-of-the-art shallow-water wave models can predict tsunami inundation correctly including extreme runup, current velocities and overland flow. Our results qualitatively suggest that for this event coastal inundation is more correlated with inundation velocities than with inundation heights, thereby providing one explanation why threshold-type modeling has substantially underpredicted coastal inundation in this and other recent events.


Important as it is for hazard mitigation planning for great or giant tsunamis such as the one attributed to the Mw ~ 9 Cascadia Subduction Zone event[Geist and Yoshioka, 1996, Satake et al., 1996] of 1700, the calculation of three-dimensional tsunami inundation flows remains formidable. Up until recently [Shuto, 1991, Liu et al., 1995, Myers and Baptista, 1995, Synolakis, 1995 Takahashi et al., 1995&1996], field-applicable inundation predictions were only possible with what is referred to as threshold-type models; the latter only propagated the tsunami wave from generation up until the 10m depth contour, and then used the wave height at that depth to infer the inundation. The field surveys of the 1992-1995 tsunamis showed that threshold models sometimes produced predictions differing by a factor ranging from five to ten compared with field observations, raising concern because of their use for certain civil defense applications, and fueling controversy as to whether the poor predictions were due to insufficient ground deformation data affecting the initialization of the calculations, or to inherent limitations in the propagation model.

Recent advances in hydrodynamics through validation of model results with large-scale laboratory experiments [Titov and Synolakis, 1993, Yeh et al., 1994 Liu et al., 1995, Briggs et al., 1995] have demonstrated excellent predictive capabilities of the state-of-the-art inundation algorithms, yet the extrapolation of these results to prototype tsunamis is not unequivocal. In the physical models, idealized and well-defined solitary waves evolve over constant depth and then over idealized bathymetry, while in nature the initial wave profile is never known in detail, and the wave propagates over complex bathymetry and climbs up topography known only to finite resolution ; the latter was not believed to be an important constraint, given that typical tsunami wave scales are larger than microtopographic features, until recent field data demonstrated otherwise [Yeh et al., 1993; Synolakis et al. 1995].

The Mw = 7.8 Hokkaido-Nansei-Oki earthquake of July 12, 1993 was a fortuitous large-scale experiment which allowed the measurement of high-quality runup data and the inference of fairly unambiguous ground deformation contours due to its proximity to seismic instrument arrays in a locale where high-resolution bathymetry and topography data exist. Hydrodynamic computations reported to date [Myers and Baptista, 1995, Satake and Tanioka, 1995, Takahashi et al., 1995] using the shallow-water wave approximation (SW) have reproduced to first-order the runup heights distribution along the coast of Hokkaido and to some extent around Okushiri, but neither the wave current velocities nor the extreme runup values observed in the front of the island; at best [ Takahashi et al., 1995], predictions for the extreme measurements have differed by a factor of two from the field data. Further, the region of most extensive devastation did not correlate well with the highest measured runup. Given that the Hokkaido-Nansei-Oki event was not a tsunami earthquake [F. Schindelé et al, 1995], i.e., an earthquake which produces a tsunami disproportionally large for its size, these differences are intriguing, and they have implications in planning for large tsunami disasters.

To explore the relative effects and magnitudes of of inundation heights and inundation-flow velocities, and to determine whether the modeling difficulties are due to inherent limitations of the SW models or their application outside their range of validity or to splash, or to resolution limits, or even whether to the notoriously difficult [Aida, 1977, Iwasaki and Mano, 1979, Hibberd and Peregrine, 1979, Shuto, 1991, Liu et al., 1991] and unstable [Synolakis, 1989] shoreline computations, we investigated this tsunami using a new inundation model.

The mathematical model

We solve the 2+1 nonlinear shallow water wave equations,

ht + (uh)x + (vh)y = 0,
ut+ uux + vuy + ghx = gdx,       (1)
vt+uvx + vvy + ghy = gdy.

Here h = h(x,y,t) + d(x,y,t), where h(x,y,t) is the amplitude, and z = -d(x,y,t) is the location of the bottom with respect to the undisturbed water surface, while u(x,y,t), v(x,y,t) are the depth-averaged velocities in the onshore x and the longshore y directions respectively; g is the acceleration of gravity. These equations are believed to model the dynamics of the runup process well. To avoid confusion, we will refer to these equations as the 2+1 NSW equations to indicate that there are two spatial propagation dimensions and one temporal one; the 1+1 NSW equations only model wave evolution in the x direction. We solve without adding dissipation terms. Although physically bottom friction does affect the wave evolution, it has been shown [Kobayashi et al., 1987, Zelt, 1991, Kobayashi, and Wurjanto, 1992, Liu et al., 1995, Raubenheimer and Guza, 1996] that inundation results are not sensitive to changes in the roughness coefficient, at least when these results are compared with laboratory data. Also, since the bottom-friction coefficient can only be inferred from related studies or by model self-calibration, we do not use it, particularly because we have found our method sufficiently robust.

We will use explicit finite-differencing and we will split the system of three 2+1 coupled equations (1) into two 1+1 quasi-linear systems. Note, that the conventional practice is to use an implicit scheme with the splitting method for solving the resulting one-dimensional equations; this approach works well for elliptic and parabolic equations, reducing the number of operations substantially compared with applying an implicit scheme applied directly on the 2+1 equations [Fletcher, 1991]. Here, (1) is a hyperbolic quasi-linear system, yet we did find it advantageous to use the splitting method in combination with an explicit FD technique. The main advantage is the use of the characteristic form of the 1+1 NSW equations, which allows the definition of of a well-posed boundary value problem for the numerical method. The characteristic form of equations also allows for an efficient FD realization [Titov and Synolakis, 1995], using the fractional step method of Yanenko [1971].

The most critical step in the inundation calculation is the shoreline computation, which involves a moving boundary interface between three different phases, liquid air and solid. The 1+1 shoreline computation was a classic problem of hydrodynamics and it is discussed in detail by Shuto [1991]. 2+1 inundation calculations have been reported recently [Shuto, 1991, Myers and Baptista, 1995, Takahashi et al., 1995 , Liu et al., 1995], but all include friction factors.

In our inundation computations, it is impossible to use the direct relationships between the Riemann invariants near the shoreline, since here the Froude number may be greater than one, and all characteristic families have the same inclination in this region. Therefore, we use approximations of the boundary values from previous space nodes. Our shoreline algorithm uses the time-dependent space step Dx(t) of the last node of the computational area. The objective is maintain the shoreline boundary point on the surface of the beach during the computation, as first suggested by Iwasaki and Mano [1979]. By adjusting the length of the last space-step every time-step, we ensure so that the shoreline point is at the intersection of the beach with the horizontal projection of the last "wet" point. The value of the velocity on the shoreline node is equal to the velocity of the previous "wet" point. We will refer to our method henceforth as VTCS-3; it has been shown in excellent agreement with data from large-scale laboratory experiments on wave runup on circular islands [Titov and Synolakis, 1995], even for mildly plunging waves, where the energy dissipation of the physical manifestation is modeled through our particular FD scheme's inherent dissipation.

Since our objective is to primarily evaluate the predictive power of the hydrodynamic model under extreme conditions and not the relative merits of different source models, we will only use the three-fault plane dipole-shaped ground deformation model of Takahashi et al. [1995], known as DCRC-17a, as our initial condition. DCRC-17 predicts a maximum uplift of 4.2m and subsidence of 1.1m, and hydrodynamic calculations using it have produced the closest qualitative agreement with the field observations to date among different source models [Myers and Baptista, 1995, Takahashi et al., 1995, Myers and Baptista, 1995, Satake and Tanioka, 1995].

Finally, in our simulations we used the digitized bathymetry as provided by the Disaster Control Research Center in Tohoku University consisting of several nested bathymetric arrays with three level of grid resolution, ranging from 450m offshore to 50m in several locations close to the initial shoreline.

Figure 1 pdf

Figure 1. Comparison of field measurements (empty circles) with VTCS-3 inundation predictions (solid circles).


Figure 1 shows the distribution of computed runup heights around Okushiri using VTCS-3 and it compares them with the field measurements [Hokkaido Tsunami Survey Group, 1993, Shuto and Matsutomi, 1995]. The computed distribution of the maximum runup heights reproduces most features of the field measurements well. Note the variation in the computed runup height over small longshore distances, verifying that the variation in the field data over the same areas is due to hydrodynamic effects and not to random observation scatter.

The extreme 31.7m measurement near Monai was computed as 29.7m. The measurement was made at the tip of a small canyon and it is undoubtedly a local effect; yet all measurements around it were consistently high with a global maximum at 31.7m, similar to the distribution of computed values shown in the figure; this measurement distribution was not present in the other simulations [Myers and Baptista, 1995, Satake and Tanioka, 1995]. Takahashi et al [1995] produced qualitatively and quantitatively correct distributions, but only predicted maximum runup of 15m in this locale; we speculate that their friction factor disproportionally affect the larger values. Most of our computed runup values are a slightly higher than measured, but then our model does not include any dissipation. The largest difference is near Inaho at the north, where the calculations predict 15m heights compared to the 10m field measurements [Hokkaido Tsunami Survey Group, 1993] ; interestingly, 12-15m field mearurements have been reported at Inaho by other japanese reconnaissance teams [Nakasuji and Takahashi, 1993]. Overall, for this extreme tsunami event the SW model appears adequate for modeling runup heights, suggesting that wave breaking and bottom friction are apparently not important to first-order for predicting inundation.

Figure 2 pdf

Figure 2. Comparison of obsrvations (stars) with inundation model using different grid resolutions and depth-threshold model predictions.

We then performed a series of numerical experiments to evaluate the effect of other aspects of numerical modeling, such as the grid resolution and the 10m-depth threshold. Figure 2 compares runup results among computations with different grid resolutions for a region of the west coast of Okushiri island where topographic data exist down to 50m-square resolution. For reference, the leading N-wave close to the shoreline [Tadepalli and Synolakis, 1995&1996] has a O (1km) wavelength , and therefore the computations shown range from 3-4 points per wavelength for the 450m grid to 20 points per wavelength for the 50m grid. This number may appear small by computational fluid dynamics standards (CFD) standards, but then as Shuto [1991] pointed out, the effects of the grid resolution are often counterintuitive, because of numerical dispersion.

The computations with the 150m grid produce qualitatively correct runup distribution, but miss the extreme runup height. The computations with the 450m grid, even though adequate in terms of hydrodynamic stability to resolve the waveform [ Titov and Synolakis, 1995], these calculations seriously underpredict the runup height, suggesting that it is the effect of the small-scale features hidden in the coarse grid and not only the effect of wavelength resolution. This is rather important, as problems arising only from wavelength resolution can be solved easily by generating denser grids (through interpolation of the coarser grid), yet these denser grids have identical topography resolutions as the coarse grids that produced them.

In figure 2, we also present results from numerical experiments we did to compare the predictions of models with 10m-contour calculation thresholds. These models are still ubiquitous [Yeh et al., 1993, Synolakis, et al. 1995, Satake and Tanioka, 1995]; and they propagate the wave using a NSW model up to the 10m contour, and then use the maximum wave height at that location to infer the maximum runup height. For identical grid resolutions, our results show that interrupting the computation at the 10m depth, a step essentially equivalent to placing reflective "wall"-type boundaries at that depth, underpredicts the runup by a factor of two, as compared with the inundation computations. Even the computed maxima at the 10m contour from the full inundation calculation differ substantially from the computed maxima from the 10m-threshold "wall" calculations, suggesting that the practice of using the waveheight at the 10m depth to infer the runup to first order needs to be re-evaluated. Clearly the wave evolves substantially as it propagates from the 10m depth up the beach to its maximum runup.

Figure 3 pdf

Figure 3. Snapshots of the computed overland flow in Aonae, at 5 min (a), 6 min (b), 8 min (c) and 12 min (d) after the earthquake.

Another challenging computation for the SW theory is overland flow, a condition referring to waves overtopping a peninsula or a narrow strip of land sandwiched between the ocean and a tidal inlet, then propagating over it and then draining on the lee side of the land. Overland flow was first described by Yeh et al. [1993], and later identified [Synolakis et al., 1995] as one of the leading causes of the heavy casualties during tsunami attack because of the high flow velocities; during tsunami runup on a sloping beach, the shoreline fronts slows down as its kinetic energy is converted into potential energy, but during overland flow over flat land the kinetic energy is only reduced by dissipation. Overland flow is usually supercritical with bore-like dynamics. We present overland flow results near Aonae in the four flow snapshots of figure 3, before, during and after the overland flow. The dynamics are very similar to those inferred from the field observations [ Hokkaido Tsunami Survey Group, 1993, Shimamato, 1995]. The computed flow velocities range from 10-19 m/s over the Aonae peninsula, consistent with the 10-18 m/sec field estimates of in Aonae [Shimamato, 1995]. These values are consistently higher compared with elsewhere around the cape and they correlate well with the devastated area.

Figure 4 pdf

Figure 4. Computed envelope of maximum tsunami heights and of maximum wave velocities over one typical topographical cross-secsion of Aonae cape.

Figure 4 presents the envelope of maximum tsunami heights and of maximum wave velocities over one typical topographic cross-section of Aonae cape. Notice that the highest inundation velocities occur at the east side of Aonae where inundation heights are smallest, and how these extreme velocities correlate at least qualitatively with the most devastated area. This is entirely consistent with the fact that the during the lee-side rundown of the overland flow, the flow depths are small, but the flow is supercritical with large flow velocities, and it underlines the risk of relying exclusively on runup heights to predict inundation.


Despite the inherent difficulties of comparing tsunami field data with model predictions, we conclude that the SW approximation and current state-of-the-art solution methods are capable of modeling quantitatively correctly the tsunami inundation, including extreme runup heights and inundation velocities. We find that even small local bathymetric structures affect the runup to first-order and that the resolution of the bathymetric data maybe more important than the grid resolution. Our results also suggest-albeit only qualitatively- that during overland flow, coastal devastation correlates stronger with inundation velocities than inundation heights. This difference along with the often substantial wave evolution from the 10m depth to the maximum runup suggest caution in interpreting threshold-depth type calculations as they may miss extreme events.

Acknowledgements. We thank Mamoru Takahashi from Kokusai Kogyou Co. for providing the aerial photograph of the Aonae Cape and Dr. Tomu Takahashi of DCRC of Tohuku University for the digitized bathymetry. We are grateful to Dr. Cliff Astill of the Natural Hazard Mitigation Program of the National Science Foundation for supporting the development and validation of our research.


Aida, I., Numerical experiments for inudnation of tsunamis, Bull. Earth., res. Inst., 52, 441-460, 1977.
Bernard, E., Mader, C., Curtis, G., Satake, K., Tsunami inundation Model study of Eurica and Crescent city, California, NOAA Technical Memorandum ERL PMEL-103, 1994.
Briggs, M. J., Synolakis, C. E., Harkins, G. S., Green, D. R., Laboratory experiments of tsunami runup on a circular island, PAGEOPH 144 (3/4), 569-593, 1995.
Geist E., and Yoshioka, S., Natural Hazards 13, 151-177, 1996.
Hibbert, S., and Peregrine, D. H., Surf and runup on a beach: a uniform bore, J. Fluid Mech. 95, 323, 1979.
Hokkaido Tsunami Survey Group, Tsunami devastates Japanese coastal Regions, EOS, Trans. AGU 74 (37), 417-432, 1993.
Iwasaki, R., and Mano, A., Two-dimensional numerical simulation of tsunami runup in the Eulerian description, Proc. 26th Conf. Coastal Eng., JSCE, 70-74, 1979.
Kobayashi, N., Otta, A.K. and Roy, I., Wave reflection and runup on rough slopes, J. Waterways, Port, Coastal and Ocean Eng. 113 (3), 282, 1987.
Kobayashi, N., and Wurjanto, A., Irregular wave setup and run-up on beaches, J. Waterways, Harbors, Port, Coastal and Ocean Eng. 118 (4), 368, 1992.
Liu, P., L-F., Synolakis, C.E., Yeh, H., Report on international workshop on long wave runup, J. Fluid Mech. 229, 675, 1991.
Liu, P., L-F., Cho, Y.-S., Briggs, M. J., Kanoglu, U., Synolakis, C.E., Runup of solitary waves on a circular island, J. Fluid Mech. 302, 259-285, 1995.
Myers, E. P. and Baptista, A. M., Finite Element modeling of the July 12, 1993 Hokkaido-Nansei-Oki tsunami Pure and Appl. Geophys. 144, (3/4), 769-802, 1995.
Nakasuji T. and Takahashi, M. Field report of the July 12, 1993 Hokkaido earthquake, Kokusai Kogyou Ltd, Tokyo T102, Japan, 80pp.
Raubenheimer, B. and Guza, R. T., Coastal Engineering in press ,1996.
Satake, K. and Tanioka, Y., Tsunami generation of the 1993 Hokkaido-Nansei-Oki earthquake Pure and Appl. Geophys. 144, (3/4), 803-822, 1995.
Satake K., et al.,Nature, 379 (6562), 246, 1996.
Schindelé, F., Reymond, D., Gaucher E., and Okal, E.A., Analysis and automatic processing in near field of eight 1992-1994 tsunamigenic earthquakes: Improvements towards real-time tsunami warning, Pure and Appl. Geophys. 144 (3/4), 381-408, 1995.
Shimamato, T., Tsutsumi, A., Kawamoto, M., Miyawaki, M., and Sato, H., Field survey report on tsunami disasters caused by the 1993 Southwest Hokkaido earthquake, Pure and Appl. Geophys., 144 (3/4), 665-692, 1995.
Shuto, N., Numerical simulation of tsunamis, in Tsunami Hazard Bernard, E. (ed), Kluwer Academic Publishers, Doldrecht, The Netherlands, 171-191, 1991.
Shuto N. and Matsutomi H., Field survey of the 1993 Hokkaido Nansei-Oki tsunami, Pure and Appl. Geophys. 144 (3/4), 649-664, 1995.
Synolakis, C.E., Discussion, J. Waterways, Harbors, Port and Coastal Eng. 115 (1), 139, 1989.
Synolakis, C. E., Imamura, F., Tinti, S., Tsuji, Y., Matsutomi, H., Cooc, B., and Usman M., The East Java tsunami of the July 4 , 1994, EOS Transactions AGU 76 (26), 257,261-262, 1995.
Synolakis, C.E., Tsunami Prediction, Science 270, 15, 1995.
Tadepalli S., and Synolakis, C.E., The run up of N-waves on sloping beaches, Proc. Roy. Soc A, 445, 99, 1995 .
Tadepalli S., and Synolakis, C.E., Model for the leading waves of tsunamis Phys. Rev. Letters 77, (10), 2141-2144, 1996.
Takahashi, To., Takahashi, Ta., Shuto, N., Imamura, F. and Ortiz, M., Source models for the 1993 Hokkaido-Nansei-Oki earthquake tsunami Pure and Appl. Geophys. 144 (3/4), 747-768, 1995.
Titov, V. V. and Synolakis, C. E., A numerical study of wave runup of the September 1, 1992 Nicaraguan tsunami, Proc. Intern. Tsunami Symposium, 627, 1993.
Titov, V. V. and Synolakis, C. E., Modeling of Breaking and Nonbreaking Long Wave Evolution and Runup using VTCS-2, J. Harbors, Waterways, Port, Coastal and Ocean Eng. 121 (6), 308-316, 1995.
Yeh, H. H.,Imamura, F., Synolakis, C. E., Tsuji, Y., Liu, P. L.-F. and Shi, S., The Flores island tsunamis, EOS Trans. AGU, 74 (33), 369,371-373, 1993.
Yanenko, H.H., The Method of Fractional Steps, Springer, New York, NY, 1971.
Yeh,H., Liu, P.L.-F., Briggs, M.and Synolakis, C. E., Propagation and amplification of tsunamis at coastal boundaries, Nature 372, 353-355, 1994.
Zelt, J.A., The runup of nonbreaking and breaking solitary waves, Coastal Engineering, Coastal Engineering 15, 205-246, 1991.
Tsunami HOME