UNITED STATES

DEPARTMENT OF COMMERCE

Carlos M. Gutierrez

Secretary

NATIONAL OCEANIC AND

ATMOSPHERIC ADMINISTRATION

VADM Conrad C. Lautenbacher, Jr.

Under Secretary for Oceans

and Atmosphere/Administrator

Office of Oceanic and 

Atmospheric Research

Richard W. Spinrad

Assistant Administrator

NOAA Technical Memorandum OAR PMEL-135

STANDARDS, CRITERIA, AND PROCEDURES FOR NOAA

EVALUATION OF TSUNAMI NUMERICAL MODELS

Joint Institute for the Study of the Atmosphere and Ocean (JISAO)

University of Washington, Seattle, WA

Pacific Marine Environmental Laboratory

Seattle, WA

Pacific Marine Environmental Laboratory

Seattle, WA

May 2007

4

2

Costas E. Synolakis1

Eddie N. Bernard2

Vasily V. Titov3

Utku Kânoğlu4

Frank I. González2

Viterbi School of Civil Engineering

University of Southern California

Los Angeles, CA

1

Department of Engineering Sciences

Middle East Technical University

Ankara, TURKEY

3
NOTICE

Mention of a commercial company or product does not constitute an endorsement by

NOAA/OAR. Use of information from this publication concerning proprietary products

or the tests of such products for publicity or advertising purposes is not authorized.

Contribution No. 3053 from NOAA/Pacific Marine Environmental Laboratory

Also available from the National Technical Information Service (NTIS)

(http://www.ntis.gov)

ii
Contents iii

Contents

1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2. Model Evaluation Standards . . . . . . . . . . . . . . . . . . 3

2.1 Basic considerations . . . . . . . . . . . . . . . . . . . . . . . 4

2.2 Analytical benchmarking . . . . . . . . . . . . . . . . . . . . 4

2.3 Laboratory benchmarking . . . . . . . . . . . . . . . . . . . . 5

2.4 Field data benchmarking . . . . . . . . . . . . . . . . . . . . 7

2.5 Scientific evaluation . . . . . . . . . . . . . . . . . . . . . . . 8

2.6 Operational evaluation . . . . . . . . . . . . . . . . . . . . . . 9

3. Criteria for Evaluating Operational Forecasting and Inundation Mapping Models . . . . . . . . . . . . . . . . . . . . . 9

3.1 Model computational time constraints . . . . . . . . . . . . . 10

3.2 Model accuracy constraints . . . . . . . . . . . . . . . . . . . 11

3.3 Model reliability and realism constraints . . . . . . . . . . . . 11

4. Summary of Model Evaluation Procedures . . . . . . . . . . 12

5. Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . 14

6. References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

Appendix A: Existing Methods for Model Validation and Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

1 Basic hydrodynamic considerations . . . . . . . . . . . . . . . 17

1.1 Mass conservation . . . . . . . . . . . . . . . . . . . 17

1.2 Convergence . . . . . . . . . . . . . . . . . . . . . . 18

2 Analytical benchmarking . . . . . . . . . . . . . . . . . . . . 18

2.1 Single wave on a simple beach . . . . . . . . . . . . 18

2.2 Solitary wave on composite beach . . . . . . . . . . 30

2.3 Subaerial landslide on simple beach . . . . . . . . . 33

3 Laboratory benchmarking . . . . . . . . . . . . . . . . . . . . 34

3.1 Solitary wave on a simple beach . . . . . . . . . . . 36

3.2 Solitary wave on a composite beach . . . . . . . . . 37

3.3 Solitary wave on a conical island . . . . . . . . . . . 39

3.4 Tsunami runup onto a complex three-dimensional

beach; Monai Valley . . . . . . . . . . . . . . . . . . 45

3.5 Tsunami generation and runup due to three-dimensional

landslide . . . . . . . . . . . . . . . . . . . . . . . . 47

4 Field benchmarking . . . . . . . . . . . . . . . . . . . . . . . 48

4.1 Okushiri Island . . . . . . . . . . . . . . . . . . . . . 48

4.2 Rat Islands tsunami . . . . . . . . . . . . . . . . . . 54

Standards, criteria, and procedures for NOAA evaluation of

tsunami numerical models

C.E. Synolakis1, E.N. Bernard2, V.V. Titov3, U. Kˆano ̆glu4, F. Gonz ́alez2

1. Introduction

The National Oceanic and Atmospheric Administration (NOAA) is the federal agency charged with mitigating tsunami hazards in the United States.

NOAA’s National Weather Service operates the two Tsunami Warning Centers (TWCs) in the U.S., and NOAA has spearheaded the national effort to

develop inundation maps for evacuation planning through the National Tsunami Hazard Mitigation Program (NTHMP). The latter was formed through

a directive of the U.S. Senate Appropriations Committee in 1994 to develop

a plan for a tsunami warning system that reduces the risk to coastal residents. Following the horrific Indian Ocean tsunami of 26 December 2004, the

U.S. expanded the role of NTHMP to serve as the organizational framework

to implement the recommendations of the National Science and Technology

council report “Tsunami Risk Reduction for the United States: A Framework for Action” (NSTC, 2005):

• Develop standardized and coordinated tsunami hazard and risk assessments for all coastal regions of the United States and its territories.

• Improve tsunami and seismic sensor data and infrastructure for better

tsunami detection and warning.

• Enhance tsunami forecast and warning capability along our coastlines (Pacific, Atlantic, Caribbean, and Gulf of Mexico) by increasing the number of Deep-ocean Assessment and Reporting of Tsunamis

(DART) buoys, tide gauges, and seismic sensors feeding real-time data

into on-line forecast models.

• Ensure interoperability between the U.S. national system and other

regional tsunami warning systems.

• Provide technical expertise and assistance, as appropriate, to facilitate

development of international tsunami and all-hazard warning systems,

including for the Indian Ocean.

1Viterbi School of Civil Engineering, University of Southern California, Los Angeles,

CA 90089-2531, USA 2NOAA/Pacific Marine Environmental Laboratory, 7600 Sand Point Way NE, Seattle,

WA 98115-6349, USA 3Joint Institute for the Study of the Atmosphere and Ocean (JISAO), University of

Washington, Box 357941, Seattle, WA 98195-4235, USA 4Department of Engineering Sciences, Middle East Technical University, 06531 Ankara,

TURKEY
2 Synolakis et al.

• Encourage data exchange and interoperability among all regional tsunami and all-hazard warning systems, such as the Intergovernmental

Oceanographic Sub-commission for the Caribbean (IOCARIBE).

• Promote development of model mitigation measures and encourage

communities to adopt construction, critical facilities protection, and

land-use planning practices to reduce the impact of future tsunamis.

• Increase outreach to all communities, including all demographics of

the at-risk population, to raise awareness, improve preparedness, and

encourage the development of tsunami preparedness plans.

• Conduct an annual review of the status of tsunami research and develop a strategic plan for tsunami research in the United States.

Since standards for modeling tools do not currently exist, and given the

increased number of states in the Atlantic, the Pacific, and territories in

the Caribbean that will need to develop tsunami mitigation plans, there

is the risk that forecast products may be produced with older or untested

methodology. This is not just a U.S. problem, but worldwide, as UNESCO’s

plans are to build local capacity for developing maps in most member nations

facing tsunami hazards. Incorrectly assessing possible inundation can be

costly both in terms of lives lost, or in unnecessary evacuations in areas

larger than warranted that may put lives at risk and reduce the credibility

of the system, even in areas that were not directly affected. Standards are

urgently needed to ensure a minimum level of quality and reliability for

forecasting and inundation products.

To calculate tsunami currents, forces and runup on coastal structures,

and inundation of coastlines, one must numerically calculate the evolution

of the tsunami wave from the deep ocean to its target coastal community. No

matter what the numerical model, both validation (the process of ensuring

that the model solves the parent equations of motion accurately) and verification (the process of ensuring that the model used represents geophysical

reality appropriately) are an essential part of the model development. Validation ensures that the model performs well in a wide range of circumstances

and is accomplished through comparison with analytical solutions. Verification ensures that the computational code performs well over a range of

geophysical problems. Many existing numerical models have been validated

through comparison with analytical solutions. Very few have also been verified with laboratory and field measurements. Even fewer have been tested

in an operational forecast setting. Numerical models that have been tested

under all three conditions represent the gold standard for both inundation

mapping and real-time forecasting.

In the last ten years, the process of model validation and verification has

shown that terminal effects of tsunamis can be described by depth-averaged

equations. Interest in these equations has arisen because comparisons with

both large-scale laboratory data and field data have demonstrated a remarkable and perhaps surprising capability to model complex evolution phenomena, and in particular the maximum runup and inundation, extremely well
NOAA tsunami numerical model standards 3

over wide ranges of tsunami waves. In the current state of knowledge, the

main uncertainty arises from the ambiguities of the initial condition, presuming that the solution methodology solves the equations of motion satisfactorily. With the increasing deployment of tsunameters that monitor the

tsunami evolution in the deep ocean and allow for real-time updates of the

initial data, this remaining uncertainty is rapidly diminishing.

While equation solvers for more comprehensive approximations of the

parent Navier-Stokes equations now exist, they are presently too computationally intensive for inundation mapping or operational forecasting, and

are generally used for free-surface flows of very limited geographical extent.

These models remain largely unvalidated over wide ranges of tsunami events.

Further, the rapid development of packaged numerical modeling tools that

allow for the production of high-end animations even by untrained users,

underscores the issue of validating specific calculation tools used in tsunami

inundation and forecasting to avoid producing mathematically correct, but

physically unrealistic solutions.

In this report, benchmark tests for validating and verifying computational tools for predicting the coastal effect of tsunamis are discussed. Stateof-the-art methods for validation of tsunami codes are reviewed. Standards

and guidelines for validation and verification are recommended for operational codes used for inundation mapping and tsunami forecasting.

2. Model Evaluation Standards

Tsunami inundation models have evolved in the last two decades through

careful and explicit validation through comparisons of their predictions with

1+1 (one directional and time) and 2+1 (two directional and time) analytical

solutions, laboratory experiments, and field measurements. (See Synolakis

and Bernard (2006) for a comprehensive discussion of the evolution of tsunami hydrodynamics.)

Numerical codes often include terms with “friction” factors in an attempt to model physical dissipation, but primarily to stabilize what is an

inherently marginal inundation computation, given the small flow depths

and large velocities. Large friction factors make calculations more stable,

but may under-predict inundation distances and runup heights. Smaller

friction factors (or just relying on the inherent numerical dissipation of any

numerical calculation) result in more realistic predictions, but they tend to

be less stable, depending on the particular differencing scheme used. Thus it

is recommended that any code used for modeling tsunami inundation at geophysical scales be tested with all three types of validation data: analytical

solutions, laboratory measurements, and field measurements.

While there is in principle no absolute certainty that a numerical code

that has performed well in all the benchmark tests will also produce realistic

inundation predictions with any given source motions, validated codes largely

reduce the level of uncertainty in their results to the uncertainty in the

geophysical initial conditions. Further, when coupled with real-time freefield tsunami measurements from tsunameters, validated codes are the only
4 Synolakis et al.

choice for realistic forecasting of inundation; the consequences of failure are

too serious to take chances with less-validated numerical procedures.

In what follows, the steps for validating 2+1 codes are discussed leaving

the details of the benchmark problems to Appendix A. Although 1+1 codes

were used two decades ago for developing early inundation maps, are generally unreliable for inundation mapping and entirely inadequate for tsunami

forecasting. Nonetheless, analytic solutions to the corresponding 1+1 hydrodynamic equations are invaluable to the process of model validation, and

several are presented in Appendix A. Specific steps are recommended for the

approval of modeling tools, their further development, and their transfer to

operations. These steps can be classified into six categories: basic hydrodynamic considerations, analytical benchmarking, laboratory benchmarking,

field data benchmarking, scientific evaluations, and operational evaluations.

2.1 Basic considerations

Two basic steps are required before benchmarking the numerical codes: mass

conservation and convergence.

Mass conservation: The first basic step in ensuring that a numerical

model accurately simulates tsunami evolution is ensuring that the model

conserves mass. While the conservation of mass equation is one of the equations of motion that are solved in any numerical procedure, cumulative numerical approximations can sometimes produce results that violate mass

conservation. This is particularly true when friction factors or smoothing to

stabilize inundation computations are used. Calculations of conservation of

mass should be such that the total initial displaced volume should be within

5% of the total displaced volume at the end of the computation, i.e., when

the initial wave is entirely reflected and offshore.

Convergence: The next basic step is checking convergence of the numerical code to a certain asymptotic limit, ideally the actual solution of the

equations solved, if one exists. The optimal locations to check convergence

are the extreme runup and rundown. The numerical predictions should be

seen to converge to a certain value, and further reductions in step sizes should

not change the computed results.

2.2 Analytical benchmarking

Why is analytical benchmarking important? To calculate tsunami currents,

forces, and runup on coastal structures and coastal inundation, one must

calculate the evolution of the tsunami wave from deep ocean to its target.

Exact solutions of the shallow water-wave equations are useful for validating the complex numerical models that are used for final design. These

often involve ad-hoc assumptions, particularly during inundation computations when grid points are introduced and withdrawn as the shoreline recedes and advances. Comparisons of numerical predictions with analytical

solutions can identify systematic errors, as when using friction factors or

dissipative terms to augment the idealized equations of motion.

Here, we only present certain common 1+1 propagation problems. The
NOAA tsunami numerical model standards 5

waves evolve over constant-depth and then over plane or composite beaches.

Even though most results are derived for idealized waveforms often used in

tsunami engineering to describe the leading wave of a tsunami, the generalization to realistic spectral distributions of geophysical tsunamis is trivial

with the closed-form integrals provided.

It is important to note that validation should take place with nonperiodic waves. During runup, individual monochromatic waves reflect with

slope-dependent phase shifts. Whereas a code may model a periodic wave

well, it may not model wave superposition well. This was a problem of

earlier shallow water-wave formulations that did not account for reflection.

While their predictions for the Carrier–Greenspan (Carrier and Greenspan,

1958) sinusoids was correct, they exhibited significant errors when modeling

solitary waves or N-waves.

Further, the analytical solutions allow for using cyber-waves of the same

scale as geophysical tsunamis. This is in contrast to laboratory experiments

that require modeling of small-scale experiments.

Analytical solutions for three problems are presented for analytical benchmarking:

Single wave on simple beach: Solitary and N-wave propagation first over

the constant-depth then sloping beach is considered. Details of the analytical solutions are given in Appendices A2.1.1, A2.1.2, and A2.1.3. Numerical

models should calculate the maximum runup of nonbreaking solitary and

N-waves within 5% of the calculated value from the analytical solution. Numerical experiments should be undertaken in a range of depths from 50 cm

to 1000 m, and the initial waves scaled appropriately, i.e., for a H/  ̃  ̃d = 0.02

initial solitary wave; at 50 cm the initial height should be 1 cm, at 1000 m

it should be 2 m. In addition, the initial value problem solution of the nonlinear shallow water-wave equations is described and runup/rundown and

velocity quantities for shoreline are presented in Appendix A2.1.4 for benchmarking. Any well-benchmarked code should produce results within 5% of

the calculated value from the analytical solution.

Solitary wave on composite beach: 1+1 models that perform well with

the solitary wave experiments must still be tested with the composite beach

geometry, for which an analytical solution (Appendix A2.2) exists, with solitary waves as inputs. Numerical predictions should not differ by more than

5% in predicting maximum runup for non-breaking waves.

Subaerial landslide on simple beach: Inundation computations are exceedingly difficult when the beach is deforming as a landslide is occurring.

Numerical predictions of the runup from an idealized landslide of translating

Gaussian shaped mass should not differ from the analytical model by more

than 10% (Appendix A2.3).

2.3 Laboratory benchmarking

Why is laboratory benchmarking important? It is quite clear from the earlier discussion that numerical methods have evolved only through careful

validation accomplished through comparisons with analytic solutions and

laboratory measurements. Long before the availability of numerical codes,
6 Synolakis et al.

physical models at small scale have been used to visualize wave phenomena

in the laboratory and then predictions were scaled to the prototype. Even

today, when designing harbors, laboratory experiments—scale model tests—

are used to confirm different flow details and validate the numerical model

used in the analysis.

For the purpose of validating inundation models, the scale differences

are not believed to be important. Numerical codes developed in the last

decade that consistently produce predictions in excellent agreement with

measurements from small-scale laboratory experiments have been shown to

model geophysical-scale tsunamis well. For example, a numerical code that

adequately models the inundation in a 1 m-deep model is also expected

to model the inundation in the 1 km-deep geophysical geometry, as the grid

sizes are adjusted accordingly and in relationship to the scale of the problem.

Scale models, in general, do not have bottom friction characteristics similar

to real ocean floors or sandy beaches, but this has proven not to be a severe

limitation. Tsunamis are such long waves that bottom friction tends to be

less important than the inertia of the motion. Friction may be important

in cases of extreme inundation, as observed during the 2004 Boxing Day

tsunami in Banda Aceh with 3 km inundation distances. However, it has

been observed that even with numerical codes that use friction factors within

reasonable limits, the predictions are not sensitive to the first order.

Experimental results from five laboratory tests are described as laboratory benchmarking:

Single wave on simple beach: Given a small number of 2+1 wave basin

laboratory measurements, 1+1 versions of the 2+1 models should be tested

with 1+1 directional laboratory models. The solitary wave experiments on

the canonical model—waves propagating over a constant-depth region and

running up a 1:20 sloping beach—should be used first. Numerical models

should calculate the maximum runup of nonbreaking solitary waves within

5% of the measured values in the laboratory. For breaking waves, the models

should produce predictions within 10% of the measured values, and they

should consistently predict the runup variation described in Appendix A3.1.

Solitary wave on composite beach: 1+1 models that perform well with

the solitary wave experiments must still be tested with the composite beach

geometry (see Appendix A3.2 for details). This additional test will ensure

that the code is stable enough for large waves that are near the breaking

limit offshore (H/  ̃  ̃d = 0.696). Numerical predictions should not differ by

more than 5% from the experimental values, and the numerical procedure

should be capable of predicting the entire runup variation.

Conical island: 2+1 dimensional calculations should be tested with the

conical island geometry of 1:4 slope (Appendix A3.3). The numerical method

should stably model two wave fronts that split in front of the island and

collide behind it. Predictions of the runup on the back of the island where

the two fronts collide should not differ by more than 20% from the laboratory

measurements.

Monai Valley: 2+1 numerical computations should then be tested with

the laboratory model of Monai Valley, Okushiri Island, Japan (Appendix

A3.4). The initial condition is a leading depression N-wave, and the entire
NOAA tsunami numerical model standards 7

simulation shows how well the code performs in a rapid sequence of withdrawal and runup. Comparison of results from different codes has shown

that the maximum runup in these experiments can be calculated within

10%, which is thus the standard.

Landslide: As discussed, landslide wave generation remains the frontier

in terms of numerical modeling, particularly for aerial slides. These involve

not only the rapid change of the seafloor, but also the impact of the slide

on the shoreline. Therefore numerical codes that will be used to model

landslide-generated tsunamis need to be tested against three-dimensional

landslide experiments given in Appendix A3.5. It is expected to estimate

maximum runup predictions within 10% of the experiments.

2.4 Field data benchmarking

Why is field data benchmarking important? Verification of a model in a

real-world setting is an important part of model validation, especially for

operational models. No analytical or laboratory data comparisons (or any

limited number of tests, for that matter) can assure robust model performance in the operational environment. Test comparisons with real-world

data provide an additional important step in the validation of a model to

perform well during operational implementation. The main challenge of

testing a model against real-world events is to overcome uncertainties of the

tsunami source. While the source of the wave is deterministic in the controlled setting of the laboratory experiment and can be usually reproduced

precisely, field data always have uncertainties in the source. For tsunamis,

deep ocean measurements are the most unambiguous data quantifying the

source of a tsunami.

Only a very limited number of high-quality deep-ocean tsunami measurements exist and these do not include data for Pacific tsunamis that

have been destructive or caused inundation (Titov et al., 2005). No DART

tsunameters existed in the Indian Ocean at the time of the megatsunami,

since DART had only been deployed in the Pacific Ocean. Satellite altimetry

measurements of the Indian Ocean tsunami do not provide enough quality

and coverage to constrain the tsunami source. Hydrodynamic inversion remains an ill-posed problem and criteria for its regularization are lacking.

Deep ocean measurements allow for more defensible inversions, since they

are not affected by local coastal effects. Several events have been recorded

by both deep-ocean and coastal gages in the Pacific that allow reasonably

constrained comparison with models. The expanded DART system array

will be providing more tsunami measurements for future events, expanding

the library of well-constrained propagation scenarios for testing.

The 10 June 1996 Andreanov Islands (Tanioka and Gonz ́alez, 1998) and

the 4 October 1994 Kuril Islands (Yeh et al., 1995) events were recorded

by several research tsunameters (without real-time data transmission) at

similar locations offshore of Alaska and the U.S. West Coast. The offshore

model scenario for the Andreanov Island event was obtained from the forecast database by inverting data from research Bottom Pressure Recorders

(BPRs) as described in Titov et al. (2005). The inversion of the Kuril Islands
8 Synolakis et al.

data was performed using BPR gages in similar locations. The 17 November

2003 Mw = 7.8 Rat Islands tsunami was recorded in real time by the DART

tsunameter system and provides a good operational test of field data.

High quality bathymetry/topography data, source definition, and field

runup measurements suggest using Okushiri Island, Japan, tsunami data as

one of the field benchmarks (Appendix A4.1). In addition, the Rat Islands,

Alaska, offshore model scenario constrained by the deep-ocean measurement

can be used as input for testing high-resolution inundation models for Hilo

Bay, Hawaii, where tide gage recording provides data for model comparison

(Appendix A4.2).

Okushiri Island: 2+1 computations of the field measurements from the

Hokkaido-Nansei-Oki tsunami around Okushiri Island, Japan (Appendix

A4.1). The initial condition to be used is DCRC-17a (Takahashi, 1996),

which is a composite fault with three segments. The bathymetry and measurements are provided in Appendix A4.1. Predictions for the maximum

runup at Aonae, Okushiri Island, Japan, should not differ by more than

20% from the measurements.

17 November 2003 Rat Islands tsunami: For operational codes, testing

should invert the tsunameter signal of the 17 November 2003 Rat Islands

tsunami to improve the initial estimate of sea surface displacement derived

from a seismic deformation model, then use the results as input to a Hilo,

Hawaii, inundation model to hindcast the tide gage record. This is the most

difficult but most realistic test for any operational model, for it involves

a forecast (now hindcast) and has to be done much faster than real time.

Here, at least the first four waves must be simulated and compared with the

observations, with amplitudes and periods accurate to within 25%, maximum

amplitude accurate to within 10%, and an arrival time error less than 3 min.

Details of the problem are given in Appendix A4.2.

2.5 Scientific evaluation

Peer-review documentation: Model validation and verification is a continuing process. Any model used for inundation mapping or operational forecasts

needs to be presented in peer-reviewed scientific journals with impact factors greater than 0.5. One or more of these publications should include

benchmark comparisons identical or similar to those described above.

Formal scientific evaluation: A formal evaluation process of individual

models needs to be established to avoid ad hoc decisions as to the suitability

of any given model. This process may include solicitation of additional

reviews of the model’s veracity by experts, or the requirement that additional

testing be performed. This process will set the standard for the best available

practice at any given time, and it will hopefully eliminate the liability to

NOAA, universities, states, engineers, and geophysicists who collaborate on

the development of inundation maps.
NOAA tsunami numerical model standards 9

2.6 Operational evaluation

The operational evaluation should be done by a test-bed consisting of research and operational parts of NOAA. While the scientific evaluation process may identify models that are realistic and computationally correct, some

models may still not be appropriate for operational inundation mapping or

forecasting applications. An additional evaluation process must be established to assess a number of model features that bear on important operational factors, such as special implementation hardware/software issues, ease

of use, computation time, etc.

Operational evaluation of candidate real-time forecast models should include the following steps:

Step 1—Meet Operational Forecasting and Inundation Mapping Requirements and Objectives: Operational requirements include: basic forecasting

computation, analysis and visualization tools; integration with operations

(vs. separate, stand-alone applications); basic data assimilation techniques;

computational resources needed to meet milestones; etc. If a candidate

model does not meet specified NOAA forecasting or inundation mapping

requirements and objectives, it should be rejected at this point.

Step 2—Meet Modular Development Requirements: Various pieces of the

forecast model must be developed in parallel based on the overall objectives

defined in Step 1.

Step 3—Meet Test Bed and Model Standards: In this step, the candidate

model is tested against operational standards, with special attention given

to the models ability to simulate previous, major tsunamis with the required

speed and accuracy. Based on these test results, forecast model development

may return to Step 2, proceed, or the candidate model may be rejected for

operational use.

Step 4—Meet Operational Testing Requirements: The candidate model

is integrated into the operational setting for testing. Potential sources are

defined and the model is tested in a forecasting mode on an operational

platform. Graphical interfaces are developed and forecast models are applied

to a few cases to test operational integration and important individual factors

such as speed, accuracy, and reliability (see Section 3, below). Operational

testing and feedback is provided by the TWCs at this point, and adjustments

are made as necessary.

Step 5—Implement Operationally: The model is fully integrated into the

operational setting and procedures to provide operational products.

3. Criteria for Evaluating Operational

Forecasting and Inundation Mapping Models

Given the accumulated experience in the tsunami community in the past 50

years, it is now possible to describe the requirements for an ideal tsunami

model. Given an earthquake fault mechanism and tsunameter data, the ideal

model should accurately predict tsunami inundation of at-risk coastlines in

a sufficiently short time. Sufficiently short is defined as the time interval
10 Synolakis et al.

between the initiation of the tsunami and the calculation of the inundation

forecast that allows for evacuation of the target communities. For example,

the State of Hawaii needs about 3 hr for a complete and orderly evacuation.

An ideal model would accurately forecast inundation at least 3 hr before the

tsunami impact is expected anywhere in Hawaii.

What are the requirements for building such an ideal model, and what

constraints need to be imposed for model selection and/or further development? The following sections define and provide rationale for these constraints.

3.1 Model computational time constraints

The forecast speed will always depend on the distance from the source to the

forecast community. However, an ideal forecast would provide a real-time

forecast which, after being transmitted to local authorities, allows at least a

few minutes before the tsunami impacts the nearest target community. This

time interval would allow sirens to trigger the evacuation of beaches and

coastal residents and give emergency personnel time to mobilize resources

and prepare for search and rescue.

Further, an ideal model would correctly predict the duration of the event.

Tsunamis are a series of waves. Tsunamis often get trapped in closed bays

or on the continental shelf, resulting in sea level oscillations that may persist for several hours. During the 1993 Okushiri Island, Japan, inundation,

bay oscillation at Aonae trapped the tsunami for over 30 min, and a large

portion of Aonae remained submerged for a large portion of this time. The

Crescent City, California, harbor oscillated for more than 4 hr following the

15 November 2006 tsunami. An operational forecast must ideally provide an

estimate of the time that it is safe for search and rescue operations to begin

without endangering the lives of responders.

Computational speed standards for inundation maps and real-time forecasts can be different. The current best-practice standard in the U.S. is

to use a combination of real-time forecast models at coarse resolution and

integrate real-time tsunameter measurements to update the initial sea surface forecast (Titov et al., 2005; Whitmore, 2003). Then a fine-resolution

inundation model, developed beforehand and placed on stand-by, ready to

be run, uses the output of the fast coarse-resolution computation to produce

a quick but physically realistic inundation forecast. The latter is referred to

as a Stand-by Inundation Model (SIM).

A long-term forecast for inundation map production is produced months

or years before a tsunami strikes; computational time is not the limiting step

here. When conducting inundation mapping, the effort is concentrated on

acquiring the best available bathymetry/topography and initial conditions

to produce the most accurate model results. The objective is to provide

guidance for evacuation planning and other hazard assessment products.

Inundation mapping is performed at the highest resolution the available

resources allow. The current standard is 50 m grid resolution, but sitespecific features may demand even finer resolution. Successful computations

at resolutions down to 5 m have been performed.
NOAA tsunami numerical model standards 11

3.2 Model accuracy constraints

The accuracy of any given model depends on how well the computational

procedure represents the correct solution of the parent equations of motion.

When exact solutions exist (as, for example, for certain cases of the linear and

nonlinear shallow water-wave equations), the determination of the accuracy

of a solution algorithm is straightforward, i.e., through comparisons of the

numerical results with the analytical predictions. Determining maximum

runup numerically within 5% of the analytical solution is now possible with

a handful of models.

For most bathymetries of geophysical interest, analytical solutions do not

exist, and it is unlikely that they will ever be determined, due to the complexity of the physical terrain. However, a few laboratory models at smaller

scale than the prototype exist: for example, the scale model of Okushiri

Island, as described earlier. The Catalina Island, Los Angeles, 2004 model

validation workshop of the National Science Foundation identified a handful of models that could predict the laboratory measurements within 10%.

While greater compliance with measurements is hoped for in the next decade,

10% accuracy is achievable now and should be considered a standard.

While a numerical model may be validated through comparisons with analytical results and laboratory data, this does not necessarily imply physical

realism for tsunamis of geophysical scales. One example is wave breaking.

While a numerical model may realistically approximate the solution of the

Navier-Stokes equations at laboratory scales, it may not do so at large scales.

Calculating the evolution of breaking waves involves calculating turbulent

shear terms and invoking turbulence closure constraints which are scale dependent. Therefore a reliability constraint needs to be applied, and this is

discussed in the next subsection.

For operational forecast models, propagation accuracy of 5% is now possible (Titov et al., 2005). For inundation models, accuracy of 10% with

respect to analytical results and laboratory data is also now possible. Both

should be considered as standards.

An associated accuracy constraint is grid resolution. This depends on

the complexity of the shoreline. On a fairly plane and very long beach such

as those of Southern California, a 100 m grid resolution may be sufficient.

The smallest offshore and onshore features likely to affect tsunami impact

on a coastal community should be reflected in the numerical grids. If a

community is fronted by a sand spit of width 100 m, at least four grid points

are needed to provide accurate resolution of the flow over the spit, implying

a resolution of less than 25 m. If a coastline is sparsely populated, a 100 m

resolution may be sufficient for satisfactory inundation maps, even when the

shoreline is complex.

3.3 Model reliability and realism constraints

Model reliability refers to how well a given model predicts inundation consistently and realistically from a geophysical point of view. Linear theory

may predict wave evolution consistently, but unfortunately not always in
12 Synolakis et al.

a realistic manner. Linear theory predicts that waves during shoaling will

keep growing in height, whereas in reality waves will eventually break, if

large enough and when exceeding threshold height-to-depth and height-towavelength ratios.

The issue of reliability is crucial. Several numerical models now widely

used include ad hoc friction factors. Given that these factors were developed

not to model the physical manifestation of frictional dissipation, but simply

to stabilize what is by its very nature a marginally stable computation,

it is not possible to know a priori how well a model that has been fairly

successful in a small number of cases performs in more general cases. For

example, a model developed and calibrated for stability with examples from

steep coastlines with less than 200 m inundation distances may not perform

equally well when employed to forecast the inundation from a tsunami that

penetrates more than 3 km inland, as the 26 December 2004 Boxing Day

tsunami did in Banda Aceh.

It is clear that any numerical model for inundation predictions needs to be

tested over a variety of scales from the laboratory to prototype to ensure both

reliability and realism. Ideally, inundation models should be continuously

tested with every new set of laboratory data or tsunami field data that

becomes available. This will also allow for their further improvement. On

the other hand, operational propagation forecasting is by its nature less

dependent on scale. Some of the existing methods to test for model accuracy,

reliability, and realism are discussed in Appendix A.

4. Summary of Model Evaluation Procedures

State-of-the-art inundation codes in use today have evolved through a painstaking process of careful validation and verification, while operational codes

have been developed through extensive verification with measurements from

real tsunamis, to the point that every new event poses a diminishing challenge. Mining this experience, procedures for approval and application of

numerical models for operational uses are proposed as follows:

1. Establishment of standards for model validation and verification;

2. Scientific evaluation of individual models;

3. Operational evaluation of individual models;

4. Development of operational applications for forecasting;

5. Procedures for transfer of technology to operations.

Only through parallel testing of models under identical conditions, as

when there is a tsunami emergency and an operational forecast is performed,

can the community determine the relative merits of different computational

formulations, an important step to further improvements in speed, accuracy,

and reliability.

Figure 1 illustrates the process of model development and evaluation for

operational use. The process of operational model development should be
NOAA tsunami numerical model standards 13

Developing NOAA

operational model

Candidate models 

Scientific community: 

Peer reviewed articles 

Conferences

Benchmarked models 

NSF workshops: 

Direct comparisons 

Published benchmarks

Operational evaluation 

NOAA testbed: 

NCTR, TWCs testing 

Operational standards

Operational models 

NOAA operations: 

NWS, TWCs use

Community modeling activity 

for

evaluated and benchmarked models 

Research & Development Applications

Refresh 

Next 

Generation 

Models 

Developing NOAA

operational model

Figure 1: Flowchart illustrating the process of NOAA operational model development.

as inclusive as possible, to ensure the active participation of the geophysical,

oceanographic, engineering, and NOAA research communities. Currently,

scientific evaluation of a tsunami model takes the form of peer-reviewed

publications and NSF tsunami model benchmark workshops. NOAA research and operations will evaluate the tsunami model operationally. The

Community Modeling Activity will solidify these evaluations into a sustainable developmental effort that will feed NOAA operations with new and

refreshed models. NOAA’s responsibility is to test benchmarked models for

operational suitability to define if a model fits the NOAA operational standards. This evaluation is done by operational and research components of

NOAA as a part of the Research-to-Application implementation plan. The

models—those that fit NOAA’s operational standards—may be implemented

for operations. The models that do not fit the standards (due to accuracy,

speed, or robustness deficiencies) will stay in the research community modeling activity for improvement and further development.

It is emphasized again that model testing must remain a continuous process. Operational products produced in real time during an actual event

must be thoroughly reviewed, and the operational models systematically

tested in hindcast mode after each tsunami strikes. The results must be

documented and reported to assist the community in developing and implementing improvements, through the identification and resolution of any

serious problems or inadequacies of the models and/or products.

While this process may appear onerous, it reflects our state of knowledge

as of December 2006, and is the only defensible methodology when human

lives are at stake.
14 Synolakis et al.

5. Acknowledgments

The authors thank Paul Whitmore (Director, West Coast and Alaska Tsunami Warning Center) for providing a draft of Section 2.6, “Operational

evaluation.” We also thank Paul and other reviewers of the manuscript for

fruitful discussions and many useful comments that shaped the document

into its current form. This work was partially funded by the U.S. National

Tsunami Hazard Mitigation Program. We thank the National Science Foundation of the United States for supporting some of the early validation exercises through benchmark testing in three individual workshops, and for

supporting the analytical studies and laboratory investigations that resulted

in the benchmark data sets discussed in this report.

6. References

Borerro, J., M. Ortiz, V.V. Titov, and C.E. Synolakis (1997): Field survey of

Mexican tsunami. Eos Trans. AGU, 78(8), 85, 87–88 (Eos Cover Article).

Briggs, M.J., C.E. Synolakis, G.S. Harkins, and D. Green (1995): Laboratory experiments of tsunami runup on a circular island. Pure Appl. Geophys., 144,

569–593.

Carrier, G.F., and H.P. Greenspan (1958): Water waves of finite amplitude on a

sloping beach. J. Fluid Mech., 17, 97–110.

Carrier, G.F. (1966): Gravity waves of water of variable depth. J. Fluid Mech., 24,

641–659.

Carrier, G.F., T.T. Wu, and H. Yeh (2003): Tsunami runup and drawdown on a

sloping beach. J. Fluid Mech., 475, 79–99.

Goring, D.G. (1978): Tsunamis—the propagation of long waves onto a shelf. W.M.

Keck Laboratory of Hydraulics and Water Resources, California Institute of

Technology, Report No. KH-R-38.

Grilli, S.T., I.A. Svenden, and R. Subrayama (1997): Breaking criterion and characteristics of solitary waves on a slope. J. Waterw. Port Coast. Ocean Eng.,

123(2), 102–112.

Hall, J.V., and J.W. Watts (1953): Laboratory investigation of the vertical rise of

solitary waves on impermeable slopes. Tech. Memo. 33, Beach Erosion Board,

U.S. Army Corps of Engineers, 14 pp.

Hammack, J.L. (1972): Tsunamis—A model for their generation and propagation.

W.M. Keck Laboratory of Hydraulics and Water Resources, California Institute

of Technology, Report No. KH-R-28.

Kˆano ̆glu, U. (1998): The runup of long waves around piecewise linear bathymetries.

Ph.D. Thesis, University of Southern California, Los Angeles, California, 90089-

2531, 273 pp.

Kˆano ̆glu, U. (2004): Nonlinear evolution and runup–rundown of long waves over a

sloping beach. J. Fluid Mech., 513, 363–372.

Kˆano ̆glu, U., and C.E. Synolakis (1998): Long wave runup on piecewise linear

topographies. J. Fluid Mech., 374, 1–28.

Kˆano ̆glu, U., and C. Synolakis (2006): Initial value problem solution of nonlinear shallow water-wave equations. Phys. Rev. Lett., 97, 148,501–148,504, doi:

10.1103/PhysRevLett.97.148501.

Keller, J.B., and H.B. Keller (1964): Water wave run-up on a beach. ONR Research

Report NONR-3828(00), Dept. of the Navy, Wash. DC, 40 pp.

Li, Y., and F. Raichlen (2000): Energy balance model for breaking solitary wave

runup. J. Waterw. Port Coast. Ocean Eng., 129(2), 47–49.
NOAA tsunami numerical model standards 15

Li, Y., and F. Raichlen (2001): Solitary wave runup on plane slopes. J. Waterw.

Port Coast. Ocean Eng., 127 (1), 33–44.

Li, Y., and F. Raichlen (2002): Non-breaking and breaking solitary run-up. J. Fluid

Mech., 456, 295–318.

Liu, P.L.-F., Y.-S. Cho, M.J. Briggs, U. Kˆano ̆glu, and C.E. Synolakis (1995): Runup

of solitary waves on a circular island. J. Fluid Mech., 320, 259–285.

Liu, P.L.-F., P. Lynett, and C.E. Synolakis (2003): Analytical solutions for forced

long waves on a sloping beach. J. Fluid Mech., 478, 101–109.

Liu, P.L.-F., C.E. Synolakis, and H.H. Yeh (1991): Report on the International

Workshop on Long Wave Runup. J. Fluid Mech., 229, 675–688.

Liu, P.L.-F., T.-R. Wu, F. Raichlen, C.E. Synolakis, and J. Borrero (2005): Runup

and rundown generated by three-dimensional sliding masses. J. Fluid Mech.,

536, 107–144.

Liu, P.L.-F., H. Yeh, and C. Synolakis (2007): Advanced Numerical Models for

Simulating Tsunami Waves and Runup. Advances in Coastal and Ocean Engineering, 10, 250 pp.

National Science and Technology Council (2005): Tsunami Risk Reduction for the

United States: A Framework for Action. A Joint Report of the Subcommittee

on Disaster Reduction and the United States Group on Earth Observations,

30 pp.

Pedersen, G., and B. Gjevik (1983): Runup of solitary waves. J. Fluid Mech., 135,

283–299.

Russel, J.S. (1845): Report on Waves. Rp. Meet. Brit. Assoc. Adv. Sci. 14th, 311–

390, John Murray, London.

Shuto, N. (1973): Shoaling and deformation of nonlinear waves. Coastal Eng. Japan,

16, 1–12.

Synolakis, C.E. (1986): The Runup of Long Waves. Ph.D. Thesis, California Institute of Technology, Pasadena, California, 91125, 228 pp.

Synolakis, C.E. (1987): The runup of solitary waves. J. Fluid Mech., 185, 523–545.

Synolakis, C.E., and J.E. Skjelbreia (1993): Evolution of maximum amplitude of

solitary waves on plane beaches. J. Waterw. Port Coast. Ocean Eng., 119(3)

323–342.

Synolakis, C.E., and E.N. Bernard (2006): Tsunami science before and after Boxing

Day 2004. Phil. Trans. R. Soc. A, 364(1845), doi: 10.1098/rsta.2006.1824,

2231–2265.

Tadepalli, S., and C.E. Synolakis (1994): The runup of N-waves on sloping beaches.

Proc. R. Soc. Lond. A. 445, 99–112.

Tadepalli, S., and C.E. Synolakis (1996): Model for the leading waves of tsunamis.

Phys. Rev. Lett., 77(10), 2141–2144.

Takahashi, T. (1996): Benchmark problem 4; the 1993 Okushiri tsunami—Data,

conditions and phenomena. In Long-Wave Runup Models, World Scientific,

384–403.

Titov, V.V., and C.E. Synolakis (1995): Modeling of breaking and nonbreaking

long-wave evolution and runup using VTCS-2. J. Waterw. Port Ocean Coast.

Eng., 121(6), 308–316.

Titov, V.V., and C.E. Synolakis (1997): Extreme inundation flows during the

Hokkaido–Nansei-Oki tsunami. Geophys. Res. Lett., 24(11), 1315–1318.

Titov, V.V., and F.I. Gonz ́alez (1997): Implementation and testing of the method of

splitting tsunami (MOST). NOAA Technical Memorandum ERL-PMEL-112,

PB98-122773, Pacific Marine Environmental Laboratory, Seattle, WA, 11 pp.

Titov, V.V., and C.E. Synolakis (1998): Numerical modeling of tidal wave runup.

J. Waterw. Port Ocean Coast. Eng., 124(4), 157–171.
16 Synolakis et al.

Titov, V.V., F.I. Gonz ́alez, E.N. Bernard, M.C. Eble, H.O. Mofjeld, J.C. Newman,

and A.J. Venturato (2005): Real-time tsunami forecasting: Challenges and

solutions. Nat. Haz., 35(1), 45–58.

Tuck, E.O., and L.S. Hwang (1972): Long wave generation on a sloping beach. J.

Fluid Mech., 51, 449–461.

Tanioka, Y., and F.I. Gonz ́alez (1998): The Aleutian earthquake of June 10, 1996

(Mw 7.9) ruptured parts of both the Andreanof and Delarof Segments. Geophys. Res. Lett., 25 (12), 2245–2248.

Whitmore, P.M. (2003): Tsunami amplitude prediction during events: A test based

on previous tsunamis. In Science of Tsunami Hazards, 21, 135–143.

Yeh H., P.L.-F. Liu, M. Briggs, and C.E. Synolakis (1994): Tsunami catastrophe

in Babi Island. Nature, 372, 6503–6508.

Yeh, H., P.L.-F. Liu, and C.E. Synolakis (1996): Long-Wave Runup Models. World

Scientific, 403 pp.

Yeh, H., V.V. Titov, V. Gusyiakov, E. Pelinovsky, V. Khramushin, and V. Kaistrenko

(1995): The 1994 Shikotan earthquake tsunami. Pure Appl. Geophys., 144(3/4),

569–593.

Zelt, J.A. (1991): The runup of breaking and nonbreaking solitary waves. Coastal

Eng., 125, 205–246.
NOAA tsunami numerical model standards 17

Appendix A: Existing Methods for Model

Validation and Verification

Benchmarking of numerical models can be classified into four categories:

Basic hydrodynamic considerations, analytical benchmarking, experimental

benchmarking, and field benchmarking. Here, specific benchmark problems

for validating and verifying computational tools for predicting the coastal

effect of tsunamis are described in detail. Some of the benchmark problems

described here were used as benchmark problems in the 1995 (Yeh et al.,

1996) and 2004 (Liu et al., 2007) Long-Wave Runup Models Workshops in

Friday Harbor, Washington, and Catalina Island, Los Angeles, California,

respectively.

1 Basic hydrodynamic considerations

Two most basic steps are required in ensuring that a numerical model works

for predicting evolution and inundation. While the first step is ensuring that

the model conserves mass, the second basic step is checking convergence of

the numerical code to a certain asymptotic limit.

1.1 Mass conservation

The conservation of mass equation is part of the equations of motion that

are solved in any numerical procedure, but cumulative numerical approximations can sometimes produce results that violate mass conservation. This is

particularly the case when friction factors are used, or smoothing to stabilize

inundation computations for breaking waves.

Conservation of mass can be checked by calculating the water volume at

the beginning and at the end of the computation. This should be done by

integrating disturbed water depth η(x, y, t) over the entire flow domain, i.e.,

if the flow domain extends from the maximum penetration during inundation

x = Xmax to the outer location of the source region XS, and y = Ymax to

YS, then total displaced volume V (t) is

V (t) =

X

S

Xmax



YS

Ymax

η(x, y, t)dxdy. (A1)

The integral of η(x, y, t) should be used instead of the integral of the

entire flow depth h(x, y, t) = η(x, y, t) + d(x, y, t)—where d(x, y, t) is the

undisturbed water depth—because the latter is likely to conceal errors in the

calculation. Typically, η d at offshore integrating h will simply produce

the entire volume of the flow domain and will mask errors. Note that testing

of the conservation of mass as above involves placing a closed domain within

reflective boundaries.

Numerical models use absorbing boundaries to allow the wave to radiate

outwards. While some loss of mass may in theory occur due to the differencing at the boundary, placing the computational boundaries far enough from

the source ensures this loss is minimal. Nonetheless, testing conservation
18 Synolakis et al.

of mass with open boundaries and monitoring the volume fluxes over the

corresponding boundaries at x = xb and y = yb is computationally possible.

The net flux needs to be added or subtracted to the total displaced volume.

The calculation of the entire volume of the wave motion of the flow domain at the beginning and end of the computation—to ensure that mass is

conserved—is very important in code validation. Once the associated mass

fluxes at the boundaries are considered, numerical errors can be highly additive, and mass invariably decreases in numerical computations. It is essential

that codes self-monitor the entire volume; if the difference between start and

finish is greater than 5%, the grid needs to be readjusted. Usually with few

changes in Δx, Δy, and Δt, the conservation of mass can be improved.

Calculations of conservation of mass should be such that the total initial

displaced volume V (t = 0) should be within 5% of the total displaced volume

at the end of the computation V (t = T) where T represents the computation

end time. It is assumed that the end of the computation is when the initial

wave is entirely reflected and offshore.

1.2 Convergence

The next basic step is checking convergence of the numerical code to a

certain asymptotic limit, presumably the actual solution of the equations

solved. The grid steps Δx and Δy need to be halved, and the time step Δt

reduced appropriately to conform with the Courant–Friedrics–Levy (CFL)

criterion. The optimal locations to check convergence are the extreme runup

and rundown. A graph needs to be prepared presenting the variation of the

calculated runup and rundown (ordinate) with the step size (abscissa). As

the step size is reduced, the numerical predictions should be seen to converge

to a certain value, and further reductions in step size should not change the

results.

2 Analytical benchmarking

The shallow water-wave (SW) equations describe the evolution of the watersurface elevation and of the depth-averaged water particle velocity of waves

with wavelengths large compared with the depth of propagation. The equations assume that the pressure distribution is hydrostatic everywhere, i.e.,

there is no variation with depth of any of the other flow variables. In this

section we present several analytic solutions to the 1+1 version of these

equations. As stated in Section 2 of this report, 1+1 models are generally

unreliable for inundation mapping and entirely inadequate for operational

tsunami forecasting, but they are invaluable to the process of testing and

validating models.

2.1 Single wave on a simple beach

The so-called canonical problem of the shallow water-wave equations is the

calculation of a long wave climbing up a sloping beach of angle β coupled to

a constant-depth region (Fig. A1). The origin of the coordinate system is at
NOAA tsunami numerical model standards 19

initial shoreline

position x = 0

t = tmax Y

X

x = Xs x = X0

L

H

t = 0

β

R

d

~

Figure A1: Definition sketch for canonical bathymetry, i.e., sloping beach connected to a constant-depth region.

the initial position of the shoreline and x increases seaward. Dimensionless

variables are introduced as follows:

x = x ̃

d



, (η, h0) = ( ̃η, h ̃0)

d

 , u = u ̃ 

gd



, and t = t

 ̃ 

d/ g

. (A2)

Here quantities with tilde are dimensional and η is the amplitude, u is the

depth-averaged horizontal velocity, h0 is the undisturbed water depth, d

 is

the depth of the constant-depth region, and g is the gravitational acceleration. The topography is described by h0(x) = x tan β when x ≤ X0 = cot β

and h0(x) = 1 when x ≥ X0 = cot β. Even though dimensionless variables

are not preferred in engineering practice, here they have distinct advantages

as everything scales simply with an offshore characteristic depth. In numerical solutions, dimensional variables are most often used.

Consider a tsunami evolution problem described by the 1 + 1 nonlinear

shallow water-wave (NSW) equations:

ht + (uh)x = 0,

ut + uux + ηx = 0, (A3)

with h(x, t) = η(x, t)+h0(x). Through elementary manipulations, neglecting

nonlinear terms, (A3) reduces to

ηtt − (ηxh0)x = 0, (A4)

an equation known as the linearized shallow water-wave (LSW) equation.

Over constant-depth h0 = 1, then

ηtt − ηxx = 0, (A5)

the classic one-dimensional wave equation.

The solution follows directly from the Fourier transform of the equation

when a boundary condition for the wave amplitude η(X1, t) is specified, i.e.,
20 Synolakis et al.

when the incident wave η(x, t) at the constant-depth region is known at some

x = X1, and can be described by a Fourier integral of the form

η(X1, t) = +∞

−∞

Φ(ω)e−iωtdω. (A6)

Consider the canonical problem of a constant-depth region of depth h0 = 1

joined to a sloping beach of angle β, with the toe of the beach at x = X0.

Synolakis (1986, 1987) matched the linear theory solution at the constantdepth with the solution over the sloping beach as derived by Keller and

Keller (1964) to derive the solution over the sloping beach for a wave with

a transform as given by (A6),

η(x, t)=2 +∞

−∞

Φ(ω)

J0(2ω

√xX0)e−iω(X0+t)

J0(2X0ω) − iJ1(2X0ω)

dω, (A7)

where X0 = cot β. This solution is only valid when 0 ≤ x ≤ X0; when x < 0,

(A4) does not reduce to Bessel’s equation. Notice that the integral (A7) can

be evaluated with standard numerical methods; however, the advantage of

this form is that it allows calculation of the solution for many physically

realistic tsunami waveforms simply by plugging in the Φ(ω) of the incoming

wave, hopefully known at some offshore location X1.

2.1.1 Solitary wave evolution and runup As discussed in Synolakis (1986, 1987) it is possible to derive exact results for the evolution and

runup of solitary waves based on linear theory. Solitary waves have long

been used as a model for the leading wave of tsunamis. Solitary waves were

first described by Russel (1845) as the great waves of translation, and consist of a single elevation wave. While capturing some of the basic physics

of tsunamis, they do not model the physical manifestation of tsunamis in

nature, which are invariably N-wave like with a leading-depression wave followed by an elevation wave. A solitary wave centered offshore at x = Xs at

t = 0 has the following surface profile,

η(x, 0) = Hsech2γ(x − Xs), (A8)

where γ = 3H/4 and H is the dimensionless wave height, i.e., H = H/  d

.

The function Φ(ω) associated with this profile is derived in Synolakis (1986)

and it is given by:

Φ(ω) = 2

3

ωcosech(αω)eiωXs , (A9)

where α = π/(2γ).

In the context of water-wave theory, the solitary wave (A8) is an exact solution of the Korteweg-de-Vries (KdV) equation; therefore, a KdV

solitary wave propagates over constant-depth without any change in shape.

The KdV theory is both dispersive and nonlinear, and solitary waves are

the only waves with this unique property of unchanging shape. However,

(A8) can be used as an initial condition for other wave theories, without,

of course, a priori expectation that the SW model will preserve the classic
NOAA tsunami numerical model standards 21

soliton properties, which include their ability to go through each other (interact) without any change in shape through nonlinear interactions. This

having been said, since the LSW is nondispersive and linear, hence all waves

propagate over constant-depth without any change in shape. However, in

the range of wave steepness and amplitudes relevant for tsunamis, it is now

well established that, at least for the 1 + 1 problem far from the shoreline,

the LSW theory, which also preserves the wave shape for propagation over

constant-depth, is quite adequate (Liu et al., 1991) and useful when the

engineering problem has simple geometry.

The derivation of the amplitude evolution for solitary waves is not as

straightforward as often assumed. When superposing sinusoids, there is a

frequency ω-dependent phase shift. It is therefore not obvious that linear superposition will produce a similar amplitude variation given this frequencydependent phase shift.

To describe the evolution of a solitary wave up a plane beach, Synolakis

(1986, 1987), substituted (A9) into (A7) to obtain

η(x, t) = 4

3

+∞

−∞

ωcosech(αω)

J0(2ω

√xX0)e−iω(X0−Xs+t)

J0(2X0ω) − iJ1(2X0ω) dω, (A10)

where, as earlier, α = π/√

3H. This integral can be evaluated directly

through contour integration. In the region where the wave evolves on the

sloping beach far off the shoreline, x is large, and (A10) becomes

η(x, t) = 4π2

3α2 (

X0

x )

1/4

+∞

n=1

(−1)n+1ne−(π/a)θ

n, (A11)

with θ = X0 − Xs − t − 2

√xX0. The maximum of the power series is

1/4, therefore the maximum local value of the wave amplitude ηmax is given

explicitly by

ηmax

H = (X0x )

1/4 = ( 1

h0

)

1/4; (A12)

this is an amplitude variation similar to Green’s law.

The region over which (A12) applies is the region of gradual shoaling; the

region of rapid shoaling is often identified with the Boussinesq result, i.e.,

ηmax ∼ h. The fact that both evolution laws may coexist was first identified

by Shuto (1973). Synolakis and Skjelbreia (1993) also present results which

show that Green’s law type evolution is valid over a wide range of slopes and

for finite-amplitude waves at least in the region of gradual shoaling.

The results (A10) can now be readily applied to derive a result for the

maximum runup of a solitary wave climbing up a sloping beach. Writing

R(t) = η(0, t), i.e., R(t) is the free-surface elevation at the initial shoreline;

in the LSW theory, the shoreline does not move beyond x = 0. The maximum value of R(t) is the maximum runup R, arguably the most important

parameter in the long-wave runup problem, and it is the maximum vertical

excursion of the shoreline. Per Synolakis (1986), from equation (A10), it can
22 Synolakis et al.

be deduced that

R(t)=8H 

+∞

n=1

(−1)n+1ne−2γ(Xs−X0−t)n

I0(4γX0n) + I1(4γX0n)

. (A13)

The series can be simplified further by using the asymptotic form for large arguments of the modified Bessel functions. The resulting series is of the form

s = ∞

n=1(−1)n+1n3/2χn; its maximum value smax occurs at χ = 0.481 =

e−0.732 with smax = 0.15173. Then the following expression results for the

maximum runup R:

R = 2.831cot βH5/4. (A14)

This result is formally correct when √

H  0.288 tan β—the assumption

implied when using the asymptotic form of the Bessel functions. Equation

(A14) was first derived by Synolakis (1986) and has since been referred

to as the runup law and shown in Fig. A2. As will be apparent in later

sections, this methodology is quite powerful to find the maximum runup and

it allows calculation of the runup of other waveforms such as N-waves, not

to mention the runup of waves evolving over piecewise-linear bathymetries.

Recent results suggest that the dependence of the runup on the slope and on

the offshore wave height in a two-dimensional problem of idealized conditions

is often quite similar to this one-dimensional power law.

The asymptotic result (A14) is valid for waves that do not break during

runup, suggesting that it is appropriate to use the qualifier nonbreaking for

waves that do not break during runup but may or may not break during rundown. The real usefulness of any asymptotic result is how well it identifies

the scaling, i.e., it can identify the solution dependence on the problem parameters; numerical solutions will invariably produce more accurate specific

predictions, but they can rarely provide useful information about the problem scaling. To check if the runup law (A14) provides the correct scaling,

Synolakis (1986, 1987) examined the classic laboratory data set of Hall and

Watts (1953), Fig. A2. That study includes both breaking and nonbreaking wave data without identifying them as such, clearly because there was

no realization of the differences; the empirical runup relationships derived

by Hall and Watts (1953) are not directly applicable when determining the

runup of nonbreaking waves. To perform a posteriori identification of those

data, the breaking criterion H < 0.479(cot β)−10/9 was used.

2.1.2 N-wave runup Most tsunami eyewitness accounts suggest that

tsunamis are N-wave like, i.e., they are dipolar, which means they appear

as a combination of a depression and an elevation wave, and frequently as

a series of N-waves, sometimes known as double N-waves. Up until the late

1990s, the solitary wave model was used exclusively to evaluate the runup

of tsunamis. The N-wave model was motivated by observations from earthquakes in Nicaragua [1 September 1993], Flores, Indonesia [12 December

1992], Okushiri, Japan [7 July 1993], East Java, Indonesia [6 June 1994],

Kuril Islands, Russia [4 October 1994], Mindoro, Philippines [14 November 1994], Manzanillo, Mexico [9 October 1995], Chimbote, Peru [3 March

1996], Papua New Guinea [17 July 1998], Vanuatu [26 November 1999], and
NOAA tsunami numerical model standards 23

1.0 10

10

1

0.1

0.01

1 : 19.85

1 : 11.43

1 : 5.67

1 : 3.73

1 : 2.75

1 : 2.14

1 : 1.00

Runup law

0.01 0.1

cot β(H/d)5/4 2.831

R/d

Figure A2: Laboratory data for maximum runup of nonbreaking waves climbing

up different beach slopes: 1:19.85 (Synolakis, 1986), ♦ 1:11.43 (Hall and Watts,

1953),  1:5.67 (Hall and Watts, 1953), ∗ 1:3.73 (Hall and Watts, 1953), 	 1:2.75

(Pedersen and Gjevik, 1983), + 1:2.14 (Hall and Watts, 1953), × 1:1.00 (Hall and

Watts, 1953). Solid line represents the runup law (A14).

Caminade, Peru [20 June 2001], all of which produced tsunami waves which

caused nearby shorelines to first recede before advancing. The most specific

description was during the 9 October 1995 Manzanillo, Mexico earthquake.

One eyewitness saw the shoreline retreat beyond a rock outcrop which was

normally submerged in over 5 m depth and at a distance of about 400 m

from the shoreline, suggesting a leading-depression wave. Before the megatsunami of 26 December 2004, this was the only photographic evidence of

these leading-depression waves (Borrero et al., 1997). Recall that the megatsunami manifested itself first with a rapid withdrawal of the water surface

in most locales east of the rupture zone.

Modeling tsunamis with solitary waves cannot possibly explain these observations, because a solitary wave is technically a leading-elevation wave.

Therefore, and to reflect the fact that tsunamigenic faulting in subduction

zones associated with both vertical uplift and subsidence of the sea bottom,

Tadepalli and Synolakis (1994) conjectured that all tsunami waves at generation have an N-wave or dipole shape. Tadepalli and Synolakis (1994, 1996)

proposed a general function as a unified model for both near-shore and farfield tsunamis as generalized N-waves. When a wave propagates with the

trough first it is referred to as a leading depression N-wave or LDN. When

the crest arrives first, it is a leading-elevation wave or LEN. Also, Tade-
24 Synolakis et al.

palli and Synolakis (1994) defined another type of N-wave of this class with

leading-elevation and depression waves of the same height and at a constant

separation distance and refer to this wave as an isosceles N-wave with a

surface profile given by

η(x, 0) = 3

√3H

2

sech2[γ(x − XN )]tanh[γ(x − XN )], (A15)

with γ = 3

2

3

4H. Wave profile (A15) is an LDN and has a maximum wave

amplitude H. The function Φ(ω) associated with this profile is derived in

Tadepalli and Synolakis (1994) using contour integration and it is given by:

Φ(ω) = 2i

3γ

ω2cosech(πω

2γ

)eiωXN . (A16)

Now, the maximum of R(t) = η(0, t) can be evaluated for LEN using the

symmetry of the profile given in (A15) and one can find that

RN−wave = 3.86cot βH5/4. (A17)

Comparing the runup of the Boussinesq solitary wave (A14) with the runup

of an isosceles N-wave, RN−wave = 1.364RSolitary. Because of the symmetry

of the profile, this is also the minimum rundown of an isosceles leadingdepression N-wave. Tadepalli and Synolakis (1994) showed that the normalized maximum runup of nonbreaking isosceles LEN is smaller than the runup

of isosceles LDN, and that both are higher than the runup of a solitary wave

with the same wave height, and the latter is known as the N-wave effect

(Fig. A3).

The two-dimensional character of the generation region limits the direct

application of the N-wave and solitary wave models. However, N-wave theory

does provide a conceptual framework for analysis and for explaining certain

field observations qualitatively.

2.1.3 Boundary value problem To solve the nonlinear set (A3) for

the sloping beach case, h0(x) = x tan β, Carrier and Greenspan (1958) introduced the hodograph transformation known as Carrier–Greenspan transformation:

u = ψσ

σ , (A18)

η = ψλ

4 − u22 , (A19)

t = cot β(

ψσ

σ − λ2

), (A20)

x = cot β(

σ2

16 − ψλ4 +

u2

2 ), (A21)

into (A3), and they derived the following second-order ordinary differential

equation:

(σψσ)σ = σψλλ. (A22)
NOAA tsunami numerical model standards 25

10–1

10–2

10–3

10–3 10–2 10–1

0.01

η 0

–0.01 0 50

x

100

0.01

η 0

–0.010 50

x

100

R/d

cot β(H/d)5/4 3.86

Figure A3: Maximum runup of isosceles N-waves and solitary wave. Top and

lower set of points are results for the maximum runup of leading-depression and

-elevation isosceles N-waves, respectively. Dotted line represents the runup of solitary wave (A14). Lower and upper insets compare a solitary wave profile to a

leading-depression and -elevation isosceles N-waves, respectively.

Note that the original Carrier–Greenspan transformation does not include

cot β in (A20) and (A21) because of a different scaling as given in (A25) in

section A2.1.4. Note also the similarity with the linear form of the shallow

water-wave equation (LSW), (ηxh0)x = ηtt. Also, notice the conservation

of difficulty. Instead of having to solve the coupled nonlinear set (A3), one

now has to solve a linear equation, but the transformation equations which

relate the transformed variables with the physical variables are nonlinear,

coupled, and implicit. Yet, a redeeming feature is that in the hodograph

plane, i.e., in the (σ, λ)-space, the shoreline is always at σ = 0. This allows for direct analytical solutions without the complications of the moving

shoreline boundary.

In general, it is quite difficult to specify initial or boundary data for the

nonlinear problem in the physical (x, t)-space coordinates without making

restrictive assumptions; a boundary condition requires specification of the

solution at (X0, ∀t), and an initial condition specification at (t0, ∀x), but, in
26 Synolakis et al.

practice, the wave approaching the beach is only known offshore for (X0 ≥

cot β, t < t0), where t0 is the time that the wave reaches the x-location X0.

Even when boundary or initial conditions are available in the (x, t)-space,

the process of deriving the equivalent conditions in the (σ, λ)-space is not

trivial.

These difficulties have restricted the use of this transformation to problems that can be reduced directly to those solved by Carrier and Greenspan (1958). Synolakis (1986) revived the Carrier–Greenspan formalism by

developing a method to specify a boundary condition including reflection.

Synolakis (1986) used the solution of the equivalent linear problem, as given

by (A7), at the seaward boundary of the beach, i.e., at x = X0 = cot β

corresponding to σ = σ0 = 4. Then, equation (A19) implies that η(X0, t) ≈

1

4ψλ(4, λ). Assuming that ψ(σ0, λ) → 0 as λ → ±∞, Synolakis (1986, 1987)

showed that the Carrier–Greenspan potential is given by

ψ(σ, λ) = −16i

X0

+∞

−∞

Φ(κ)

κ

J0(σκX0/2)e−iκX0(1− λ

2 )

J0(2κX0) − iJ1(2κX0) dκ. (A23)

Even though the solution now can be obtained in the (σ, λ)-space using

(A18) and (A19) and can be converted to the solution in the (x, t)-space

through (A20) and (A21), the problem with this transformation is deriving

a solution for a particular-time t∗ or at a particular-location x∗. Synolakis

(1986) and later Kˆano ̆glu (2004) evaluated the solution either for given t∗ or

at given x∗ using the Newton–Raphson iteration algorithms, respectively;

λi+1 = λi −

	

t

∗ − t(λ)

∂

∂λ [t∗ − t(λ)]
















λi

or σi+1 = σi −

	

x∗ − x(σ)

∂

∂σ [x∗ − x(σ)]
















σi

.

(A24a, b)

Here t(λ) and x(σ) are given with (A20) and (A21), respectively.

An astonishing feature of the NSW is that the predictions for the maximum runup are identical to those of the LSW, when identical boundary

conditions are specified at X0 = cot β. The maximum runup according to

LSW is the maximum value attained by the wave amplitude at the initial position of the shoreline, while the maximum runup is given by the maximum

value of the amplitude at the evolving shoreline η(xs, λ), where xs is the

x-coordinate of the shoreline tip and corresponds to σ = 0. Carrier (1966)—

without reflection—and Synolakis (1987)—including reflection—have shown

that the linear and nonlinear theory produce mathematically identical predictions. Nonlinear analytical evolution of solitary wave with H = 0.0185

over a sloping beach of 1:19.85 is given in Fig. A4.

2.1.4 Initial value problem The nonlinear evolution of a wave over

a sloping beach is theoretically and numerically challenging due to the moving boundary singularity. Yet, it is important to have a good estimate of

the shoreline velocity and associated runup–rundown motion, since they are

crucial for the planning of coastal flooding and of coastal structures. As

explained in the previous section, Synolakis (1987) solved this problem as a

boundary value problem considering canonical bathymetry. Kˆano ̆glu (2004)
NOAA tsunami numerical model standards 27

t = 25

t = 35

t = 45

t = 55

t = 65

x

0.04

0.02

0

–0.02

0

0.04

0.02

–0.02

0.08

0.06

0.04

0.02

0

–0.02

0.08

0.06

0.04

0.02

0

–0.02

0.10

0.06

0.04

0.02

0

–0.02

–0.04

–2 0 2 4 6 8 10 12 14 16 18 20

–2 0 2 4 6 8 10 12 14 16 18 20

–2 0 2 4 6 8 10 12 14 16 18 20

–2 0 2 4 6 8 10 12 14 16 18 20

–2 0 2 4 6 8 10 12 14 16 18 20

η

Figure A4: Time evolution of H = 0.0185 initial wave over a sloping beach with cot β = 19.85 from t = 25

to 65 with 10 increments. Constant depth-segment starts at X0 = 19.85. While markers show experimental

results of Synolakis (1986, 1987), solid lines show nonlinear analytical solution of Synolakis (1986, 1987).
28 Synolakis et al.

η(x,t) 

h(x) 

X 

ηs(λ) 

β

Y 

Figure A5: Definition sketch for the initial value problem.

solved nonlinear evolution of any given wave form over a sloping beach as an

initial value problem (Fig. A5). It is proposed that any initial waveform can

first be represented in the transform space using the linearized form of the

Carrier–Greenspan transformation for the spatial variable, then the nonlinear evolutions of these initial waveforms can be directly evaluated. Later,

Kˆano ̆glu and Synolakis (2006) solved the similar problem considering a more

general initial condition, i.e., initial wave with velocity.

Kˆano ̆glu (2004) considers NSW equations (A3) with slightly different

nondimensionalization than (A2), i.e., using the reference length l as a scaling parameter, the dimensionless variables are introduced as

x = x

l

, (h, η) = (

h, η)

ltan β

, u = u 

g ltan β

, and t = t



l/(g tan β)

. (A25)

Using the original Carrier–Greenspan transformation—without cot β coefficient in (A20) and (A21)—it is possible to reduce the NSW equations to

the following single second-order linear equation the same as (A22):

σφλλ − (σφσ)σ = 0, (A26)

using the Riemann invariants of the hyperbolic system (Carrier and Greenspan, 1958). The Carrier–Greenspan transformation not only reduces the

nonlinear shallow water-wave equations into a second-order linear equation,

but also fixes the instantaneous shoreline to σ = 0 in the (σ, λ)-space as

explained previously. Furthermore, a bounded solution at the shoreline combined with a given initial condition in terms of a wave profile at λ = 0 in the

(σ, λ)-space, η(σ, 0) implies the following solution in the transform space,

φ(σ, λ) = −

∞

0

∞

0

1

ω

ξ2Φ(ξ)J0(ωσ)J1(ωξ) sin(ωλ)dωdξ, (A27)

where Φ(σ) = uλ(σ, 0) = 4ησ(σ, 0)/σ. Further, given the initial waveform
NOAA tsunami numerical model standards 29

η(σ, 0), the evolution of the water-surface elevation is now given by

η(σ, λ) = 1

4

φλ − 1

2

u2

= −1

4

∞

0

ξ2Φ(ξ)

∞

0

J0(ωσ)J1(ωξ) cos(ωλ)dω




dξ



− 1

2

∞

0

ξ2Φ(ξ)

∞

0

J1(ωσ)

σ

J1(ωξ) sin(ωλ)dω




dξ

2

,(A28)

where, again, Φ(σ)=4 ησ(σ, 0)/σ.

Since it is important for coastal planning, simple expressions for shoreline

runup–rundown motion and velocity are useful. Considering the shoreline

corresponds to σ = 0 in the (σ, λ)-space, (A28) reduces to the following

equation for the runup–rundown motion:

ηs(λ) = η(0, λ) = 1

4

φλ − 1

2

u2

s

= −1

4

∞

0

ξ2Φ(ξ)

∞

0

J1(ωξ) cos(ωλ)dω




dξ



− 1

2

∞

0

ξ2Φ(ξ)

∞

0

1

2

ωJ1(ωξ) sin(ωλ)dω




dξ

2

. (A29)

Here us and ηs represent shoreline velocity and motion, respectively. The

singularity of the u = φσ/σ at the shoreline (σ = 0) is removed with the

consideration of the limσ→0 [J1(ωσ)/σ] = 1

2ω.

The difficulty of deriving an initial condition in the (σ, λ)-space is overcome by simply using the linearized form of the hodograph transformation

for a spatial variable in the definition of initial condition. It is proposed that

any initial waveform can first be represented in the transform space using

the linearized form of the Carrier–Greenspan transformation for the spatial variable ((A21) without cot β coefficient), then the nonlinear evolutions

of these initial waveforms can be directly evaluated. Once an initial value

problem is specified in the (x, t)-space as η(x, 0), the linearized hodograph

transformation x ∼= 1

16σ2 is used directly to define the initial waveform in the

(σ, λ)-space, η( 1

16σ2, 0). Thus Φ(σ)=4ησ( 116σ2, 0)/σ is found, and φ(σ, λ)

follows directly through a simple integration, as in (A28). Then it becomes

possible to investigate any realistic initial waveform such as Gaussian and

N-wave shapes employed in Carrier et al. (2003) and the isosceles and general N-waves defined by Tadepalli and Synolakis (1994). Again, solution in

the physical space can be found using the Newton–Raphson algorithm proposed by Synolakis (1987) and later used by Kˆano ̆glu (2004), as presented

in (A24a, b).

The shoreline runup-rundown motion and velocity are presented for one

of the initial profiles given by Carrier et al. (2003):

η(x, 0) = H1 exp(−c1(x − x1)

2) − H2 exp(−c2(x − x2)2). (A30)

The following initial profile can be obtained in the transform space after

using the linearized form of the transformation for the spatial variable:

η(σ, 0) ≈ H1 exp(− 1

256

c1(σ2 − σ2

1)

2) − H2 exp(− 1

256

c2(σ2 − σ2

2)

2), (A31)
30 Synolakis et al.

0 2 4 6 8

−0.02

−0.01

0

0.01

(a)

x

η

0 1 2 3 4 5

–0.06

–0.03

0

0.03

0.06

(b)

t

ηs

0 1 2 3 4 5

−0.3

−0.15

0

0.15

0.3

(c)

t

u

s

Figure A6: (a) The leading-depression initial waveform (A30) presented by Carrier et al. (2003) with H1

= 0.006, c1 = 0.4444, x1 = 4.1209, H2 = 0.018, c2 = 4.0, and x2 = 1.6384 (solid line) compared with the one

resulting from approximation (evaluated through (A28)) (dots), (b) shoreline wave height, and (c) shoreline

velocity.

which leads to the definition of the initial condition Φ(σ):

Φ(σ) = − 1

16

H1 c1 (σ2 − σ2

1) exp(− 1

256

c1(σ2 − σ2

1)

2)

+

1

16

H2 c2 (σ2 − σ2

2) exp(− 1

256

c2(σ2 − σ2

2)

2). (A32)

Figure A6a compares the initial waveforms defined in the physical space

as in (A30) with the one resulting from the approximation of it, i.e., (calculated through (A28)). The linearized form of the spatial variable in the

definition of the initial waveforms gives satisfactory comparison. Figures

A6b and A6c present the shoreline runup–rundown motion and velocity, respectively, calculated from equation (A29) using the corresponding parts.

It should be added that the solution presented here cannot be evaluated

when the Jacobian of the transformation, J = xσ tλ − xλ tσ, breaks down.

Even though the transformation might become singular at certain points,

the solution can still be obtained at other points, since local integration can

be performed without prior knowledge of the dependent variables, unlike in

numerical methods. This feature is discussed in detail in Synolakis (1987)

and Carrier et al. (2003), and is not explained further in here.

2.2 Solitary wave on composite beach

Most topographies of engineering interest can be approximated by piecewiselinear segments allowing the use of LSW to determine approximate analytical

results for the wave runup of more complicated waveforms, in closed form.

In principle, fairly complex bathymetry can be represented through a combination of positive and/or negative segments and constant-depth segments.

Solutions of the LSW (A4) at each segment can be matched analytically at

the transition points between the segments, and then the overall amplification factor and reflected waves can be determined, analytically. We consider

here an application to check this assertion with a complex topography consisting of three segments and a vertical wall (Fig. A7). Laboratory data
NOAA tsunami numerical model standards 31

VERTICAL

WALL

WAVE

MAKER

15 m 4.4 m 2.9 m 0.9 m

23.2 m

d = 21.8 cm

d = 18.8 cm

1/13

1/150

1/53

 WAVE

 GAGES

1 2 3 4 5 6 7 8 9 10

~

~

Figure A7: Definition sketch for Revere Beach. Not to scale. Transition points

and segments are numbered from shoreward to seaward, i.e., segment numbers 1, 2,

3, and 4 show 1/13, 1/150, 1/53, and constant-depth segments, transition points 1,

2, and 3 show transition points between 1/13 and 1/150, 1/150 and 1/53, and 1/53

and constant-depth segment, respectively.

exist for this topography from a U.S. Army Corps of Engineers, Coastal Engineering Research Center, Vicksburg, Mississippi experiment of wave runup

on a model of Revere Beach, Massachusetts. This beach profile and the laboratory data are discussed in greater detail in Yeh et al. (1996), Kˆano ̆glu

(1998), and Kˆano ̆glu and Synolakis (1998). The profile of the model consists

of three piecewise-linear slopes of 1:13 (1st segment), 1:150 (2nd segment),

and 1:53 (3rd segment) from shoreward to seaward. At the shoreline there is

a vertical wall. In the laboratory experiments, to evaluate the overtopping

of the seawall, the wavemaker was located 23.22 m and tests were done at

two depths, at 18.8 cm and at 21.8 cm.

Kˆano ̆glu (1998) and Kˆano ̆glu and Synolakis (1998) associated constantdepth segment of depth hr with the matrix:

Dpr =

⎡

⎣

e

− iωxp √hr e

iωxp √hr

ie− iωx

√ p

hr −ie

iωx

√ p

hr

⎤

⎦ , (A33)

and the linearly varying depth region with positive slope mr is associated

with the matrix:

Qpr =

	 J0(2ω

hp/mr) Y0(2ωhp/mr)

J1(2ω

hp/mr) Y1(2ωhp/mr)



. (A34)

Above, in equations (A33) and (A34), the first subscript p identifies the

transition point, and the second subscript r identifies the segment, i.e., if a

segment has two transition points, there are two associated 2 × 2 matrices.

Here nondimensional quantities are defined as in (A2) using the depth of the

4th segment h4 as the characteristic length scale.

Using (A33) and (A34), continuities of water surface elevation and its

spatial derivative at the transition points lead to the following matrix equa-
32 Synolakis et al.

tions:

Q11V1 = Q12V2, Q22V2 = Q23V3, and Q33V3 = D34V4, (A35)

for each transition point. Here column vectors Vr = [Vr1, Vr2]

T identify

unknowns for each segment. One additional condition must be considered at the vertical wall, i.e., perfect reflection. Boundary condition at

the wall requires V1 to be a column vector considering perfect reflection,

i.e., Qw = 

1, −J1(2ω

√hw/m1)/Y1(2ω√hw/m1)

T . V1 = QwV11 with hw is

the dimensionless water depth at the wall. Combining (A35) with perfect

reflection gives

QwV11 = Q−1

11 Q12Q−122 Q23Q−133 D34V4. (A36)

To evaluate time histories of the surface elevation, the following integral

must be evaluated

ηr(x, t) =

+

∞

−∞

Φ(ω) 1

V41

⎧

⎨

⎩

Vr1J0(

2ω

√hr(x)

mr ) + Vr2Y0(

2ω

√hr(x)

mr )

V41e

− √iωx

h4 + V42e

√iωx

h4

⎫

⎬

⎭

e−iωtdω,

(A37)

where r = 1, 2, 3 and V41 and V42 represents incident and reflected wave

amplitudes, respectively. The temporal variation of the water surface elevation can be calculated through (A37) evaluating the unknowns in terms of

V41 through (A36) and (A35).

More interestingly, following analysis of Synolakis (1986) the maximum

runup of solitary waves over Revere Beach can be calculated analytically

given Φ(ω) as in (A9). The amplitude at the shoreline is

R(t) = η(0, t) = −(4/3) m1

π

√hw

+∞

−∞

cosech(αω)eiω(Xs−X3−t)

φ(ω) + iψ(ω) dω. (A38)

Here φ(ω) + iψ(ω) is a complicated but determinable quite easily from expression (A35) and (A36) in terms of Bessel functions of zero and first order, using symbolic mathematical tools. Kˆano ̆glu (1998) conjectured that

φ(z) + iψ(z) is an entire function in the upper half plane, and derived the

Laurent series expansion and its asymptotic form as

R(t)=8h−1/4 w H 

+∞

n=1

(−1)n+1nX n, (A39)

where X = e

− π

α (Xs−X3−t+2{

√h1−√hw

m1 +

√h2−√h1

m2 +

√h3−√h2

m3 })

. This is a power

series of the form (−1)n+1nX n and its maximum is equal to 1/4. Therefore

the maximum runup for solitary waves propagating up Revere Beach is given

by the runup law,

R = 2h−1/4 w H. (A40)
NOAA tsunami numerical model standards 33

10

1

0.1

0.01

0.001

0.001 0.01 0.1 1 10

R

H

Figure A8: Comparison of maximum runup values for the linear analytical solution

(A39) and laboratory results for two different depths, i.e., d

 = 18.8 cm and d =

21.8 cm.

The runup law above suggests that the maximum runup only depends on

the depth at the seawall fronting the beach, and it does not depend on any

of the three slopes in front of the seawall. The runup law (A40) predicts the

nonbreaking data surprisingly well (Fig. A8).

Note that no Jacobian regularization conditions as yet exist for wave

evolution on composite beaches, and this theory can only be applied with

caution. For lack of a better method, it is recommended that a limiting

height be determined for applying this methodology on a composite beach,

by finding the limiting height for the single beach with slope equal to the

least steep slope in the composite topography under consideration. This

practice will produce adequate validity thresholds for bathymetries such as

Revere Beach, but it should not be used when abrupt changes in bathymetry

exist. In most cases, for the range of wave-steepness and beach slopes relevant in tsunami studies, this method will give good results with fairly small

computational effort.

2.3 Subaerial landslide on simple beach

Liu et al. (2003) considered tsunami generation by a moving slide on a uniformly sloping beach, using the forced linear shallow water-wave equation as

in Tuck and Hwang (1972);

∂2η

∂t2 − tanβμ

∂

∂x(x

∂η

∂x) = ∂2h0∂t2 , (A41)
34 Synolakis et al.

in the dimensionless form. Here z = η(x, t) is the free surface elevation and

z = −h(x, t) the sea floor. On a sloping beach h(x, t) = H(x)−h0(x, t) with

H(x) = x tan β/μ; tan β is the beach slope, h0(x, t) is the time-dependent

perturbation of the sea floor with respect to the uniformly sloping beach,

δ

 and L are the maximum vertical thickness of the sliding mass and its

horizontal length, respectively, and μ = δ/ L. Normalization is done as

(η, h, H, h0) = (η,  

h, H,  h0)

δ

 , x = x

L , and t = t



δ/ g/μ 

. (A42)

Here again quantities with tilde represent dimensional quantities. The focus is on thin slides where μ = δ/ L 1. Since the shallow water-wave

assumption requires that tan β 1, then tan β/μ ∼ O(1).

Consider a translating Gaussian shaped mass, initially at the shoreline

as in Liu et al. (2003), i.e., the seafloor can be described by h0(x, t) =

exp − (ξ − t)

2



with ξ = 2μ x/ tan β. Once in motion, the mass moves

at constant acceleration. Solution of (A41) under this translating Gaussian

shape is given as

η(ξ, t) = ∞

0

J0(ρξ)ρ



a(ρ) cos(ρt) + 1

ρ

b(ρ) sin(ρt)




dρ +

1

3

(h0 − ξ

∂h0

∂ξ ),

(A43)

in Liu et al. (2003) with

a(ρ) = ∞

0

ςJ0(ρς)



−1

3

(h0 − ξ

∂h0

∂ξ )













t=0, ξ→ς

dς, (A44)

b(ρ) = ∞

0

ςJ0(ρς)



1

3

(2∂h0

∂t + ς

∂2h0

∂ξ∂t)













t=0, ξ→ς

dς. (A45)

Two specific examples are presented in Figs. A9 and A10. Maximum runup

estimates, maximum wave height at x = 0, from the analytical solution

(A43) and the numerical solution of Liu et al. (2003) are given in Fig. A11.

3 Laboratory benchmarking

As a preamble, and before describing the benchmark laboratory tests, waves

are generated in the laboratory by moving vertical paddles. While other generation methods exist, hydraulic paddles allow for the precise and repeatable

specification of arbitrary trajectories. In one-dimensional wave channels—

also known as wave flumes and sometimes wave tanks—there is usually only

one paddle, perpendicular to the water surface. As the paddle displaces horizontally, free-surface waves are generated. In two dimensional wave basins,

there are often multiple paddles that can move horizontally. The experiments to be described were conducted in wave basins with multiple paddles.

This kind of generator is known as the snake generator. When the paddles move independently, its motion resembles the ophite motion of a snake.
NOAA tsunami numerical model standards 35

1.0 1.5 2.0

0 2 4 6 0246 x x

t = 0.1 t = 0.5

t = 1 t = 1.5

1

0

–1

–2

–3

–4

1

0

–1

–2

–3

–4

0

–2

–4

–6

–8

–10

–12

0

–2

–4

–6

–8

–10

–12

0 0.5 0 0.5 1.0 1.5 2.0

y

y

Figure A9: Spatial snapshots of the analytical solution at four different times for a beach slope, β = 5◦,

and landslide aspect ratio, μ = 0.05 (tan β/μ = 1.75). The slide mass is indicated by the light shaded area,

the solid beach slope by the black region, and η by the solid line (Liu et al., 2003).

t = 0.25 t = 1

t = 1.5 t = 1.5

0.3

0.2

0.1

0

–0.1

0.8

0.6

0.4

0.2

0

0.6

0.5

0.4

0.3

0.2

0.1

0

–0.1

0

–5

–10

–15

–20

01234 0246

0 2 4 6 8 10 0246

x x

y

y

Figure A10: Spatial snapshots of the analytical solution at three different times for a beach slope, β = 10◦,

and μ = 0.05 (tan β/μ = 3.5). The slide mass is indicated by the light shaded area, and the solid beach

slope by the black region. In the lower right panel, the water depth and landslide profile is shown at t = 1.5,

the same time as shown in the lower left plot (Liu et al., 2003).
36 Synolakis et al.

–0.2 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6

log(tan β/μ)

Rmax/δ ~ ~

0.55

0.50

0.45

0.40

0.35

0.30

0.25

0.20

Figure A11: Maximum runup as a function of log(tan β/μ). The analytical solutions are shown by the solid line, and the various symbols are from NLSW simulations of Liu et al. (2003), corresponding to different slopes ranging from 2◦ to

20◦.

They are also known as Directional Wave Spectrum Generators (DWSGs).

In the U.S., such DWSGs exist at the Coastal and Engineering Laboratory

of the U.S. Army Corps of Engineers in Vicksburg, Mississippi, at the Texas

A&M University in College Station, Texas, and at Oregon State University

in Corvallis, Oregon.

Here five sets of experiments will be described in detail for laboratory

benchmarking data: solitary wave experiments on a 1:20 sloping beach (Synolakis, 1987), solitary wave runup over a composite beach (Kˆano ̆glu, 1998;

Kˆano ̆glu and Synolakis, 1998), conical island experiments (Liu et al., 1995;

Kˆano ̆glu, 1998; Kˆano ̆glu and Synolakis, 1998), and tsunami runup onto a

complex three-dimensional beach (Takahashi, 1996), and tsunami generation and runup due to three-dimensional landslide (Liu et al., 2005).

3.1 Solitary wave on a simple beach

In this set of experiments, the 31.73 m-long, 60.96 cm-deep, and 39.97 cmwide California Institute of Technology, Pasadena, California wave tank was

used with water at varying depths. The tank is described by Hammack

(1972), Goring (1978), and Synolakis (1986, 1987). The bottom of the

tank consisted of painted stainless steel plates. Carriage rails run along

the whole length of the tank, permitting the arbitrary movement of instrument carriages. A ramp was installed at one end of the tank to model the

bathymetry of the canonical problem of a constant-depth region adjoining

a sloping beach. The ramp had a slope of 1:19.85. The ramp was sealed to

the tank side walls. The toe of the ramp was distant 14.95 m from the rest

position of the piston generator used to generate waves.
NOAA tsunami numerical model standards 37

A total of more than 40 experiments with solitary waves running up the

sloping beach were performed, with wave depths ranging from 6.25 cm to

38.32 cm. Solitary waves are uniquely defined by their maximum height H

to depth d

 ratio and the depth, i.e., H/  d and d are sufficient to specify the

wave. H/  d

ranged from 0.021 to 0.626. Breaking occurs when H/  d > 0.045,

for this particular beach. This is the same set of experiments used to validate

the maximum runup analytical predictions presented in section A2.1.1.

Initial location, Xs in the analytical consideration (A8), changes with

different wave heights. The reason Xs distance varies is that solitary waves

of different heights have different effective “wavelengths.” A measure of

the “wavelength” of a solitary wave is the distance between the point xf

on the front and the tail xt where the local height is 1% of the maximum, i.e., η(xf , t = 0) = η(xt, t = 0) = (H/  d

)/100. The distance Xs is

at an offshore location where only 5% of the solitary wave is already over

the beach, so that scaling can work. Therefore, in the laboratory experiments initial wave heights are identified at a point Xs = X0 + (L/2) where

L/2 = (1/γ) arccosh √20 with γ =



3(H/  d

)/4. In the laboratory, even solitary waves can dissipate. If the wave height is measured far offshore and

used as an initial condition for non-dissipative numerical models, the comparisons will be less meaningful, as the solitary wave will slightly change as it

propagates in the laboratory. By keeping the same relative offshore distance

for defining the initial condition, meaningful comparisons are assured.

While only 10 wave gages were used in each experimental run, the generation was extremely repeatable. As experiments were repeated for each wave

height, the wave gages were moved to different locations, and the same H/  d



wave was generated again until a sufficient number of data points existed to

resolve the entire wave profile. In Synolakis (1987) two different comparisons

are presented: one is the amplitude variation at specific x-locations, and the

second is the amplitude variation at specific t-times, the latter resembling

the image that a photograph from the side with a large depth of field and

angular range would show.

This set of laboratory data has been used extensively for code validation:

refer to Synolakis (1987), Zelt (1991), Titov and Synolakis (1995, 1997,

1998), Titov and Gonz ́alez (1997), Grilli et al. (1997), Li and Raichlen (2000,

2001, 2002). In particular, the data sets for the H/  d

 = 0.0185 (Fig. A4)

nonbreaking and H/  d

 = 0.3 (Fig. A12) breaking solitary waves seem the

most often used and most appropriate for code validation.

3.2 Solitary wave on a composite beach

Revere Beach is located approximately 6 miles northeast of Boston in the

City of Revere, Massachusetts. To address beach erosion and severe flooding

problems, a physical model of the beach was constructed at the Coastal

Engineering Laboratory of the U.S. Army Corps of Engineers, Vicksburg,

Mississippi facility, earlier known as Coastal Engineering Research Center.

The model consists of three piecewise-linear slopes of 1:53, 1:150, and

1:13 from seaward to shoreward. At the shoreline there is a vertical wall

(Fig. A7). In the laboratory, to evaluate the overtopping of the seawall,
38 Synolakis et al.

η

x

t = 10

t = 15

t = 20

t = 25

t = 30

0.40

0.30

0.20

0.10

0

–0.10

0.50

0.40

0.30

0.20

0.10

0

–0.10

0.50

0.40

0.30

0.20

0.10

0

–0.10

0.50

0.40

0.30

0.20

0.10

0

–0.10

0.50

0.40

0.30

0.20

0.10

0

–0.10

0.50

–10 0 20 10

–10 0 20 10

–10 0 20 10

–10 0 20 10

–10 0 20 10

Figure A12: Time evolution of H = 0.3 initial wave over a sloping beach with cot β = 19.85 from t = 10

to 30 with 5 increments. Constant-depth segment starts at X0 = 19.85. Markers show a different realization

of experimental results of Synolakis (1987).
NOAA tsunami numerical model standards 39

the wavemaker was located at 23.22 m and tests were done at two depths,

at 18.8 cm and at 21.8 cm. Water surface elevations were measured with

10 wave gages, and LDV measurements were performed at two different

locations. The maximum runup on the vertical wall was measured visually.

In the experiments, solitary waves of different heights were generated,

and the runup variation for solitary waves striking a vertical wall was determined experimentally. The results were used as benchmark data to validate

1+1 numerical codes in the 1995 Friday Harbor, Seattle, Washington workshop, as discussed in Yeh et al. (1996), and by Kˆano ̆glu (1998) and Kˆano ̆glu

and Synolakis (1998) in validating their analytical formulation.

In terms of specific measurements, time histories of the water surface

elevations exist for H/  d

= 0.038 at x = 12.22 m, H/  d= 0.259 at x = 13.96 m,

and for H/  d

 = 0.681 at x = 14.37 m. These are the locations where the

solitary waves of the given H/  d

are centered for the same reason explained for

the solitary wave experiments over a 1:19.85 beach. In all cases, d

= 21.8 cm.

x = 0 is at one end of the wave flume, i.e., the initial shoreline where the

vertical wall is at x = 23.22 m. Other water-surface elevation time series

were measured midway in each sloping segment, i.e., at x = 17.22 m and

d

= 17.7 cm, at x = 20.86 m and d= 12.5 cm and x = 22.80 m and d= 8.1 cm.

Time histories of water surface elevations are given in Figs. A13, A14, and

A15 for three different wave heights. The measured maximum runup values

are given in Table A1 for the three cases presented here. Maximum runup

measurements for the whole experiments are presented in Fig. A8.

Table A1: Maximum runup measurements on

the laboratory model of Revere Beach.

H/  d

 R R/d

0.038 2.7 cm 0.13

0.259 45.7 cm 2.10

0.681 27.4 cm 1.26

When modeling these experiments, care is needed in calculations in the

near-wall region. As the depth goes to zero, the wave breaks closer to the

shoreline. The breaking wave collapses on the wall, and the air void explodes

producing a splash which cannot be modeled effectively, except with very

high resolution codes.

3.3 Solitary wave on a conical island

Motivated by the catastrophe in Babi Island, Indonesia (Yeh et al., 1994),

during the 1992 Flores Island tsunami, large-scale laboratory experiments

were performed at Coastal Engineering Research Center, Vicksburg, Mississippi, in a 30 m-wide, 25 m-long, and 60 cm-deep wave basin (Fig. A16).

Waves were realistically created in the tank by a horizontal wave generator

with 60 different paddles each 46 cm-wide and moving independently. These

experiments provided runup observations for validating numerical models

and supplemented comparisons with analytical results (Kˆano ̆glu and Synolakis, 1998).
40 Synolakis et al.

0.08

0.04

0.00

0.08

0.04

0.00

0.08

0.04

0.00

0 25 50 75 100 125 150

0 25 50 75 100 125 150

0 25 50 75 100 125 150

η

t

Figure A13: Comparison of the time histories of the free surface elevations midway

in each sloping segment for the analytical solution of Kˆano ̆glu and Synolakis (1998)

and the laboratory data for H/  ̃  ̃d = 0.038 solitary wave.

0.30

0.20

0.10

0.00

0.30

0.20

0.10

0.00

0.30

0.20

0.10

0.00

0 25 50 75 100 125

0 25 50 75 100 125

0 25 50 75 100 125

η

t

Figure A14: The laboratory data for the time histories of the free surface elevations midway in each sloping segment for H/  ̃  ̃

d = 0.259 solitary wave.

0.60

0.40

0.20

0.00

0.60

0.40

0.20

0.00

0.60

0.40

0.20

0.00

0 25 50 75 100 125

0 25 50 75 100 125

0 25 50 75 100 125

η

t

Figure A15: The laboratory data for the time histories of the free surface elevations midway in each sloping segment for H/  ̃  ̃

d = 0.681 solitary wave.
NOAA tsunami numerical model standards 41

Figure A16: View of conical island (top) and basin (bottom).

The detailed experiments are described elsewhere in greater detail (Liu

et al., 1995; Briggs et al., 1995; Kˆano ̆glu, 1998; Kˆano ̆glu and Synolakis,

1998). Briefly, a Directional Spectral Wave Generator (DSWG), located

at x = 12.96 m from the island, generated waves with an initial solitary

wave-like profile. The 27.42 m-long DSWG consists of sixty 46 cm-wide and

76 cm-high individual paddles, each of which can be driven independently.

This allowed performance of experiments with different wave crest lengths.

However, the cases presented here were performed using all the paddles.

Experimental results for different wave crest lengths are given in Briggs et

al. (1995) and Kˆano ̆glu (1998).

In the physical model, a 62.5 cm-high, 7.2 m toe-diameter, and 2.2 m

crest-diameter circular island with a 1:4 slope was located in the basin

(Fig. A17). Experiments were conducted at two different water depths,

32 cm and 42 cm, but presented here with dimensionless solitary wave heights

H/  d

 equal to 0.045, 0.091, and 0.181 at 32 cm. Each experiment was repeated at least twice and maximum runup heights around the perimeter of

the island were measured at 24 locations. Wavemaker signals were presented

in Fig. A19 for these cases to allow direct implementation of these solitary

waves as a wavemaker motion in the numerical models. Water-surface time

histories were measured with 27 wave gages located around the perimeter of

the island (Fig. A18). However, here, time histories of the surface elevation

around the circular island are given at four locations, i.e., in the front of the

island at the toe and gages closest to the shoreline with the numbers 9, 16,

and 22 located at the 0◦, 90◦, and 180◦ radial lines, respectively (Figs. 20

and 21).

These experiments were used as benchmark tests for validating 2+1 nu-
42 Synolakis et al.

 220

720

32 42 62.5

464

384

Figure A17: Definition sketch for conical island. All dimensions are in cm. Not

to scale.

 DIRECTIONAL SPECTRAL WAVE GENERATOR

 9

 7 8

 6

 13 15 17 19 21

10 12 14 16 18 20

 24

27

23

25

y

 

x

TOE OF THE

ISLAND

1800 TRANSECT

00 TRANSECT

900 TRANSECT

2700 TRANSECT

1 2 3 4

22

11

 5

26

....

 

 

 ..

..

. . . .

 

..

.... .... .... 

.

Figure A18: Schematic showing gage locations around the conical island. Not to

scale.

5.0

0.0

–5.0

–10.0

–15.0

8.0

4.0

0.0

–4.0

–8.0

–12.0

–16.0

15.0

10.0

5.0

0.0

–5.0

–10.0

–15.0

0 100 200 300 400

0 100 200 300 400

0 100 200 300 400

H = 0.045

H = 0.091

H = 0.181

Command (cm)

t

Figure A19: Wavemaker motions for the generation of H/  ̃  ̃d = 0.045, 0.091, and

0.181 solitary waves. Target wave heights are given in the insets.
NOAA tsunami numerical model standards 43 0.08 0.04 0.00 –0.04 70 60 50 40 30 20 10 0 Gage 6 0.08 0.04 0.00 –0.04 70 60 50 40 30 20 10 0 Gage 9 0.08 0.04 0.00 –0.04 70 60 50 40 30 20 10 0 Gage 16 0.08 0.04 0.00 –0.04 70 60 50 40 30 20 10 0

t

Gage 22

η

Figure A20: Laboratory data for the time histories ofsurface elevation for a  ̃H/  ̃

d = 0.045 solitary wave at four

gages and  ̃

d = 32 cm. Gage 6 is located at the toe of

the conical island on 0◦ radial line. Gages 9, 16, and 22

are the gages closest to the shoreline on the 0◦, 90◦, and

180◦ radial lines, respectively. Initial wave is defined halfwavelength (L/2, i.e., gages 1 to 4) away from the toe of

the conical island.

0.150.100.050.00

–0.05–0.10

70 60 50 40 30 20 10 0

Gage 6

0.150.100.050.00

–0.05–0.10

70 60 50 40 30 20 10 0

Gage 9

0.150.100.050.00

–0.05–0.10

70 60 50 40 30 20 10 0

Gage 16

0.150.100.050.00

–0.05–0.10

70 60 50 40 30 20 10 0

t

Gage 22

η

Figure A21: Laboratory data for the time histories ofsurface elevation for a  ̃H/  ̃

d = 0.091 solitary wave at four

gages and  ̃

d = 32 cm. Gage 6 is located at the toe of

the conical island on 0◦ radial line. Gages 9, 16, and 22

are the gages closest to the shoreline on the 0◦, 90◦, and

180◦ radial lines, respectively. Initial wave is defined halfwavelength (L/2, i.e., gages 1 to 4) away from the toe of

the conical island.
44 Synolakis et al. 0.30 0.20 0.10 0.00 –0.10 –0.20 70 60 50 40 30 20 10 0 Gage 6 0.30 0.20 0.10 0.00 –0.10 –0.20 70 60 50 40 30 20 10 0 Gage 9 0.30 0.20 0.10 0.00 –0.10 –0.20 70 60 50 40 30 20 10 0 Gage 16 0.30 0.20 0.10 0.00 –0.10 –0.20 70 60 50 40 30 20 10 0 t Gage 22 ηFigure A22: Laboratory data for the time histories of surface elevation for a  ̃H/  ̃d = 0.181 solitary wave at four gages and  ̃d = 32 cm. Gage 6 is located at the toe of the conical island on 0◦ radial line. Gages 9, 16, and 22 are the gages closest to the shoreline on the 0◦, 90◦, and 180◦ radial lines, respectively. Initial wave is defined half- wavelength (L/2, i.e., gages 1 to 4) away from the toe of the conical island. 0.10 0.08 0.06 0.04 0.30 0.25 0.20 0.15 0.10 0.50 0.40 0.30 0.20 0.10 0 60 120 180 240 300 360 0 60 120 180 240 300 360 0 60 120 180 240 300 360 H = 0.045 H = 0.091 H = 0.181 θ RFigure A23: Maximum runup heights from the labora- tory data for three solitary waves  ̃H/  ̃d = 0.045, 0.091, and 0.181 and  ̃d = 32 cm. Initial waves are defined half- wavelength (L/2, i.e., gages 1 to 4) away from the toe of the conical island.
NOAA tsunami numerical model standards 45

Figure A24: Bathymetric profile for experimental setup for Monai Valley experiment.

merical codes in the 1995 Friday Harbor, Seattle, Washington Long-Wave

Runup Models Workshop (Yeh et al., 1996).

3.4 Tsunami runup onto a complex three-dimensional beach; Monai

Valley

The Hokkaido-Nansei-Oki tsunami of 1993 that struck Okushiri Island, Japan,

with extreme runup height of 30 m and currents of the order of 10–18 m/sec

was a disaster, but provided fortuitous high-quality data. The extreme tsunami runup mark was discovered at the tip of a very narrow gulley within

a small cove at Monai. High resolution seafloor bathymetry existed before

the event and when coupled with bathymetric surveys following the event

allowed meaningful identification of the seafloor deformation.

A 1/400 laboratory model of Monai was constructed in a 205 m-long, 6 mdeep, and 3.5 m-wide tank at Central Research Institute for Electric Power

Industry (CRIEPI) in Abiko, Japan and partly shown in Fig. A24. The

laboratory setup closely resembles the actual bathymetry. The incident wave

from offshore, at the water depth d

= 13.5 cm is known. There are reflective

vertical sidewalls at y = 0 and 3.5 m (Fig. A25). The entire computational

area is 5.448 m × 3.402 m, and the recommended time step and grid sizes

for numerical simulations are Δ

x = Δy = 1.4 cm and Δt = 0.05 sec.

The input wave is a LDN with a leading-depression height of −2.5 mm

with a crest of 1.6 cm behind it (Fig. A26). Waves were measured at

thirteen locations and complete time histories are given at three locations,

i.e., (x,  y) = (4.521, 1.196), (4.521, 1.696), and (4.521, 2.196) in meters

(Fig. A27). These experiments were used in the 2004 Catalina Island, Los

Angeles, California NSF Long-Wave Runup Models Workshop.
46 Synolakis et al.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

0

0.5

1

1.5

2

2.5

3

x (m) 

y (m) Input wave 

Okushiri

Island 

Monai

Valley

Figure A25: Computational area for Monai Valley experiment.

0 5 10 15 20 25 −0.02

−0.01

0

0.01

0.02

t (sec)

η (m)

Figure A26: Initial wave profile for Monai Valley experiment.

−2

0

5

Gage 1

−2

0

5

η (cm)

0 20 40 60 80 100 120 140 160 180 200 −2

0

5

t (sec)

Gage 2

Gage 3

Figure A27: Time series of surface elevation at three different gages. Gages 1,

2, and 3 are located at (4.521, 1.196), (4.521, 1.696), and (4.521, 2.196) meters,

respectively.
NOAA tsunami numerical model standards 47

Figure A28: Picture of part of the experimental setup.

3.5 Tsunami generation and runup due to three-dimensional landslide

As discussed, landslide wave generation remains the frontier in terms of

numerical modeling, particularly for aerial slides. These involve not only

the rapid change of the seafloor, but also the impact of the slide on the

shoreline. Details and more experimental results can be found in Liu et al.

(2005).

Large-scale experiments have been conducted in a wave tank with a

104 m-long, 3.7 m-wide, and 4.6 m-deep wave channel with a plane slope

(1:2) located at one end of the tank; part of the experimental setup is shown

in Fig. A28. A solid wedge was used to model the landslide. The triangular

face had a horizontal length of b = 91 cm, a vertical face with a height of

a = 45.5 cm, and a width of w = 61 cm (Fig. A29). The horizontal surface of

the wedge was initially positioned either a small distance above or below the

still water level to reproduce a subaerial or submarine landslide. The block

was released from rest, abruptly moving downslope under gravity, rolling on

specially designed wheels (with low friction bearings) riding on aluminum

strips with shallow grooves inset into the slope. The wedge was instrumented

with an accelerometer to accurately define the acceleration-time history and

a position indicator to independently determine the velocity and position

time histories. Wedge positions are given in Fig. A30 for the two cases

presented here.

A sufficient number of wave gages were used to determine the seaward

propagating waves, the waves propagating to either side of the wedge, and

for the submerged case, the water surface-time history over the wedge. In

addition, the time history of the runup on the slope was accurately measured.

Time histories of the surface elevations and runup measurements for two

cases are presented in Figs. A31 and A32 for the submerged cases with
48 Synolakis et al.

Table A2: Fault parameters of DCRC-17a.

Subfault North Central South

Width (km) 25 25 25

Length (km) 90 30 24.5

Strike (deg) 188 175 163

Dip angle (deg) 35 60 60

Slip angle (deg) 80 105 105

Depth (km) 10 5 5

Dislocation (m) 5.71 2.50 12.00

Origin (◦) 139.40◦E 41.13◦N 139.25◦E 42.34◦N 139.30◦E 42.10◦N

M0 (×1027 dyn cm) 3.85 0.56 2.21

Δ = −0.025 m and Δ = −0.1 m, respectively. While a total of more

than 50 experiments with moving wedges, hemispheres, and rectangles were

conducted, the wedge experiments were used as benchmark tests in the 2004

Catalina Island, Los Angeles, California workshop.

4 Field benchmarking

Verification of a model in a real-world setting is an important part of an

operational model validation. No analytical or laboratory data comparisons

(or any limited number of tests, for that matter) can assure robust model

performance in the operational environment. Test comparisons with realworld data provide an additional important step in the validation of a model

to perform well during operational implementation.

4.1 Okushiri Island

Okushiri Island, Japan, data was the benchmark problem for the 2nd International Long-Wave Runup Models Workshop, Friday Harbor, Washington

(Takahashi, 1996). The magnitude Ms = 7.8 Hokkaido Nansei-Oki, Japan,

earthquake occurred on 12 July 1993 with the hypocenter located at 37 km

depth at 42.76◦N and 139.32◦E off the southwestern coast of Hokkaido,

Japan. The Disaster Control Research Center (DCRC), Japan, digitized

the bathymetric and topographic data from several sources (Fig. A33).

DCRC constructed initial wave profile DCRC-17a with 4.9 m and 2.2 m

uplift and 1.1 m subsidence (depression) as shown in Fig. A34 and given

in Table A2. There are several observations which need to be explained by

numerical modeling:

• Arrival of the first wave to Aonae 5 min after the earthquake should be

estimated with the numerical model. Also, a numerical model should

reveal two waves at Aonae approximately 10 min apart; while the first

wave came from the west, the second wave came from the east. In

addition, two tide gage records at Iwanai and Esashi given in Fig. A35

need to be estimated.

• Maximum runup distribution around Okushiri Island should compare

well with the field measurements (Fig. A36). High runup height at
NOAA tsunami numerical model standards 49

Figure A29: Schematic of the experimental setup.

0 1 2 3 4

0

250

500

Δ = −0.025 m

t (sec)

Wedge location (cm)

0 1 2 3 4

0

250

500

Δ = − 0.1 m

t (sec)

Wedge location (cm)

Figure A30: Time histories of the block motion for the submerged case with

Δ = −0.025 m and −0.1 m.
50 Synolakis et al.

0 2 4 6 8 10

−0.1

−0.05

0

0.05

Gage 1

t (sec)

η (m)

0 2 4 6 8 10

−0.1

−0.05

0

0.05

Gage 2

t (sec)

η (m)

0 2 4 6 8 10

−0.05

0

0.05

Runup gage 2

t (sec)

R (m)

0 2 4 6 8 10

−0.05

0

0.05

Runup gage 3

t (sec)

R (m)

Figure A31: Time histories of the surface elevation and runup measurements for

the submerged case with Δ = −0.025 m.

0 2 4 6 8 10

−0.15

−0.1

−0.05

0

0.05

Gage 1

t (sec)

η (m)

0 2 4 6 8 10

−0.15

−0.1

−0.05

0

0.05

Gage 2

t (sec)

η (m)

0 2 4 6 8 10

−0.08

−0.04

0

0.04

Runup gage 2

t (sec)

R (m)

0 2 4 6 8 10

−0.08

−0.04

0

0.04

Runup gage 3

t (sec)

R (m)

Figure A32: Time histories of the surface elevation and runup measurements for

the submerged case with Δ = −0.1 m.
NOAA tsunami numerical model standards 51

REGION C23

(Δ x = 50 m)

REGION C1 (Δ x = 50 m)

REGION B1 (Δ x = 150 m)

Hokkaido

Okushiri Is.

Honshu

Okushiri Is.

REGION A (Δ x = 450 m)

REGION B

(Δ x = 150 m)

REGION B3 (Δ x = 150 m)

REGION B2 (Δ x = 150 m)

Figure A33: Bathymetric data provided by the Disaster Control Research Center,

Japan.

Hamatsumae, located to the east of Aonae, needs to be explained since

Hamatsumae is sheltered against the direct attack of the tsunami by

the Aonae point. Also, topography does not suggest any tsunami amplification mechanism at this location.

• The highest runup of 31.7 m in a valley north of Monai needs to be

approximated with the numerical model (Fig. A37).
52 Synolakis et al.

Iwanai

Okushiri

Is.

Oshima

 Peninsula

Esashi

upheaval

subsidence

0 50 km

+2.2 m

+4.9 m

–1.1 m

Figure A34: Fault plane DCRC-17a constructed by the Disaster Control Research

Center, Japan.

Iwanai Esashi

t (min) t (min)

η (m) η (m)

1

0

–1

1

0

–1

0 60 120 180 240 300 360 0 10 20 30 40

Figure A35: Tide gage data at Iwanai and Esashi.
NOAA tsunami numerical model standards 53

Figure A36: Maximum runup measurements around Okushiri Island.

Matsue

Hamatsumae

Hoyaishi Point

Monai

505

381

686

837

1234

1360

1921

1870

1584

1618

1164

1027

1071

706

1130

1131

1084

1767

2531

2319

2328

3170

1588

1619

1395

1401

1374

0 500 1000 1500 m

Figure A37: Detailed runup distribution at Hamatsume and north of Monai.
54 Synolakis et al.

4.2 Rat Islands tsunami

Observation data from the expanding tsunami observation network will continue to provide more data for model verification. NOAA’s National Geophysical Data Center (http://www.ngdc.noaa.gov/seg/hazard/tsu.

shtml), NOAA’s Tsunami Warning Centers (http://wcatwc.arh.noaa.

gov/), and NOAA’s Center for Tsunami Research websites provide updated

information on the latest tsunami data. Here, we present the data used for

the first real-time model forecast test as an example of data use for model

verification.

The magnitude Mw = 7.8 earthquake was located on the shelf near Rat

Islands, Alaska, on 17 November 2003 and generated a tsunami. This tsunami provided good data for testing operational models, since the tsunami

was detected by three tsunameters located along the Aleutian Trench and

was recorded at many coastal locations (Titov et al., 2005). This was the

first real-time tsunami detection by the newly developed tsunameter system

(DART). In addition, for the first time, tsunami model predictions were obtained during the tsunami propagation, before the waves had reached many

coastlines. Here, the combined use of tsunami propagation and inundation

models is required for simulation of tsunami dynamics from generation to

inundation. The test requires matching the propagation model data with

the DART recording to constrain the tsunami source model (Fig. A38). If a

finite-difference method on a structured grid is used, several nested numerical

grids would allow “telescoping” from a coarse-resolution propagation model

into a high-resolution inundation model with a model grid of at least 50 m

resolution. If an unstructured grid method is used, a single grid may include

enough resolution near the coast. The data-constrained propagation model

should drive the high-resolution inundation model of Hilo Harbor. The inundation model is to reproduce the tide gage record at Hilo (Fig. A39). Since

this benchmarking is required for the forecasting models, it is essential to

model 4 hr of Hilo Harbor tsunami dynamics in 10 min of computational

time.
NOAA tsunami numerical model standards 55

0.0

1.0

2.0

3.0

4.0

5.0 cm

maximum

offshore

amplit udes

(comput ed)

DART locat ions

DART record

(de-t ided)

D171; max = 2cm

D165; max = 0.5cm

D157; max < 0.2cm

Hilo

Figure A38: Propagation of the 17 November 2003 Rat Islands tsunami. Star indicates epicenter location of the earthquake. Yellow dots are locations of DART systems. White lines near the

DART locations show recorded tsunami signal (detided) at corresponding tsunameter, arrows indicate tsunami arrival on the recordings. Filled colors show example of computed maximum tsunami

amplitudes of a model propagation scenario.

Hilo t ide-gage

Hilo bat hymet ry

November 17, 2003 Rat Is. t sunami at Hilo 0.2

0.1

0.0

-0.1

-0.2

amplit ude, m

3 4 5 6 7 8 9 10

hours aft er eart hquake (06:43:07 UTC, November 17, 2003)

Figure A39: Location of Hilo tide gage and the recording of the 17 November 2003 Rat Islands

tsunami.