Tsunami science has evolved differently from research on other extreme natural hazards, primarily because of the unavailability until recently of instrumental recordings of tsunamis in the open ocean. Here, the progress towards developing tsunami inundation modelling tools for use in inundation forecasting is discussed historically from the perspective of hydrodynamics. The state-of-knowledge before the 26 December 2004 tsunami is described. Remaining aspects for future research are identified. One, validated inundation models need to be further developed through benchmark testing and instrumental tsunameter measurements and standards for operational codes need to be established. Two, a methodology is needed to better quantify short-duration impact forces on structures. Three, the mapping of vulnerable continental margins to identify unrecognized hazards must proceed expeditiously, along with palaeotsunami research to establish repeat intervals. Four, the development of better coupling between deforming seafloor motions and model initialization needs further refinement. Five, in an era of global citizenship, more comprehensive educational efforts on tsunami hazard mitigation are necessary worldwide.
The evaluation of the terminal effects of natural hazards remains one of the holy grails of geophysical research. Uncertainties arise in the specification of the initial condition, the calculation of the evolution of the geophysical process and finally assessing the repeat interval of the phenomenon. By some measures, tsunamis are more extreme hazards than earthquakes, hurricanes and tornadoes. They occur less frequently, and, except possibly in Japan, historic records are unsystematic. Smaller tsunamis have highly localized impact, and tsunamis in the past centuries probably were under-reported and unrecognized, as the recent once-a-year, on average, worldwide incidence suggests. While in earthquake and atmospheric sciences, inferring gross estimates of the repeat period preceded any understanding of the basic physical phenomena, the evolution of tsunami hydrodynamics is an example of extreme natural hazard research where progress in developing modelling tools preceded developments in probabilistic hazard assessment.
The reason has been the lack of tsunami measurements in the open ocean, only accomplished in the later 1990s, by the National Oceanic and Atmospheric Administration's (NOAA) tsunameters, discussed by Bernard et al. (2006). Without instrumental recordings, there has been no widely applicable tsunami magnitude scale, but only qualitative intensity scales, e.g. the six point Sieberg–Ambraseys. Tsunami science has evolved through a combination of near hits and near misses and worse. While the basic governing equations have been known for over 150 years, without instrumental recordings of tsunami signatures, even the appropriate idealized initial conditions for model initialization were not well understood. Until very recently, the only available recordings were from tide gages, which are known to only represent arrival times and possibly the character of the first wave reliably. Tide gauges are located in sheltered regions inside harbours; once the basin is forced, they just record the basin excitation, not the offshore tsunami waveform. Tsunameters are free-field instruments, generally less dependent on location response.
The grand synthesis had to await the development of sophisticated modelling tools, high-resolution laboratory experiments in the 1980s and the field survey results of the 1990s, which served as crude proxies to free-field recordings and allowed for validation of the numerical models. Tsunami forecasting was only made possible through the availability of instrumental tsunameter recordings, which also allowed for closure, e.g. Titov et al. (2005a,b). Starting with hindsight, the process would have been entirely different, but as is most often the case in science, incremental Brownian-moving advances in the entire web of knowledge are eventually integrated. This was largely the case just before the horrific catastrophe of 26 December 2004, whose only silver lining has been the impetus to pursue the asymptotic limit even more vigorously.
Tsunamis are long waves of small steepness generated by impulsive geophysical events of the seafloor and of the coastline, such as earthquakes, submarine and aerial mass failures, volcanic eruptions and asteroid impacts. Tsunamis undergo substantial deformations as they propagate over bathymetry or when evolving over the continental shelf. The determination of the inundation and run-up of tsunamis and forces on coastal structures is one of the quintessential problems in tsunami hazard mitigation. Operational and mitigation issues as well as several relevant tsunami statistics are discussed by Bernard et al. (2006). For a detailed summary of field observation of the 2004 megatsunami refer to Synolakis (2006). Here we will present a short history of the evolution of tsunami hydrodynamics and draw some conclusions on lessons learned and unlearned, and then identify the remaining challenges. We will not comment on complementary harbour excitation from tsunamis; for this, refer to Raichlen (1966) or Synolakis (2002).
2. A short history of tsunami hydrodynamics before 1980
The study of the propagation, run-up and inundation of tsunamis mirrors to a certain extent the evolution of coastal engineering and in some cases has driven it. Early efforts were restricted to one propagation direction, with subsequent attempts directed toward solving the more physically relevant two-dimensional problems. Interestingly, but perhaps not entirely surprisingly, given the differences in cultures between applied mathematicians and hydraulic engineers, the laboratory studies proceeded independently from advances in analytical solutions, and both almost entirely separately from geophysics. The initial model for tsunamis evolved from a bore to a solitary wave, a single wave with the particular property that it does not change its shape when propagating over constant depth. Without any direct measurements of tsunamis in the open sea, there was little realization as to the character or the scale of tsunamis, with few attempts to couple the seafloor motion with the generated wave.
Despite the tsunami catastrophes following the 1946 Aleutian tsunami locally at Scotch Cap, Unimak, Alaska, and at Hilo, Hawaii, and the inundation of Petrapavlosk following the 1952 Kamchatkan tsunami, studies of impulsively generated waves had focused on explosion-triggered water waves, and many remained classified for decades. Numerical studies were non-existent. Laboratory and analytical studies of long waves in the 1950s contained only cursory references to tsunamis. Yet, they did set the foundations for the development of the basic numerical modelling tools.
As a preamble, and to introduce the nomenclature, the equations that govern the motion of tsunamis, as all incompressible fluids, are the Navier–Stokes (NS) equations. These three equations reflect Newton's conservation of momentum, in three directions, and they equate the density times the material acceleration of an infinitesimal fluid particle to the gradients of pressure and viscous stresses, and external body force fields. Along with the equation of conservation of mass, they form four coupled nonlinear equations for the three velocity components and pressure. They are notoriously difficult to solve for free-surface geophysical flows, and approximations are introduced depending on the particular scale of the problem considered. When calculating the evolution of tectonic tsunamis, the most common approximation is the elimination of viscous stresses and of any flow gradients in the vertical direction, given that tsunamis are long waves.
Then, the NS equations are averaged over the flow depth, from the seafloor to the free surface, and the pressure is eliminated directly as a dependent variable and the wave amplitude η introduced through the hydrostatic approximation for the pressure p=ρg(total depth)=ρg(η+h). In the lowest order of approximation, over a flat seafloor, the three resulting equations reduce to one, , the classic wave equation, with the wave celerity. When the seafloor varies, as is always the case, the resulting equation is known as the linear shallow water wave (LSW) equation. In the next level of approximation, nonlinear terms are kept, and the resulting equations are known as the nonlinear shallow water wave (NSW) equations. In the next level, dispersive terms are kept, and the resulting equations are known as the Boussinesq equations. When the viscous stresses are eliminated from the NS equations, the resulting forms are known as the Euler equations and sometimes as potential flow models. Depth-averaged models that study evolution in one propagation direction are referred to as 1+1, and similarly 2+1 models refer to two-directional evolution. For a comprehensive exposition of the orders of approximation through a perturbation expansion, refer to Peregrine (1967).
The first study of a single long wave approaching a sloping beach appears to be that of Hall & Watts (1953). They were the first to employ what is now known as the 1+1 canonical problem of one-dimensional long-wave theory, i.e. long waves generated over a constant depth region, evolving over constant depth, and then climbing a sloping beach, as shown in figure 1. The analogy is of waves propagating over a flat ocean and then climbing up a uniformly sloping continental shelf. The water waves they generated empirically in the laboratory resembled the classic shape of solitary waves, as described by Russel (1845). Using dimensional analysis, they normalized the offshore height H of the wave with the offshore depth d and related the normalized run-up to the run-up R normalized by the depth d. They concluded that , where α(β) and f(β) were empirical coefficients dependent on the beach slope. They did not distinguish between breaking and non-breaking waves, thus haunting investigators up to three decades later.
Wiegel (1955) investigated landslide generated waves in the laboratory. He modelled landslides using a wedge-shaped box sliding down a plane beach, thus establishing the standard model of a moving block still used in many landslide wave studies, but only measured the waves propagating offshore. He concluded that only about 1% of the energy in the landslide motion is ‘converted’ to energy in the generated waves.
While researching hodograph type transformations in gas dynamics, Carrier & Greenspan (1958) derived a seminal 1+1 nonlinear transformation (C+G, for short) that transformed the two NSW equations into a single linear Bessel-type equation solvable with standard methods. The transformation involved changing the physical variables x, t into the set σ, λ, with the peculiarity that, in the σ, λ space, the evolving shoreline in (x, t) remains fixed at σ=0. Their results were for monochromatic waves approaching from infinity. Carrier & Greenspan (1958) did not back calculate the run-up in terms of the physical variables (x, t), and their work remained largely unrecognized and unused outside applied mathematics for decades.
Fortuitously, in the same year, Whitham (1958), while working on shock wave propagation through non-uniform flow regions, suggested a method to circumvent the difficulty of evaluating differential relationships at a discontinuity. He applied the differential relationships to be satisfied by the flow variables on the shock wave, immediately behind it, a practice known in hydrodynamics as the Whitham bore rule. Whitham (1958) then derived approximate formulae of the variation in height and the strength of the bore, a relationship later re-affirmed numerically. Shen & Meyer (1963), extending earlier results, were able to calculate an approximate bore path after the bore collapsed. While not recognized by the time, the bore rule proved useful in the development of inundation algorithms.
Using the LSW equations, Keller & Keller (1964) first derived a Bessel-type equation for monochromatic wave 1+1 evolution over a sloping beach, which peculiarly was the same as the one derived through the G+C transformation for the NSW equations. They then matched this inner solution to an outer solution of periodic waves over constant depth at the transition point, i.e. the toe of the beach, and derived an amplification factor for monochromatic waves.
By the time the Great Chilean tsunami of 1960—triggered by the largest earthquake ever instrumentally measured—struck the Chilean coastline, inundated Hilo, Hawaii, once again, and then struck Japan, little progress beyond calculating estimated arrival times had been accomplished. In the study of the Amorgos, Greece tsunami of 1956, Ambraseys (1960) had demonstrated how simple linear wave theory can be used to derive isochrones, lines of equal travel time from the source. A similar analysis was probably used to estimate arrival times in Hilo in 1960, and there was a warning. Yet, casualties were not avoided, as the first wave was not the largest, triggering a false sense of security, quite reminiscent of what was observed during the 26 December 2004 megatsunami. Most contemporaneous tsunami observations from 1960 were interpreted as tsunamis resembling turbulent bores, fluid structures notorious for the difficulty in calculating their evolution, thus further discouraging any quantitative analysis.
One fundamental observation linking the source with inundation, and never published as such, was George Plafker's. While surveying the aftermath of the 1964 Great Alaskan tsunami, he proposed that the maximum run-up locally does not exceed twice the height of deformed seafloor offshore. This is now referred to as the Plafker rule, and has largely withstood the test of time, remaining unrecognized for three decades (Synolakis & Okal 2005).
Otherwise during the 1960s, tsunami studies focused on refining earlier results. Raichlen (1966) described the problem of harbour resonance for a rectangular geometry and synthesized earlier results within a systematic framework, still used today. In resonance problems, it is the period of the incoming tsunami that is of primary interest, hence progress was possible without precise offshore recordings, or shoreline computations. Tide gauge recordings existed allowing improvement of harbour models much more rapidly than inundation models. Camfield & Street (1966) undertook more comprehensive studies of the evolution of solitary-like waves in the laboratory. Carrier (1966) matched an inner solution using the C+G formulation of the NSW equations to an outer solution at the transition point of the canonical problem. While Carrier boldly attempted to calculate the run-up of a fairly arbitrary disturbance coming from offshore, his analysis did not include reflection. The work did not receive attention until two decades later, as the inversion from the transformed space appeared formidable.
The main advance in tsunami hydrodynamics in the 1960s came with the work of Peregrine (1966) who was able to calculate the evolution of a 1+1 undular bore over constant depth. Peregrine defined for the first time a bore as a transition between two different uniform flows of water, and wrote ‘in its most common form, a bore is a turbulent, breaking zone of water whose length is a few times the depth’. Continuing, Peregrine (1967) derived the different approximations of depth-averaged equations, as widely used today, thus systematically establishing the basic mathematical theory of tsunamis.
By the end of the 1960s, the analytical results remained disjointed to the laboratory experiments on periodic waves. There was still little realization of what tsunamis resembled, beyond that they were long water waves. In the 1970s, there were several substantial advances, including attempts to link the seafloor motion to tsunami evolution.
The worldwide hazard posed by near-field earthquakes in many coastal zones was not fully appreciated until the 1990s, and the paradigm was a transcontinental tsunami striking Japan, Hawaii or the west coast of the US. Most studies in this decade focused on the waves propagating far field away from the source, without really examining the near-field evolution. Ben-Menahem & Rosenman (1972) used linear theory to calculate the two-directional radiation pattern from a moving source. They were able to demonstrate, that given that the velocity of fault rupture was more than 10 times greater than the speed of propagation of tsunamis even in deep water, tsunami energy was shown to radiate primarily at right angles to a rupturing fault. They derived elegant results showing the dependence of the directivity on the seafloor rupture and wave celerity, and proceeded to explain the path followed by the 1964 Great Chilean tsunami. One example of such calculations is shown in figure 2, where simply on the basis of directivity, it was possible to quickly eliminate a shorter source for the Boxing Day tsunami proposed in its immediate aftermath.
Another kind of coupling of the seafloor with the free surface was investigated in 1+1 propagation by Tuck & Hwang (1972). They derived a solution to a forced wave equation for a deforming sloping beach, and suggested that it represented a fault rupture on a continental shelf. No inundation results were shown as their attention was focused on the waves propagating offshore, as their initial condition was an exponential decaying from the shoreline seaward. While the Ben-Menahem & Rosenman (1972) study has been widely used in seismology, as it allowed for a quick but crude evaluation of the tsunami path, the latter study was not rediscovered until the 1990s.
Up until the late 1960s, waves were generated in the laboratory using a horizontally moving vertical paddle. In a landmark study, Hammack (1972) used a novel generation method in the 30 m long, 76 cm wide wave tank at Caltech. One end of the channel had a short section that could be impulsively raised or down-dropped, to model the seafloor motion that triggers tsunamis. Hammack measured very precisely the evolution of the resulting waveforms over constant depth. He used the laboratory measurements to validate Peregrine's (1966) solution method for the nonlinear and dispersive Korteweg–de Vries equation, one form of the depth-averaged equation of motion. Hammack then related the initial wave to the wave motion at large distances, by using the inverse-scattering (INS) algorithm. The INS methodology developed by Gardner et al. (1967) for solving certain nonlinear equations with soliton solutions, predicted that fairly arbitrary initial conditions would generate a sequence of solitary waves at infinity, also allowing for the calculation of the leading soliton height.
Once it was clear, or so it seemed, that arbitrary disturbances would eventually fission into series of solitary waves, it was not necessary to fully understand the initial condition. If far enough away, the source details were irrelevant. The focus shifted to evaluating the transformation of solitary waves as they propagate over the continental shelf. This was accomplished by Goring (1978) who not only developed an algorithm for driving a laboratory wavemaker to generate precise solitary waves, but also developed a finite-element solution to the 1+1 Boussinesq equations and validated it with the laboratory data, thus confirming earlier theoretical predictions by Madsen & Mei (1969), suggestive of soliton fission at transitions in depth.
At about the same time, Houston & Garcia (1974a,b) in an attempt to develop flooding maps for the Federal Emergency Management Agency (FEMA) calculated the 2+1 deep-ocean propagation of tsunamis from Chile and Alaska evolving towards the west coast of the United States. They used a hybrid method and first solved a linear form of the spherical long-wave equations to propagate the tsunami from the source to the edge of the continental shelf, using a finite-difference model; at the continental shelf, they matched the outer amplitudes to an ‘equivalent’ sinusoid and determined a simple amplification to determine wave heights. Houston (1978) further refined the method to hindcast the tsunami response of the Hawaiian Islands, using steady-state models. Bernard & Vestano (1977) had earlier determined the transient excitation of the islands.
The Houston & Garcia work was groundbreaking, also because Houston & Garcia attempted to use more geophysically realistic initial conditions. They argued that the only reliable data for defining source characteristics at that time were from the 1964 Alaskan and the 1960 Chilean earthquakes. Based on these data, they approximated the initial ground deformation by a 1000 km long hypothetical uplift mass of ellipsoidal shape, with a 10 m maximum vertical uplift. They then divided the Aleutian trench into 12 segments and calculated the wave evolution from each segment, and repeated the procedure for tsunamis from the Peru–Chile trench, convolved their results with tidal cycles and thus developed the first ever probabilistic tsunami hazard analysis. While their predictions for the offshore wave height were repeatedly interpreted as run-up, and were later shown to differ with detailed run-up computations, in some cases, by a factor of two (e.g. Synolakis et al. 1997a,b), the scope of the work was pioneering.
Hibberd & Peregrine (1979) brought back the focus to coastal effects. They proposed an entirely new shoreline algorithm to calculate the evolution of strong bores over a sloping beach and then computed run-up using the NSW equations. This was a spectacular advance. The wave may break during the run-up over the dryland, bores collapse, and the moving front is thin with high initial velocities, a situation which usually results in instabilities in numerical computations. This last run-up step had been entirely missing in earlier numerical calculations of 1+1 propagation. Models either ignored it altogether and interrupted the computation at some threshold offshore location, from which point inundation was inferred—hence their name as threshold models—or attempted shoreline computations irreproducibly. Run-up calculations are still missing even in most current 2+1 propagation models owing to the perceived complexity when applied over physical bathymetries. This is quite regrettable, as some studies have suggested that run-up inferences from threshold models may differ from models with shoreline algorithms by factors up to two, and sometimes even more. The Hibberd & Peregrine evolution computation was ad hoc and conceived for bores that evolve more rapidly than non-breaking waves; this analysis allowed for the development of all shoreline algorithms used since then.
By the end of the decade, there was still little realization of the scales of tsunamis of geophysical interest, or the differences in the run-up variation for breaking and non-breaking waves. Even the term run-up promulgated since Hibberd & Peregrine's work, up to this time, tsunami vertical rise (and occasionally tsunami run-up) had been used interchangeably with tsunami height. These issues were addressed in the 1980s, and many of the earlier results were synthesized into powerful analytical and numerical tools, validated with laboratory experiments, setting thus the stage for the advances in the now widely used tsunami forecast models, as developed in the 1990s.
3. Tsunami hydrodynamics in the 1980s
Early in the decade of the 1980s, and motivated from Hibberd & Peregrine's analysis, numerical methods for calculating run-up started to rapidly develop. The focus remained on solitary waves as initial conditions, except one approach pursued by geophysicists and yet to be reconciled in classic hydrodynamic practice.
Ward (1980) introduced the concept of tsunami waves as a particular branch of the spheroidal family of normal modes, or free oscillations, of the Earth. This formalism is traceable to Love's (1911) famous monograph for the computation of the eigenfunctions of the Earth's modes. Ward asserted that it is directly applicable to the case of tsunamis, as long as the Earth model includes an oceanic layer, and the computation is carried in the full six-dimensional space allowing for gravitational effects. As discussed later by Okal (1982, 1988), this allows the direct and seamless handling of the coupling between the ocean water layer and any vertically heterogeneous solid Earth structure. While the normal modes methodology is not limited to long water waves, it is by nature linear. It also requires a spherically symmetric structure involving an ocean covering the whole planet with a uniform depth d. It remains largely unused in inundation modelling, and its veracity is not widely recognized outside geophysics.
Given the extensive work in the 1970s on the evolution of solitary waves, the prevailing paradigm as to the leading wave of tsunamis, their run-up was revisited 30 years later by Pedersen & Gjevik (1983). They numerically solved a 1+1 Boussinesq equation in Lagrangian coordinates. For a given solitary wave, they compared their results with analytical predictions for a sinusoid that best approximated the shape of a solitary wave far from its tails. Based on this sinusoid, they derived a breaking criterion for the run-down of the wave and attempted comparisons with their own laboratory data. Kim et al. (1983) also solved a potential flow 1+1 equation and compared their results for the maximum run-up with some of the Hall & Watts (1953) data, with varying degrees of agreement.
By the later 1980s, in a series of Ph.D. dissertations at Caltech, significant advances were achieved in exploring the characteristics of breaking and non-breaking solitary waves and their run-up. Synolakis (1986, 1987) solved the initial-value problem (IVP) of the NSW for the evolution and run-up of solitary waves first propagating over constant depth and then evolving over a sloping beach. To predict the run-up of wave of the classic Boussinesq solitary wave shape, he first used the Keller & Keller (1964) formalism to derive a boundary condition at the transition point between the sloping beach and constant depth region. He then linearized the C+G transformation at the transition point, arguing that far off the beach nonlinear effects were expected to be small. Synolakis (1986) then proceeded to derive a new solution of the C+G equation for general offshore initial conditions and calculated the evolution of a solitary wave with nonlinear theory, through the entire run-up and run-down process. As opposed to Carrier (1966), he was able to back calculate the evolution of the wave in the physical (x, t) space, and successfully compared the analytical predictions with laboratory measurements in the 30 m long, 76 cm wide Caltech channel.
This methodology was used to demonstrate the run-up invariance between maximum run-up results from the nonlinear and linear forms of the Shallow Water wave equations. This invariance was first suggested by Carrier (1966), without accounting for reflection off the beach which can drastically change the wave shape. Synolakis was able to demonstrate that for the same initial wave offshore, and despite the fact that linear and nonlinear theory predicted different wave evolutions nearshore, their predictions for the maximum run-up were mathematically identical. The key advance was the development of a mathematical formalism that allowed for direct contour integration of the resulting Hankel transforms and the derivation of McLaurin series that could be summed up asymptotically, in parameter ranges of geophysical interest. For example, for the run-up of solitary waves of offshore height H, propagating over depth d and then climbing up a beach of angle β, with profile , he derived asymptotically that(3.1)where the coefficient and the exponent are exact. By calculating the Jacobian of the C+G transformation, he derived a breaking criterion to identify among his own run-up data and Hall & Watts' (1953), the non-breaking solitary waves and compared the predictions from the asymptotic result (3.1) to the data, as shown in figure 3. It appeared, that at least the problem of scaling solitary wave run-up was solved and the differences in functional variations for the run-up of breaking and non-breaking waves were identified. Most later work with numerical solutions for solitary waves has used (3.1) for validation.
Another result whose usefulness was not entirely recognized until the wide availability of amateur videos of the 26 December 2004 megatsunami was the acceleration of the wavefront past the initial shoreline. As shown in the evolution of the front of the solitary wave in figure 3, the wave first slows down as it climbs up the sloping beach, due to the reduction in depth, as the wave celerity , where h(x) is the changing depth. As the wave hits the initial shoreline, it slows down considerably, but accelerates again, before decelerating to its ultimate onland penetration point, the maximum wave run-up. In a video taken near the Grand Mosque in Aceh, one can infer that the wavefront first moved at speeds less than 8 km h−1, then accelerated to 35 km h−1. The same phenomenon is probably responsible for the mesmerization of victims during tsunami attacks, first noted in series of photographs of the 1946 Aleutian tsunami approaching Hilo, Hawaii, and noted again in countless photographs and videos from the 2004 megatsunami. The wavefront appears slow as it approaches the shoreline, leading to a sense of false security, it appears as if one can outrun it, but then the wavefront accelerates rapidly as the main disturbance arrives.
Skjelbreia (1987) measured flow velocities in breaking solitary waves generated directly on a sloping beach using a laser Doppler velocimeter. The data allowed for testing of approximate long-wave relationships for the horizontal and vertical velocities with the wave amplitude. Away from the breaking point, even steep waves were shown to have an approximately uniform horizontal velocity structure, up to the boundary layer, as anticipated by SW theory. In the same laboratory of Fred Raichlen at Caltech, Zelt (1991) developed a Lagrangian solution of the 2+1 Boussinesq equation and calculated the run-up distribution of a solitary wave entering a parabolic basin.
4. The 1990 Catalina long-wave run-up workshop
The decade of the 1990s was launched with two landmark scientific meetings, one in Novosibirsk, Russia in 1989 and another in Twin Harbors, Catalina Island in the US in the 1990. While the Novosibirsk symposium organized primarily by Gusiakov (1990) included mostly geophysicists and applied mathematicians, the Catalina workshop brought together applied mathematicians and coastal engineers and was geared more towards understanding the state-of-the-art of tsunami hydrodynamics.
In the National Science Foundation (NSF) workshop in Catalina, key existing results were summarized and discussed (Liu et al. 1991). The emphasis was on examining how to best validate numerical and analytical tools under development. E. Pelinofsky demonstrated, albeit with a different transformation, the run-up invariance between linear and nonlinear theories. C. E. Synolakis discussed the generalization of Green's law (1833), a classic result for the amplitude evolution of periodic waves on sloping beaches, to solitary waves; far from the shoreline, ηmax the maximum height of solitary waves of height H evolves with the inverse quarter power of the depth, even with reflection included, as predicted by Green's law, which did not include reflection.
D. H. Peregrine presented novel computations of colliding breaking waves and inferred that the associated fluid accelerations when a breaking wave collapses on a vertical wall can reach 1000g, and noted that the associated high pressures were such that compressibility effects remained unimportant, at least for this type of impact. Peregrine then presented numerical results of the evolution of groups of bores whose run-down was interacting with the run-up of the subsequent front, a fiercely feared computation. C. C. Mei presented new results with a multiple-scale perturbation theory to study harbour resonance. H. Yeh presented laboratory measurements of a bore climbing a sloping beach, taken with a high-speed camera, and observed that the bore collapse was not as anticipated in earlier analytical predictions. Ramsden & Raichlen (1990) presented results for the impact pressure of a bore on a vertical wall.
In terms of numerical modelling, V. V. Titov presented the Novosibirsk's Computing Center's finite-difference 2+1 algorithm that allowed for efficient calculations of the far-field evolution of tsunami-looking waves, a computation that whose propagation algorithm became the basis of the code MOST (Method of Splitting Tsunami; Titov & Synolakis 1998) now used by NOAA for tsunami inundation forecasting. N. Shuto presented a spectacular animation of a long wave flooding rectangular shaped structures on dry land which, however, did not include actual inundation calculations. Z. Kowalik presented a 1+1 computation for bore run-up. S. Grili and I. A. Svendsen presented what they referred to as exact numerical computation using the boundary-element method for 1+1 long-wave evolution, up to its maximum run-up but no further, solving Euler's equations. P. L.-F. Liu presented an animation of a three-dimensional potential flow solution for a wave sloshing in a rectangular basin.
F. Gonzalez presented the first ever measurement of waves with a period greater than 5 min in the open ocean at a resolution of 1 mm, with a pressure gauge at 6000 m depth. This instrument was the basis of what is now known as tsunameters or Deep Ocean Assessment and Recording of Tsunamis (DART) buoys (Bernard et al. 2006). E. N. Bernard emphasized the link between modelling and civil defence strategies and suggested the importance of both site-specific and source-specific inundation maps that would depict the flooding hazard from scenario events.
With the exception of the latter presentation, it would have been clear to anyone not involved in tsunami research that despite the advances and the scientific talent involved with them, the understanding of real tsunamis was not altogether entirely different from 20 years earlier with improvements only in the least realistic problems, which could not rapidly lead to tsunami forecasting or inundation mapping. The main conclusion of the Catalina workshop was that the run-up of a single non-breaking wave could be computed analytically and numerically. The NSW model was determined to be adequate for the applications of geophysical interest. There was consensus of a pervasive need for laboratory data for long waves propagating in two directions to allow further progress in computational models particularly for 2+1 run-up computations. There was no realization of the importance of the seafloor/water-surface coupling, and the solitary wave model remained the standard. All this was to change.
5. Tsunami hydrodynamics in the 1990s
The decade of the 1990s showed rapid progress towards a more realistic understanding of the entire tsunami evolution process from generation to run-up. Not only did large-scale two-directional laboratory data became available, but also tsunamis started to be fortuitously reported at a rate of almost one per year, debunking their extreme hazard status, at least for moderate events. Dr Cliff Astill of the Natural Hazards Mitigation Program of NSF aggressively encouraged field measurements in the aftermath of tsunami catastrophes to generate datasets to be used for code validation. The field surveys and model-validation workshops brought applied mathematicians, coastal engineers, geologists and geophysicists closer together, and advances and results were absorbed and implemented across fields with increasing vigour (Synolakis & Okal 2005). NOAA started deploying tsunameters, and the first measurement of a tsunami in the open ocean was made possible. By the end of the decade, some regions in the US and Japan had inundation maps derived with the state-of-the-art modelling tools based on actual tsunami scenario events and with physically realistic initial conditions. Further, various paradoxes associated with earlier results were explained. Finally, with the Papua New Guinea (PNG) tsunami of 1998, the importance of landslide-triggered tsunamis was realized and a significant effort materialized to better understand landslide waves within the context of long-wave theory.
The tsunami decade was launched with the 1 September 1992 Nicaraguan tsunami. It was the first of 15 surveys carried by the International Tsunami Survey Team (ITST). With more than 160 people killed, this was the first major tsunami to strike in 9 years. This and its clear nature as a tsunami earthquake (Kanamori 1972), motivated a number of detailed surveys and studies. A tsunami earthquake is seafloor motion whose resulting tsunami has far greater amplitude than expected from conventional magnitudes, while its characteristic slow motion is often not reported felt, as was the case with the M0=3.4×1027dyn cm parent earthquake in Nicaragua. Only the shoreline recession was a precursor to coastal residents, who regrettably did not identify it as such, similar to what happened in the 2004 Boxing Day megatsunami, far field from the source. Run-up values along the Nicaraguan coastline ranged up to 11 m in El Transito, the most devastated locale, yet in the adjoining Playa Hermosa, even the beach umbrellas had been left standing, as widely noted then in local newspapers.
Without due regard to the nature of signals from tide gauges, and using established tools for inverting seismic wave time-histories, seismologists attempted hydrodynamic inversions to better define the seafloor motion, given that little is known as to the quantitative differences in the distribution of seafloor displacement of such events compared with more conventional earthquakes. The underlying assumption was that the seafloor/free-surface coupling is instantaneous, so once the initial seafloor surface was determined through inversion, the seafloor deformation becomes explicit. A less appropriate assumption in some inversion models was that the water wave propagation is as linear as seismic wave propagation, discounting the pervasive nonlinearities that inadvertently arise in the extreme nearshore. There was no discussion as to how many records were needed to regularize the inverse problem which otherwise is mathematically ill posed. This remains an unresolved problem.
In the absence of more sophisticated source models, the seafloor or onland ground displacement of non-tsunami earthquakes is computed through a series of analytical formulae, after Okada (1985). They represent the displacement field on the surface of an elastic half-space, when a dislocation of given direction and size is introduced at a given epicentral depth. The deformation is then transferred to the free surface. This differs from the normal modes approach, which allows the direct handling of the coupling between an ocean water layer of constant depth and any vertically heterogeneous solid Earth structure. To match qualitatively the patterns of observed run-up and tide gauge recordings, complex multi-segment fault models for 1992 Nicaragua were proposed through hydrodynamic inversion. Even so, predictions using seafloor displacements derived through the inversions and using 2+1 threshold type models differed by factors up to 10 from run-up measurements. Since threshold models do not calculate run-up, and given the uncertainty in defining the parent seafloor motion, a debate raged for several years, whether it was deficiencies in the hydrodynamics or the seismologic models that were at fault.
Without any time to absorb the impact of this event in scientific terms, the 12 December 1992, Flores, Indonesia M0=5.1×1027 dyn cm earthquake triggered a large tsunami. Run-up heights ranged from 5 m along the north coast of Flores across from the epicentre, to 26 m in Rangrioko in northeastern Flores. Simulations using the state-of-the-art code of the 1980s now known as TUNAMI-N2, a derivative of the Goto & Shuto (1983) model, again underpredicted the measured run-up by a factor of up to 7.
Flores 1992 is mostly remembered by the catastrophe in Babi Island and overland flow in Wurhing. Babi is volcanic, of conical shape with a diameter at the shoreline of about 2 km. It is situated about 5 km directly north of Flores, in between the epicentral region and Flores Island. While the tsunami attacked from the north, most of the inundation was in the south, normally protected from wind waves. The two fishing villages in Babi were completely inundated and run-up values ranged up to 7 m. Survivors described particularly gruesome scenes with human remains left dangling on trees. Wuhring, a 400 m long, 200 m wide densely populated peninsula, was completely inundated with overland flow of 3 m depth with more deaths per 100 persons than anywhere else other than Babi. Photographs from Wuhring are quite reminiscent of images from Banda Aceh, including the surviving mosque (Yeh et al. 1993).
At the time little attention was paid to Wuhring and the overland flow depth was noted and referred to as run-up, the differentiation between the two was not introduced until later in the decade. However, Babi received considerably more attention. Fortuitously and in response to the Catalina workshop, large-scale 2+1 laboratory experiments had been under way at a 27 m wide, 30 m long and 60 cm deep wave basin of what is now known as the Coastal and Engineering Laboratory. A series of laboratory tests were undertaken by Briggs et al. (1994, 1995). Solitary waves were created by a horizontal wave generator with 60 different 45 cm wide paddles moving independently. Given the Babi Island controversy, a laboratory model of a 7.2 m base-diameter conical island was constructed with a slope angle of 14°.
The experiments demonstrated that once the wave hits the side of the island across from the generator, the crest splits into two waves (see figure 4). These two waves move with the crest perpendicular to the shoreline, propagate around the island and collide behind it, in a spectacular demonstration of constructive interference, (Yeh et al. 1994). Neither long waves of smaller crest length nor periodic storm waves do so. This appears to be the first laboratory visualization to an observed tsunami catastrophe.
The circular island data measurements were then used to refine shoreline 2+1 algorithms under development. Preliminary modelling results were published by Yeh et al. (1994), and more comprehensive analyses by Liu et al. (1995) and Titov & Synolakis (1995, 1998), where the code used by NOAA for real-time forecasting MOST was most comprehensively presented. The experimental data formed the basis of one of the four benchmark problems used in a subsequent NSF workshop for inter-model comparison and code validation.
As experiments were still progressing, the 12 July 1993 Hokkaido-Nansei-Oki tsunami hit northern Japan and devastated the island of Okushiri west of Hokkaido. The parent earthquake had an M0=4.7×1027 dyn cm and field measurements identified extreme run-up values of 31 m at a river gully near Monai. The wave overflowed a 7 m high protection wall and destroyed the town of Aonae situated at the southern most tip of the island. Excellent bathymetry data existed before the tsunami, and hydrographic surveys were able to measure the deformed seafloor to high resolution, leading to a more realistic initial condition, used since as another benchmark test for model validation along with the field measurements. Results from a classic comparison are shown in figure 5.
In 1994, the paradigm of a solitary wave as an appropriate initial condition for analytical tsunami modelling—still widely used today for developing benchmark solutions—was challenged for the first time, and Kajiura's (1963) work revisited, unknowingly. By noting that the first arrival in Nicaragua and Flores were waves that caused the shoreline to recede before advancing, Tadepalli & Synolakis (1994) proposed a model for the leading wave of tsunamis with an N-wave shape, in analogy to dipole waves in gas dynamics. To facilitate asymptotic analysis, they chose a certain class of dipolar waveforms with amplitude given by , where a, b define the locations of the crest and the trough, γs is a wavenumber parameter that defines the spread of the wave by analogy to the solitary wave profile , and is a scaling parameter to allow for direct comparison with solitary waves. Tadepalli & Synolakis (1994) used the methodology of Synolakis (1987) to evolve these N-waves in the canonical geometry and found that N-waves with the trough first run-up higher than N-waves with the crest arriving first. They named the latter leading elevation N-waves (LEN), and the former leading depression N-waves (LDN). They derived asymptotic results for different families of N-waves and showed that LEN waves run-up higher than the equivalent solitary waves.
Their work at the time was controversial, for it was not reconcilable with the solitary wave paradigm. LDN waves were believed to be hydrodynamically unstable, the crest was supposed to quickly overtake the trough, while in laboratory realizations of LEN waves the crest quickly separated from the tail, as expected. The additional reports of the LDN waves striking the south coast of Java, Mindoro Island and Shitokan Island, following three more tsunamigenic earthquakes in 1994, did little to settle the controversy. Initial shoreline recessions were dismissed by many as unconfirmed reports of ‘untrained’ eyewitnesses.
In the summer of 1995, the follow-up of the Catalina workshop materialized at Friday Harbor, Washington. The focus here was validation of computational codes to predict run-up and inundation. The organizers had prepared four benchmark problems with laboratory data on the propagation of edge waves, solitary wave run-up, conical island run-up and the field run-up measurements from the 1993 Okushiri tsunami. They provided the initial conditions but no run-up measurements, except as published already. As perhaps expected, there were several 1+1 computations that predicted the run-up of solitary waves, but only three models were able to reproduce more than one of the 2+1 problems, TUNAMI-N2, Liu et al.'s (1995) code, and VTCS-3, a code that became known as MOST. Further, Kanoglu & Synolakis (1995, 1998) were able to compute the run-up around a conical island using analytical tools.
A weak 2+1 solution to the NSW was proposed by Brocchini & Peregrine (1996). Using a similarity transformation to relate the longshore coordinate y with t, they were able to derive a simple expression for the horizontal velocity which further transformed the 2+1 problem into a solvable 1+1 Synolakis (1987) formulation. Even though valid for small angles of incidence, as is the more realistic case, their formulation remains largely unexplored, possibly because of a perceived difficulty in numerically translating their initial condition. Nonetheless, it remains the only 2+1 analytical solution for benchmark testing.
On 9 October 1995, an M0=1.15×1028 dyn cm earthquake struck Manzanillo, Mexico and generated a moderate tsunami. The parent earthquake was the largest to strike the Mexican coastline since the 1932 series, yet run-up ranged mostly up to 4 m with extreme values of 11 m on steep coastal cliffs. The ITST survey was able to acquire two series of photographs from eyewitnesses. In the first, a family was sitting in a veranda of their hotel in the steep coastal cliffs on the south end of la Manzanilla. As soon as they felt the ground shake, they noted the Manzanillo Bay emptying. They took a photograph, believed to be the first documented observation of a leading-depression N-wave causing shoreline recession. In the other series of photographs, taken by the local tortilla maker in a cart at the square in la Manzanilla, one can see three men running away from the tsunami, which reached a maximum run-up of 2 m on the mild beach fronting the town. What is remarkable is that despite the thin advancing tsunami front, and its small overall size, the eyewitnesses could not outrun it.
The photographs motivated the long overdue reconciliation of receding shorelines with the prevailing INS-based theory and the resulting solitary wave paradigm. Tadepalli & Synolakis (1996) forced the LSW equation with a step-function motion to imitate fault rupture, and showed that N-waves were one particular exact solution for the LSW. Using the INS theory to determine when leading solitary waves would emerge from LDN waves generated from their step-function seafloor rupture, they found that N-waves of geophysical scales would have to evolve more than once around the perimeter of the Earth, to reach close to the idealized infinity of INS theory. LDN waves of steepness less than 0.0001 were found stable over 1000 km propagation distances.
Tadepalli & Synolakis (1994, 1996) suggested that because of the nature of the dipolar seafloor deformation unzipping in subduction zones, LDN waves would strike the adjacent shoreline of the subsiding plate, while LEN waves would move towards the open ocean. Using their asymptotic formulae on a representation of the initial 1992 Nicaraguan tsunami as derived from more evolved seismic models, and a simplified form of the long Nicaraguan shoreline, they showed that they could predict the observed run-up analytically within a factor of two. Geist (1998) confirmed these inferences by demonstrating how subduction zone events generated N-waves and developed approximate formulae based on the asymptotic results Tadepalli & Synolakis (1994).
The event, however, that shook tsunami science on a scale analogous to the 2004 megatsunami, was the 17 July 1998 PNG tsunami. Despite the relatively small size of the parent earthquake with M0=3.7×1026 dyn cm, this tsunami resulted in over 2100 fatalities, officially surpassed in the twentieth century only by the 1933 Sanriku, Japan tsunami. The field survey (Kawata et al. 1999) confirmed exceptional run-up heights, reaching 15 m, but concentrated on a 25 km stretch of coastline outside which the effects of the tsunami were minimal. As documented in detail in Synolakis et al. (2002), the combination of excessive amplitude and concentration of the run-up was quickly recognized as incompatible with any simulation model based on the excitation of the tsunami by a seismic dislocation; in addition, the earthquake did not feature a slow source comparable to those of documented ‘tsunami earthquakes’. Finally, eyewitness reports generally indicated that the tsunami had arrived at least 10 min later than predicted by all acceptable models of propagation.
Based on the identification of a hydroacoustic signal recorded near Wake Island featuring an anomalous combination of amplitude and duration, Okal (2003) proposed that the tsunami had been generated by an underwater landslide, itself triggered by the earthquake with a delay of 13 min. A number of hydrographic surveys, using seismic refraction and remotely operated submersibles later identified a 4 km3 slump contained in a submarine bowl-shaped amphitheatre located 25 km off the coast. This was used as the source of the tsunami in several successful numerical run-up simulations.
The identification of a submarine landslide as the source of the 1998 tsunami resulted in renewed sensitivity for the hazards created by underwater landslides. Before then, the understanding of the mechanics of landslide tsunamis was so lacking, that a simple arithmetic error of two orders of magnitude in an empirical formula had remained unnoticed, leading to a substantial underprediction of the height of the leading wave from a landslide off Palos Verdes, California. With the wrong 0.15 m estimate, many had dismissed the local landslide hazard altogether. The PNG analysis allowed for the resolution of an earlier bitter dispute concerning the landslide trigger of the 1994 Skagway, Alaska tsunami, see figure 6.
The 17 August 1999 Izmit, Turkey and 26 November 1999 Vanuatu earthquakes are believed to have triggered mass movements that generated the larger than otherwise expected observed tsunamis. These events, along with the 13 September 1999 Fatu Hiva, Marquesas aseismic tsunami generated by an aerial slump due to the collapse of weathered cliff, further focused interest on landslide tsunamis, as summarized by Bardet et al. (2003). In the aftermath, marine geology surveys that were published suggested that some continental margins in California, Oregon and even the east Coast of the US and elsewhere were especially susceptible to submarine landslides due to the combination of offshore faults and nearshore submarine canyons with stored sediment, as well as bottom material on relatively steep slopes. As a result, the level of hazard posed by relatively moderate earthquakes (typically at the magnitude 6 level) started to be re-examined upwards, on a case-by-case basis (Borrero et al. 2001, 2004). Figure 7 shows their evolution of a landslide tsunami from a source off Palos Verdes, California.
One vexing question raised by the PNG landslide trigger was whether the extreme coastal evolution of landslide tsunamis could be effectively modelled with the NSW equations, or whether the Boussinesq equations would have to be invoked. Landslide waves are steeper and shorter than tectonic tsunamis and hence disperse more rapidly. Lynett et al. (2003) compared Boussinesq and NSW models with the field measurements and concluded that while the nearshore evolution was different, the run-up predictions did not appreciably differ. In the process, they also differentiated run-up and flow depth measurements in the field observations, Before 2000, the term run-up, the elevation of the most inward penetration of the wave was used to refer to the flow depth as well, introducing further challenges in the interpretations of older field datasets.
6. Progress just before the megatsunami
By 2000, some claimed that as many as one-third of the tsunamis ever historically reported were due to landslides, and earlier events were re-interpreted to account for what appeared then as anomalous. Gutenberg's (1939) aphorism that ‘submarine landslides are to be considered at least as one of the chief causes, if not indeed the major cause of tsunamis’ was rediscovered with a vengeance.
Fortuitously, again, the M0=2.8×1027dyn cm Wewak, Papuan New Guinea earthquake of 8 September 2002, occurred 135 km northwest of the 1998 event and allowed for less speculative revisionism. The moment of the 2002 event was seven times larger than in 1998, yet the tsunami generated was much smaller, with no associated tsunamigenic landslide. The field survey of Borrero et al. (2003) revealed a strikingly different run-up distribution from the 1998 event, figure 8. Okal & Synolakis (2004) fitted Gaussian distributions to the respective longshore run-up measurements. For tectonic (dislocation-triggered) tsunamis, they proposed two dimensionless parameters I1=b/Δu and I2=b/a, where b is a measure of the width of the run-up distribution in the longshore direction; Δu, a measure of the seafloor vertical displacement; and a, a measure of the maximum height of the run-up distribution. For landslides, the invariant I1 is replaced with I3=b/ηmin, where ηmin is the height of the leading dipole wave in the LDN. I1, I2 and I3 behave as invariants; while they characterize the tsunami source landslide versus dislocation, they are largely independent of the exact parameters describing the sources, and can be used as discriminants to identify the nature of the tsunami source. Figure 8 also shows how I1 varies with I2 for 37 simulations of tectonic sources, I3 with I2 for 37 landslide sources, with field data for validation. It is clear that landslide and tectonic sources can be identified on the basis of I2.
A new methodology for modelling landslide tsunamis was pursued by Ward (2001). If the equations of motion are not depth averaged, but only linearized, the resulting formulation involves solving Laplace's equation. Ward re-derived Green's functions and proceeded to describe complex landslide sources through their linear superposition, as discussed earlier by Mei (1983). One advantage is that the formalism is by definition dispersive. A disadvantage is that for some of the extreme waves considered in later work, e.g. volcanic collapse waves of 900 m height in 2000 m depth, these waves are highly nonlinear. Nonlinear effects start to become important in tsunami evolution when the height-to-depth ratio exceeds 0.1, and they rapidly reduce the wave height of larger and steeper waves. Occasionally, dispersion and nonlinearity balance and the resulting wave does not evolve rapidly, but this remains to be established for waves of the size and steepness often considered by Ward in subsequent work.
When available, historical archives can be used successfully, as demonstrated by the brilliant reconstruction of the effects of the 1700 Cascadia tsunami in Japan (Satake et al. 1996; Atwater et al. 2005). In their absence, and to better identify the source mechanism of anomalous events such as the 1946 Aleutian tsunami or the 1956 Amorgos tsunami, interest in surveys of historic events has renewed. The 1946 tsunami was triggered by a deceptively small M0=9×1028 dyn cm earthquake, yet eradicated Scotch Cap in the near field, killed 159 people in Hawaii, far field, and wrought significant damage in the Marquesas, with some reports that it was noted as far as Antarctica. During several field surveys in 1999–2002, Okal & Synolakis (2004) were able to compile a database of 86 far-field and near-field run-up measurements through the interview of 69 witnesses aged 59–89 years, with two fundamental results. In the near field, the maximum run-up was revised upward to 42 m at Scotch Cap, rapidly decayed along the coast of Unimak Island. Their measurements require the involvement of a major underwater landslide, as they indicate an aspect ratio of the longshore run-up distribution I2=6.7×10−4 and a ratio of maximum run-up to seismic slip I1=3.4, both in excess of those theoretically acceptable for any dislocation. In the far field, and based principally on their field measurements at Juan Fernandez Island, they documented a very pronounced directivity, which cannot be reconciled with a landslide source, and requires a substantial dislocation. The 1946 controversy appears to start being reconciled.
In an effort to better understand the generation and run-up of waves from submarine and subaerial slides, Raichlen & Synolakis (2004) conducted large-scale experiments in a 104 m long, 3.7 m deep and 4.6 m wide wave channel with a plane slope (1 : 2) located at one end of the tank. Following Wiegel's (1955) experiments, they first used a 91 cm long, 61 cm wide with a 46 cm high vertical face wedge block, and then semi-spherical and rectangular boxes of equal volume. By varying the weight of the blocks, they were able to vary the initial acceleration, while the initial position varied from totally aerial to totally submerged. The experiments revealed a three-dimensional depression forming over the wedge as motion initiates. The depression forms the leading portion of the LDN propagating towards shore and running up the slope, while an LEN wave propagates offshore. For the aspect ratio of their boxes, they observed that wave generation became rapidly inefficient as the submergence of the blocks exceeded one block height. These experiments have since been used as another benchmark test for model validation.
The forced linear long-wave equation was revisited by Liu et al. (2003) with a sech-shaped block sliding from the initial shoreline. A new exact solution useful for code-validation was developed as shown in figure 9. Numerically, near the initial shoreline, landslides pose a vexing computational challenge, for not only the shoreline retreats and advances, but also the seafloor is deforming at a rate that cannot be ignored in simulations, as when calculating the evolution of tectonic tsunamis. The analytic solution allowed for a direct comparison of predictions and guided shoreline algorithm improvement. The variation of the 1+1 run-up with the triggering fault generation mechanism based on one-dimensional approximations was accomplished by Tinti & Tonini (2005). Physically, realistic fault ruptures on a uniformly sloping beach were found to result in amplification factors in the range of 1–2, implying that even in 1+1 propagation, the Okal & Synolakis (2004) parameter I1<2, precisely what the empirical Plafker rule suggests.
The 1+1 IVP of the NSW was also revisited by Carrier et al. (2003) with a hodograph transformation for a uniformly sloping beach, in contrast to the canonical problem, whose existing IVP solution is for an initial wave at constant depth and then evolving over a sloping beach. The results showed how N-waves not only have different run-up, but different onland velocity distributions. The Carrier et al. (2003) solution includes integrals with singularities and ad hoc linearizations of the transformation, issues addressed shortly thereafter by Kanoglu (2004).
Before briefly describing the surprises and challenges presented so far by the 26 December 2004 event, it is helpful to summarize the state-of-the-art. The last-to-date intermodel and code validation took place again on Catalina Island in 2004. Four benchmark problems were presented. One, the prediction of the shoreline motion from the wedge experimental data; two, the prediction of the run-up distribution measured in a laboratory model of the locale where extreme run-up occurred during the 1993 Okushiri tsunami; three, the prediction of the shoreline evolution of thin landslides and comparison with the Liu et al. (2003) exact solution; and four, the prediction of the of the run-up to be compared with the IVP solution of Carrier et al. (2003).
An entirely novel model of solving the benchmark problems using smooth particle hydrodynamics (SPH) was presented by R.A Dalrymple. SPH describes a fluid by replacing its continuum properties with locally (smoothed) quantities at discrete Lagrangian locations and hence it is meshless. However, only the Synolakis's (1987) and Carrier et al. (2003) benchmark IVP solutions were shown to be amenable to most models presented, whether solving the SW, Boussinesq or Euler's equations, as was the case in Friday Harbor in 1995. A few models were able to predict reasonable results compared with the predictions of the Liu et al. (2003) exact solution for thin landslides. Practically, all 2+1 directional models presented had attempted the Okushiri simulation with varying degrees of success. MOST, TUNAMI-N2 and Lynett et al.'s (2003) methodology appeared to have an edge, with the Boussinesq solution at greater computational cost. Only the code developed by Liu et al. (2005a) appeared capable of modelling the wedge landslide data, but with extreme computational complexity. Figure 10 shows stills from their laboratory experiments, the numerical results and comparisons of measurements with predictions for the shoreline motion. It is quite clear that operational and research applications require different models, but that both require comprehensive validation of the numerical codes.
Overall, substantial progress had been achieved in 50 years of tsunami science. Integrated tsunami inundation forecasts integrating seismic and tsunameter recordings in real time and leading to cancellation of warnings have been attempted successfully (Bernard et al. 2006). Based on a single reading off a NOAA tsunameter in the North Pacific for the 17 November 2003 tsunami, Titov et al. (2005a) scaled their precomputed scenario event that most closely matched the parent earthquake to evaluate the leading wave height off Hilo, Hawaii. They then proceeded with a real-time NSW computation that resulted in the first tsunami forecast and a warning was cancelled. Their comparison of the forecast to the measured wave is the new golden standard for operational models.
The code most used in inundation mapping in the US has been validated through extensive series of benchmark tests. Validation remains an ongoing process, with every new event posing a new, but diminishing challenge. A first generation of inundation maps based mostly on validated 2+1 models now exists for the US and Japan and in some places in South America. A high-resolution probabilistic tsunami hazard map is now available for Seaside, Oregon, integrating near-field and far-field hazards and showing the geographical distribution of different levels of flooding with associated probabilities, including high impact zones. The latter are based on estimates of the normalized momentum flux, and account for the fact that the flow depth is not the only predictor of damage. The impact zones reflect the variation of simple Froude number type parameter V2/gD, where V is the magnitude of the overland flow velocity and D the overland flow depth.
If the tsunami community appeared at first perplexed in the aftermath of the Boxing Day tsunami, it was not due to the failure of hydrodynamic paradigms, but because of the unprecedented loss of life, the worst possible surprise.
7. Lessons learned and unlearned
A recurring question among all pre-2004 Sumatra tsunami scientists has been what else could had been done from the basic science viewpoint to have reduced the horrific disaster. Given that the hazard had not been realized earlier, this deconstructionism is mute. Had the hazard been known, the basic science existed to have produced inundation maps depicting possible flooding zones from a variant of the 26 December 2004 event, everywhere across the Indian Ocean. Widely media-published simulations were rapidly prepared in the immediate aftermath, without the benefit of field measurements for guidance, and eventually published in the scientific literature with little changes (e.g. Smith et al. 2005; Titov et al. 2005b; Geist et al. 2006b). More importantly they appear in overall agreement with satellite observations, underscoring how at least the far-field impact was predictable. Even a crude tsunami warning without tsunameter data might had been taken seriously. Education campaigns might had been undertaken to better acquaint people everywhere with what was a possible, however unlikely, worst-case flooding event. A yet easier question is which specific lessons learned and advances in our understanding of tsunamis might have been useful and how the Boxing Day tsunami changed any of it, if at all.
Possibly the most glaring surprise was the hundreds of pictures in Phuket of tourists just casually watching the onslaught of the tsunami within 100 m of the coastline. Undoubtedly most perished. The images of inundation were practically indistinguishable from photos from the 15 events since 1992. The response of local residents and tourists in 2004 was unfamiliar, at least for post-1998 tsunamis. Neither the pre-existing tsunami folklore, recorded in numerous popular-science books since 1946, nor the telling Manzanillo pictures and their impact in changing the scientific paradigm of tsunamis had reached the wider public, worldwide. Whether the better communication of established scientific results on tsunami precursors was an issue for scientists, media, governments or non-governmental organizations, this is yet unresolved. Mercifully all four stakeholders have made significant progress to avert another catastrophe, as evidenced by the response to the 28 March 2004 Sumatran event; it is estimated that about 75% of the aid pledges for the Boxing Day tsunami have already materialized. If true, this is an astounding figure.
Less surprising as a concept, but not so in terms of impact, were the effects of human modifications and poor land use in enhancing tsunami risk. By their very nature, long waves were not supposed to be so exquisitely affected by features several orders of magnitude smaller than their wavelength. At one extreme, tides flood irrespective of small-scale human modifications; at the other end, storm waves are very responsive to the same features. Many tsunamis resemble a fast receding and then fast moving tide, yet their interaction with coastlines resembles more that of storm waves, particularly on gentle beaches, where the wavefront breaks.
The reef fronting the devastated El Transito during the Nicaraguan 1992 event had an opening to allow for easier navigation, hence its rapid development as a fishing village. The adjacent Playa Hermosa that was largely spared did not. During the 1993 tsunami at Aonae, a manmade dune and about 50 concrete wave protectors channelled the tsunami into the populated portions of the town, while protecting the unpopulated areas. In Sri Lanka, the Sumudra Devi, a passenger train out of Colombo, was derailed and overturned by the tsunami killing more than 1000 (Liu et al. 2005b). In the immediate fronting area, significant coral mining had occurred, related to tourism development (Fernando et al. 2005). In Patong Beach, Thailand, a less than 60 cm high seawall separating the beach from the road reduced impact velocities (Dalrymple & Kriebel 2005). Mangroves were observed to have protected coastal communities in southeastern India (Danielsen et al. 2005).
In the 1994 East Javan and 1996 Peruvian tsunamis, it had been observed that coastal dunes limited the amount of tsunami penetration, although then there had been no human settlements to be impacted. In Karon Beach, Thailand, a low sand dune protected the area behind it. In Yala, Sri Lanka, a resort, for the purpose of better scenic views, had removed some of the dune seaward of the hotel. The hotel was entirely razed to the ground. Substantially larger water elevations and greater damage observations were found in the hotel grounds as compared to neighbouring areas, located behind unaltered dunes. Earlier analytical work was suggestive that in many cases the last topographic slope long waves encountered as they attacked composite beaches affected the run-up to first order, and it was further inferred that other small-scale features do so as well (Synolakis et al. 1997b).
Low lying coastlines have been known to be particularly vulnerable to tsunami attacks, as first observed in Wuhring during the 1992 Flores tsunami. In the context of inundation mapping in California, it had been observed that areas that get severely flooded during El Nino events, such as Seal Beach, feature the longest inundation distances when exposed to scenario tsunamis. Areas where the tsunami can attack from two sides are very prone to severe inundation even from small tsunamis, as observed in 1994 in East Java and in 1996 in both Biak and Peru. Despite the Boxing Day tsunami being the fifth to hit Indonesia in a period of 13 years, the population was unprepared and was decimated in the low lying area between Banda Aceh and Longhka in northern Sumatra, see figure 11.
All over the Indian Ocean, observations of tsunami damage to buildings were unsurprisingly similar. In numerous locations in Sri Lanka, churches and Buddhist temples were left standing. Closer to the beach in Aceh, the only conspicuously standing structures in the devastated wasteland of coastal areas were mosques. In both locales, residents credited divine intervention, but another explanation is that places of worship are more carefully constructed to a higher standard than the surrounding non-engineered buildings. Also, the open architecture of mosques in southeast Asia, featuring multiple circular columns on the ground floor with no surrounding walls, allows the tsunami to flow freely through it, with little impact on the superstructure above them; few if any knew to evacuate to upper floors of surviving structures. Even in Vilufushi, the most severely impacted atoll in the Maldivian archipelago, the local school and city hall, both two-storey structures, survived with no damage (Fritz et al. 2006). Yet, people in surrounding non-engineered houses perished. While the further survival of buildings left standing from the shaking might not be known a priori, when there is no other choice, vertical evacuation is the only choice.
In terms of basic hydrodynamic science, to date, there were four noteworthy phenomena unobserved earlier. One, the manifestation of the Boxing Day tsunami as an LEN in coastlines west of the subduction zone and as an LDN east of the subduction zone. As shown in figure 12, this behaviour was implicit in the work of Tadepalli & Synolakis (1994, 1996) but had not been observed in nature. Two, the waveguide type effect from mid-ocean ridges that appears to have funnelled the megatsunami away from the tip of Africa (Titov et al. 2005,b). Three, the sparing of the Maldives, an archipelago with coral atolls which rise to no more than 2 m at best from mean water level. The islands rise from the seafloor as pillared structures, and there was no significant wave amplification. While the reef fronting the islands determined the extend of inundation, there is little question that the Maldives experienced a tsunami with heights closer to the free-field tsunami height that elsewhere, as could had been inferred from earlier analytic work by Lautenbacher (1979) and Kanoglu & Synolakis (1998). Four, the comparison between 28 March 2005 and 26 December 2004 tsunamis which suggests that the smaller water depth and the presence of the islands of Nias and Simeulue reduced the effective water mass set into motion during the later event, thereby drastically reducing the size and the impact of the generated tsunami. (See figure of Arcas and Synolakis in Kerr (2005) and Geist et al. (2006a,b).)
8. The asymptotic limit
What remains in terms of emergency management and warning is discussed by Bernard et al. (2006) and Geist et al. (2006b). Clearly, tsunameter measurements are not only needed for real-time forecasting and warning, but also to quantify what is still a very inexact science as far as developing tsunami magnitude scales, a most basic step. The 1+1 and 2+1 tsunami wave evolution is by now well understood, with uncertainties arising only from the seafloor–water wave motion coupling that initializes the models. Futile as the exercise of prescribing the progress of science is, the glaring show-stoppers are obvious.
One, tsunami numerical models must continue to improve through a combination of testing with benchmark laboratory data, instrumental tsunameter recordings and field inundation measurements. Specific emphasis needs to be given to landslide initial conditions and sediment layers, a problem identified almost 10 years earlier (Synolakis et al. 1997a,b). The potential hazard from a tsunami source should never again be underestimated by two orders of magnitude, because of the use of untested empirical numerical models. As elsewhere in science, just because a numerical model produces results, it does not imply they are meaningful or always physically realistic. This was painfully demonstrated in the 2005 Fall Meeting of the American Geophysical Union where inundation results were presented for Cascadia with estimates differing by a factor of four from measurements from sedimentologic studies. These `new' unreviewed predictions adversely affected for a few months the credibility and work emergency managers in the Pacific Northwest, who had to explain to their citizenry the differences between peer-reviewed science and the alternatives. Worse, the same model was used to predict future inundation in southern Sumatra. There needs to be a differentiation between research and operational modelling, and hopefully eventually there will not need to be. Until then, validation standards for inundation mapping and tsunami forecasting need to be established, there is now more than ever a greater risk of over-reaction than vice-versa.
Two, the forces on structures need to be better understood, particularly since tsunami floods are really debris floods. Existing methodologies are based on riverine flooding results, and there is little understanding of impulsive or impact loads of small duration or of the uplift pressures from splashing tsunamis. Given the survival of places of worship, and their likely use as shelters in the future, whether planned or not, comprehensive simulations need to establish how safe they really are. Further, most coastal dwellings and marine structures are only designed to survive extreme storm loads. Tsunami engineering stands where earthquake engineering was in 1971, when the San Fernando earthquake in California triggered massive revisions of building codes worldwide. This needs to change.
Three, continental margins in seismic areas around the world need to be mapped at high resolution to allow for better offshore hazard identification. There should never be another incidence of an unrecognized megathrust fault, or an unknown submarine amphitheatre capable of producing a massive slump threatening a nearby populated coastline. Palaeotsunami studies must continue with increased vigour to allow for better identification of the repeat interval of tsunamis to better define risk. Even the Chixulub impact left its trace in the sediment record. Experience suggests that for really extreme events of the magnitude contemplated recently for volcanic flank collapses in the Atlantic, absence of evidence of extensive tsunami inundation is often evidence of absence.
Four, our understanding of the coupling of the seafloor motion with the initial tsunami wave remains marginal in all but the simplest cases. While in the Boxing Day tsunami existing paradigms performed well to first order for far-field propagation, and they appear adequate for tsunami forecasting, they may not perform equally well for nearsource forecasting of extreme events, such as seafloor ruptures over 1000 km long, or submarine landslides. The perceived uncertainty has led to numerous highly speculative mechanisms being proposed for the 2004 seafloor deformation and how it transferred energy to the tsunami, interestingly all claiming excellent agreement with the same dataset that the results of simpler models agree with. Seafloor deformations from slumps and solid slides need to be studied time-coupled with the sea-surface waves they trigger to develop confidence in initializing numerical models to predict the inundation of nearsource coastlines.
Five, the shortcomings in the population and emergency response observed underscore the urgent need for a worldwide educational effort on tsunami hazards mitigation. This was most recently underscored with the 2006 Tonga and Kythira, Greece earthquakes. In both events, the shaking was widely felt in coastal regions, but it did not trigger any spontaneous evacuations, as the residents reportedly were expecting formal warnings. The latter were unlikely to have been timely, it took longer for the earthquake source motion to be characterized by the ITWC and the Observatory of Athens (respectively) than the travel time of the tsunami to nearby coastlines, had one occurred. Due to deep epicentral depths it did not. Further simply educating local populations at risk is not enough. In an era of global citizenship, it is important that everyone can identify the precursors of a tsunami attack and knows to evacuate to high ground or inland as quickly as possible, and if not, how to more safely vertically evacuate to well-built surviving structures.
One contribution of 20 to a Discussion Meeting Issue ‘Extreme natural hazards’.
- © 2006 The Royal Society