NUREG/CR-7223

Tsunami Hazard Assessment: 

Best Modeling Practices and 

State-of-the-Art Technology

Office of Nuclear Regulatory Research

NUREG/CR-7223

Tsunami Hazard Assessment: 

Best Modeling Practices and 

State-of-the-Art Technology

Office of Nuclear Regulatory Research

Manuscript Completed: March 2015 

Date Published: December 2016 

Prepared by: Patrick Lynett1, Yong Wei2, and Diego Arcas2 

 

1Department of Civil and Environmental Engineering

University of Southern California

224D Kaprielian Hall, Los Angeles, CA

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

University of Washington

3737 Brooklyn Ave., N.E., Seattle, WA

Prepared for:

Structural, Geotechnical and Seismic Engineering Branch

Division of Engineering

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555-0001

Rasool Anooshehpoor, NRC Project Manager

NRC Job Code: V6160

iii

ABSTRACT

The present work reviews relevant tsunami dynamics and some of the different mathematical 

models used to describe the physical processes of tsunami propagation in deep and shallow 

water and inundation onto dry land. This report presents a general methodology that 

recognizes some of the unique challenges in conducting tsunami hazard assessment for the site 

of a Nuclear Power Plant (NPP). These challenges relate to the fact that both the peak and 

minimum values of metrics computed in tsunami hazard assessments can negatively affect 

NPPs. For instance, the maximum estimated damage of tsunami impact for most manmade 

structures is generally based on peak values of specific metrics such as wave elevation, flow 

depth, current speed, and flow specific momentum. Particularly for a NPP, minimum values of 

these metrics, such as minimum wave elevation, negative (away from the shore) current speed 

and draw-down are equally dangerous. Conservative tsunami modeling should reflect these 

peculiarities and utilize parameter settings that will maximize or minimize a particular metric in 

order to investigate its effect on the NPP. This report makes recommendations regarding best 

modeling practices for tsunami impact on NPPs, and presents considerations on the most 

appropriate mathematical models for specific situations. Finally, the present work discusses a

strategy to facilitate the identification of a Probable Maximum Tsunami (PMT) using tsunami 

forecast modeling tools currently in operation at the NOAA Tsunami Warning Centers and 

suggests a methodology for probabilistic tsunami hazard assessment.

v

FOREWORD

The 2004 Sumatra tsunami in the Indian Ocean raised the level of concern for an extreme 

tsunami-initiated event, which could potentially exceed the dimensions of all of the recorded 

events taken into consideration in the design basis for those NPPs. Subsequent to the 2004 

Sumatra tsunami, the U. S. Nuclear Regulatory Commission (NRC) conducted an in-depth 

review of past tsunami evaluations and guidelines for nuclear power plants (NPPs) along the 

coastlines. The NRC staff concluded that these facilities are adequately protected.

Consequently, the NRC sponsored a series of research projects at several research institutes to 

further the staff’s understanding of tsunamis and their potential sources so that quantitative 

tsunami wave criteria can be available to assess the tsunami hazards for any prospective 

coastal NPP site. The research institutes were the National Oceanic and Atmospheric 

Administration (NOAA), the U.S. Geologic Survey (USGS), and the Joint Institute for the Study 

of the Atmosphere and Ocean (JISAO) at the University of Washington, Seattle.

This NUREG/CR, prepared by JISAO, provides recommendations for best modeling practices 

for tsunami impact on NPPs, as well as considerations on the most appropriate mathematical 

models for specific situations. This report discusses a strategy to facilitate the identification of a 

Probable Maximum Tsunami (PMT) using tsunami forecast modeling tools currently in operation 

at the NOAA Tsunami Warning Centers and suggests a methodology for probabilistic tsunami 

hazard assessment.

This report reviews relevant tsunami dynamics and some of the different mathematical models 

used to describe the physical processes of tsunami propagation in deep and shallow water and 

coastal inundation. Determination of the PMT requires a comprehensive review and 

investigation of all relevant tsunami sources, followed by a site-specific, state-of-the-art 

hydrodynamic assessment of the tsunami and its effects on a NPP. This report classifies the 

general type of analyses employed as either deterministic or probabilistic. As tsunami hazard is 

difficult to characterize at very low recurrence levels, the deterministic approach is more

commonly used in NPP tsunami hazard assessment.

This NUREG/CR will provide the NRC staff with the means and criteria to assess evaluations 

and analyses for the tsunami design for nuclear facilities provided by the licensees. This 

information will permit the NRC staff to: (1) confirm that adequate levels of safety are 

maintained; (2) improve the effectiveness and efficiency of the review processes; and (3) 

support the staff's technical decisions in a reasonably conservative and realistic manner thereby 

increasing public confidence in the staff’s actions.

vii

TABLE OF CONTENTS

ABSTRACT...............................................................................................................................iii

FOREWORD ............................................................................................................................. v

TABLE OF CONTENTS .......................................................................................................... vii

LIST OF FIGURES ................................................................................................................... ix

LIST OF TABLES..................................................................................................................... xi

EXECUTIVE SUMMARY ........................................................................................................ xiii

ACKNOWLEDGMENTS..........................................................................................................xv

ABBREVIATIONS ..................................................................................................................xvii

1. INTRODUCTION..........................................................................................................1-1

2. BRIEF REVIEW OF TSUNAMI MECHANICS ..............................................................2-1

3. TSUNAMI MODELS: THEORY....................................................................................3-1

Linear Wave Theory..........................................................................................3-1

Nonlinear Shallow Water (Long Wave) Modeling..............................................3-4

High-Order Approaches: Boussinesq-type Model..............................................3-6

High-Order Approaches: Navier-Stokes Modeling.............................................3-8

4. TSUNAMI MODELS: NUMERICAL SOLUTION METHODS........................................4-1

Numerical Solution Methods .............................................................................4-1

Coupled and Hybrid Techniques.......................................................................4-1

Moving Shoreline Algorithms ............................................................................4-2

5. APPLICATION OF NUMERICAL MODELS.................................................................5-1

Specification of Initial Condition ........................................................................5-1

Dynamic Bottom Boundary Conditions..............................................................5-1

Topographic Grid Creation................................................................................5-2

Grid and Time Steps .........................................................................................5-4

Bottom Friction and Turbulence Closure ...........................................................5-4

Antecedent Water Level....................................................................................5-4

Associated Effects ............................................................................................5-5

6. INTERPRETATION AND ANALYSIS OF NUMERICAL RESULTS .............................6-1

Numerical Convergence....................................................................................6-1

Types of Output ................................................................................................6-1

6.2.1 Instantaneous Snapshots of Ocean Elevation and Speed .....................6-1

6.2.2 Time Series ...........................................................................................6-2

6.2.3 Maximum and Minimum Surface Plots ..................................................6-2

6.2.4 Tabular Summary of Results .................................................................6-3

7. EVALUATING TSUNAMI HAZARDS NEAR CRITICAL FACILITIES ..........................7-1

Deterministic Approach.....................................................................................7-1

Probabilistic Approaches...................................................................................7-5

8. PRESENTATION OF TSUNAMI HAZARD ANALYSIS................................................8-1
viii

9. A PRACTICAL APPROACH TO DETERMINING PMT................................................9-1

Background.......................................................................................................9-1

NCTR’s Forecasting Methodology ....................................................................9-2

NCTR’s Propagation Database.........................................................................9-3

Real-time Simulation Computations.................................................................9-5

9.4.1 High-resolution Digital Elevation Models ...............................................9-6

9.4.2 Coverage of the Forecast Models..........................................................9-6

9.4.3 Development of NCTR’s Forecast Models ...............................................9-7

Application to Hazard Assessment....................................................................9-9

Additional Recommendations ......................................................................... 9-11

10. ASSESSMENT OF EARTHQUAKE-GENERATED TSUNAMI INUNDATION 

BASED ON PROBABILISTIC OFFSHORE TSUNAMI HEIGHT...............................10-1

Background and Objectives.............................................................................10-1

Methodology and Procedure ...........................................................................10-2

10.2.1 PTHA Offshore Maximum Tsunami Amplitude.....................................10-2

10.2.2 Tsunami Inundation Model and ComMIT.............................................10-3

10.2.3 Tsunami Source Inversion and Reconstruction....................................10-4

10.2.4 Procedure of Obtaining Probabilistic Tsunami Inundation....................10-5

10.2.5 Relevance for NPPs ............................................................................10-6

Case Study .....................................................................................................10-7

10.3.1 Study Area...........................................................................................10-7

10.3.2 Development of the 2500-year Tsunami Inundation Zone from the 

Probabilistic Offshore Tsunami Height.................................................10-8

11. CONCLUSIONS.........................................................................................................11-1

12. REFERENCES...........................................................................................................12-1

13. GLOSSARY ...............................................................................................................13-1
ix

LIST OF FIGURES

Figure 3.1 Characteristics for water waves in various regimes............................................3-2

Figure 3.2 Schematic of wave transformation through refraction, where the 

change in the blue shading indicates the depth transition ..................................3-3

Figure 5.1 Various types of movements, as a function of material class (taken 

from Cruden and Varnes, 1996.)........................................................................5-3

Figure 6.1 Snapshot of the ocean surface elevation as the tsunami travels across 

the Atlantic .........................................................................................................6-1

Figure 6.2 Snapshot of ocean surface elevation in the nearshore, high-resolution 

grid.....................................................................................................................6-2

Figure 6.3 Maximum elevation (left) and speed (right) predicted by a high 

resolution simulation ..........................................................................................6-3

Figure 7.1 Initial tsunami condition for an earthquake along the Caribbean 

Subduction Zone using the deterministic approach ............................................7-2

Figure 7.2 Initial condition for a tsunami generated by a submarine landslide. 

The "before" landslide profile is given by the solid yellow surface and 

the "after" landslide profile is given by the black dashed line. The 

difference in these two profiles results in the initial sea surface 

elevation profile as shown in solid blue. .............................................................7-3

Figure 7.3 Procedural flowchart for a deterministic evaluation of the PMT ..........................7-7

Figure 7.4 Example procedure tree for the inclusion of slide variability and 

uncertainty into a PTHA ....................................................................................7-8

Figure 8.1 Example of source location map including shoreline and bathymetric 

features..............................................................................................................8-1

Figure 8.2 Sample figure showing bathymetry local to the site of interest ...........................8-2

Figure 8.3 Sample figure showing initial surface conditions and location of the 

transect AAAA′........................................................................................................8-3

Figure 9.1 Typical ground deformation as computed by Okada’s expressions. An 

area of uplift (hot colors) to the North contiguous to an area of 

subsidence (cold colors) to the South can be observed. ....................................9-4

Figure 9.2 Distribution of tsunami maximum amplitude at a random coastal 

location in the Pacific Ocean from each of the potential scenarios 

extracted from NCTR’s propagation data base. A subset of all the 

scenarios represented can be selected and simulated at high-
x

resolution to arrive at a final PMT for the specific site under 

investigation.....................................................................................................9-11

Figure 10.1 2500-year PTHA source disaggregation for a site (118.36°W, 

34.136°N) in California, where the blue bars denote the source 

contributions (%) to the site indicated by the red circle.....................................10-4

Figure 10.2 Flow chart of the methodology to assess tsunami inundation based on 

probabilistic offshore tsunami height ................................................................10-6

Figure 10.3 Aerial photo overlooking Monterey Harbor .......................................................10-8

Figure 10.4 Locations at Monerey Bay, California, where the offshore tsunami 

heights were obtained (from Thio et al., 2010). The circles are areas 

of large water depth discrepancies between Thio et al. (2010) and 

NCTR’s tsunami model. .................................................................................10-10

Figure 10.5 Source disaggregation of 2,500-year offshore tsunami heights at 

Monterey, California courtesy of Thio et al. (2010). The vertical axis 

indicates how much each of the sources contributes in percentage, 

probabilistically, to the tsunami impact at Monterey........................................10-11

Figure 10.6 Source selection in Alaska-Aleutian using ComMIT .......................................10-12

Figure 10.7 Source selection in Kuril-Kamchatka using ComMIT......................................10-12

Figure 10.8 Comparison of the 2,500-year offshore tsunami heights for the M9.5 

Alaska- Aleutian scenario...............................................................................10-13

Figure 10.9 Spectrum analysis of the tsunami waves at all 16 locations where Thio 

et al. (2010) results are available ...................................................................10-15

Figure 10.10 Comparison of the 2500-year offshore tsunami heights for the M9.7 

Kuril- Kamchatka scenario ...........................................................................10-16

Figure 10.11 Computed maximum tsunami wave amplitude in the Pacific due to 

the M9.5 Alaska-Aleutian source....................................................................10-18

Figure 10.12 Computed maximum tsunami water level along Monterey Bay’s 

coastline due to the M9.5 Alaska-Aleutian source ..........................................10-19

Figure 10.13 Computed maximum tsunami flow speed along Monterey Bay’s 

coastline due to the M9.5 Alaska-Aleutian source ..........................................10-20

Figure 10.14 The 2500-year tsunami inundation zone for Monterey Bay, California............10-21
xi

LIST OF TABLES

Table 9.1 List of seismic parameters corresponding to the sea floor deformation 

shown in Figure 9-1............................................................................................9-4

Table 9.2 List of historical events used in the validation of tsunami forecast 

models in the Pacific Ocean by the NOAA Center for Tsunami 

Research ...........................................................................................................9-8

Table 10.1 The 2,500-year tsunami heights at 16 offshore locations obtained by 

Thio et al. (2010). The water depths at these locations are compared 

with those extracted from NCTR’s tsunami model..........................................10-9

Table 10.2 Comparison and errors of the 2500-year offshore tsunami heights 

between ComMIT and Thio et al. (2010) for the M9.5 Alaska-Aleutian 

sources ..........................................................................................................10-14

Table 10.3 Comparison and errors of the 2500-year offshore tsunami heights 

between ComMIT and Thio et al. (2010) for the M9.7 Kuril-Kamvhatka 

source............................................................................................................10-17

xiii

EXECUTIVE SUMMARY

Determination of the Probable Maximum Tsunami (PMT) requires a comprehensive review and 

investigation of all relevant tsunami sources, followed by a site-specific, state-of-the-art 

hydrodynamic assessment of the tsunami and its effects on a Nuclear Power Plant (NPP). PMT

analyses are generally classified as either deterministic or probabilistic. The NPP tsunami 

hazard assessment commonly uses the deterministic approach because tsunami hazard is 

difficult to characterize at very low recurrence levels.

The deterministic approach for site-specific tsunami hazard assessment should employ a 

‘hierarchal’ method, which initially uses extremely conservative assumptions in the tsunami 

hazard analysis. If, after the initial analysis, it is concluded that the PMT might represent the 

Probable Maximum Flood (PMF), these assumptions are removed after careful justification. 

This method of analysis shows how the various conservative assumptions affect the final PMT 

assessment. Although not the most efficient, this approach is systematic and clear; which is

important in justifying a PMT evaluation approach.

Earthquakes and landslides generated from both local and distant locations dominate potential 

tsunami sources. A complete inventory of possible sources from available journal papers, 

conference proceedings and technical reports should be compiled for the NPP site under 

investigation, including geological and geophysical descriptions of the possible sources 

discussed in the Pacific Marine Environmental Laboratory (PMEL) guidelines PMEL-136. With 

the deterministic approach, each source should be described in a physically reasonable yet 

highly conservative manner. While large earthquakes are considered the most common 

tsunami source, the more common sources for the PMT are local landslides, which are 

conceivably capable of generating very large wave heights over a small area. When modeling 

the tsunami generated from a submarine landslide, a range of methods should be utilized, from 

highly conservative to realistic.

For sites where the deterministic PMT is a realistic competitor for the PMF, a Probabilistic 

Tsunami Hazard Assessment (PTHA) will provide a more useful and appropriate result 

compared with other flood hazard analyses. With the completion and publication of the Seaside 

Study (González et al., 2009), PTHA has been established as a viable approach for the 

assessment of tsunami hazards. A major hurdle in the application of PTHA for the evaluation of 

the PMT is the incorporation of landslide-generated tsunamis into the analysis, or more 

specifically, the incorporation of local sources with significant heterogeneity in their spatial 

geometry and time evolution. Driven by a lack of understanding of how small-scale spatial and 

temporal details generate and evolve, if this epistemic uncertainty can be quantified, as it can be 

to a reasonable extent for earthquake sources, then proceeding with a PTHA is justified. For 

landslides, however, this is a great challenge.

xv

ACKNOWLEDGMENTS

The research described in this NUREG/CR was sponsored by the United States Nuclear 

Regulatory Commission (U.S. NRC), Office of Regulatory Research. Rasool Anooshehpoor 

provided guidance during this project and during preparation of this NUREG/CR. The authors 

acknowledge the work of researchers of the NOAA Pacific Marine Environmental Laboratory 

(PMEL) that was used for this report (PMEL Contribution number 4510). This research was 

partially funded by the Joint Institute for the Study of the Atmosphere and Ocean (JISAO) under 

NOAA Cooperative Agreement NA15OAR4320063 (2015-2020), JISAO Contribution No. 2712. 

The authors also acknowledge the NRC technical review team that provided invaluable 

comments and suggestions for improvement of the report presentation. In particular, we would 

like to thank Bhasker Tripathi with NMSS, Yong Li with the NRR, and Michelle Bensi, Henry 

Jones, and Lyle Hibler, with the NRO.

xvii

ABBREVIATIONS

ADCIRC Advanced Circulation

COBRAS Cornell Breaking Wave and Structure

COMCOT Cornell Multi-grid Coupled Tsunami Model

COULWAVE Cornell University Long and Intermediate Wave Model

ComMIT Community Model Interface for Tsunami

DART Deep-ocean Assessment and Reporting of Tsunamis

DEM Digital Elevation Model

FD Finite Difference 

FE Finite Element 

FV Finite Volume 

FUNWAVE Fully Nonlinear Boussinesq Wave Model

GEBCO General Bathymetric Chart of the Oceans

GPS Global Positioning System

JISAO Joint Institute for the Study of the Atmosphere and Ocean

JSCE Japan Society of Civil Engineers, the Nuclear Civil Engineering

LES Large Eddy Simulation

MHW Mean High Water

MOST Method of Splitting Tsunamis

NSW Non-dispersive non-linear Shallow Water

NCTR NOAA Center for Tsunami Research

NOAA National Oceanic and Atmospheric Administration

NPP Nuclear Power Plant

NRC Nuclear Regulatory Commission 

OpenFOAM Open source Field Operation And Manipulation

PMEL Pacific Marine Environmental Laboratory

PMF Probable Maximum Flood

PMT Probable Maximum Tsunami

PTHA Probabilistic Tsunami Hazard Assessment

RANS Reynolds-Averaged Navier-Stokes

RG Regulatory Guide

SIFT Short-term Inundation Forecast Tool

TEPCO Tokyo Electric Power Company

THA Tsunami Hazard Assessment

TUNAMI Tohoku University’s Numerical Analysis Model for Investigation

USGS United States Geological Survey

UTC Coordinated Universal Time (French, Temps Universel Coordonné)

VOF Volume of Fluid

1-1

1. INTRODUCTION

It is the goal of this document to provide background and guidance for the use of numerical

models in site-specific tsunami hazard assessment (THA). This document builds off the effort

presented in the 2007 National Oceanic and Atmospheric Administration (NOAA) Pacific Marine 

Environmental Laboratory (PMEL) Technical Memorandum PMEL-136, “Scientific and Technical 

Issues in Tsunami Hazard Assessment of Nuclear Power Plants.” PMEL-136 provides a 

conceptual framework for the assessment of tsunami hazard near nuclear power plants (NPPs) 

and the specific purpose of this document is to provide a more substantive THA method. While 

future THA methods will most certainly change, the information provided in this report yields 

specific guidance for those undertaking a THA for a NPP.

General Design Criterion 2 in Appendix A of 10 CFR Part 50 states that structures, systems, 

and components (SSCs) important to safety at NPPs must be designed to withstand the effects 

of natural phenomena such as tsunamis without loss of capability to perform their intended 

safety functions. The design bases for these SSCs should reflect appropriate consideration of 

the most severe of the natural phenomena historically reported for the site and surrounding 

area. The design bases should also have sufficient margin to account for the limited accuracy, 

quantity, and period of time in which the historical data were accumulated. These design bases 

may be (a) restraints derived from generally accepted "state of the art" practices for achieving

functional goals, or (b) requirements derived from an analysis (based on calculation or 

experiments, or c) both of the effects of a postulated accident for which an SSC must meet its 

functional goals. The methodology presented in this document is consistent with such 

guidelines. A hierarchal screening approach is detailed here, where different levels of 

conservatism are used depending on the modeled site impacts.

This document focuses on tsunamis generated by landslides because the guidelines discussed 

above and the current and limited understanding of landslide hazards, tsunamis generated by 

potential large, local landslides will be involved and often control the Probable Maximum 

Tsunami (PMT). The methodologies outlined in this document are applicable to any locations 

subject to impulsively generated waves, specifically tsunami created by earthquakes, submarine 

landslides and subaerial landslides. Thus, while this document will often refer to “ocean” waves, 

the associated information is relevant to waves in lakes and rivers, where waves might be 

generated by a local cliff or bluff failure. Furthermore, the term “coastal” used throughout this 

document refers to any boundary between a water body and land.

This document first presents a general background of the important physical process that a 

tsunami undergoes from generation to inundation, followed by the mathematical theories 

employed in the study of ocean waves. Next, the document explains the different types of 

numerical solution methods, from the conventional fixed-grid, finite-difference methods to the 

more recent meshless methods. Approaches used to approximate the motion of the shoreline 

during tsunami inundation are also given. This document also discusses proper model 

configuration, such as creating input bathymetry grids and initial conditions, and the methods of 

analyzing and interpreting the model output. Finally, this document provides reasonable 

methods for THA near a NPP, with a focus on the current deterministic approaches for 

determining a PMT and guidance regarding the level of detail expected for the technical 

presentation of a THA.

2-1

2. BRIEF REVIEW OF TSUNAMI MECHANICS

A tsunami in the deep ocean has a long wavelength and travels quickly (~1000 km/hr). As the 

wave reaches shallow water near the coastline, the tsunami begins to shoal. The speed at 

which a long wave such as a tsunami moves, or celerity, is a function of the local water depth. 

The shallower the depth, the slower the wave moves. A tsunami, with its long wavelength, 

experiences different water depths at any given instant as it travels up a slope; the front of the 

wave, the portion of the tsunami closest to the shoreline, will generally be in the shallower water 

and move slower. The back of the tsunami, on the other hand, will be in deeper water, moving 

faster than the front, causing the wavelength to shorten. As the wave energy is spatially 

compressed into a smaller region, the height of the wave increases. Consequently, despite 

having a height of less than a meter in the open ocean, the tsunami height over land can easily 

exceed several meters. With this increase in wave height comes a more dynamic and complex 

phenomenon as non-linear and dispersive effects become important.

The present discussion will consider a large earthquake-generated tsunami, such as the 2004 

Indian Ocean or the 2011 Tohoku, Japan, earthquake. Other impulsive waves, such as those 

generated by landslide, are more difficult to generalize and will be discussed further at the 

end of this section.

Any wave condition, whether it is a tsunami or a typical wind wave in the ocean, can be 

mathematically described as a superposition, or summation, of a series of sine or cosine waves 

with different amplitudes and wavelengths. For example, with the right choice of individual sine 

waves it is possible to construct even the idealized tsunami: a single solution, or soliton, a wave 

in which the effects of dispersion and nonlinearity compensate for each other, resulting in the 

wave propagating as a single pulse. For a dispersive wave, the various sine wave components 

will have different wave speeds and the wave will disperse as the faster moving components 

separate from the slower ones. For a non-dispersive wave, all the components move at the 

same speed, resulting in no dispersion in the direction of propagation, or spreading of the 

tsunami wave energy. Because the tsunami wave energy will not disperse but will remain in a 

compact pulse, tsunamis can be devastating across such a large spatial region.

The dispersion described above is called "frequency dispersion" as it is primarily dependent on 

the frequency of the wave component. There is another type of dispersion called "amplitude 

dispersion," which is a function of the nonlinearity of the wave, and is usually discussed

under the framework of linear versus nonlinear waves. For tsunamis, the importance of 

nonlinear effects is inferred through the ratio of the tsunami height to the water depth. When 

this ratio is small, such as in the open ocean, the wave is linear; but in shallow waters, the ratio 

is of the order of unity and the wave is no longer linear. The linear/nonlinear nomenclature is

not an intuitive physical description of the waves, but comes from the equations describing the 

tsunami motion, described later in this report. Accounting for this nonlinear effect finds that the 

wave speed is a function of both the local depth and the wave height. Specifically, for two wave 

components of the same period but different amplitudes, the component with the larger 

amplitude will have a slightly greater wave speed. Except for the case of wave fission, the 

nonlinear effect of amplitude dispersion does not disperse tsunami energy with an end result of 

diminishing nearshore impact; instead the nonlinear effect will focus wave energy at the front, 

often leading to a powerful breaking bore.
2-2

Open ocean propagation of a linear, non-dispersive tsunami is a relatively uncomplicated

process that translates wave energy across basins subject to wave speed changes that are a 

function of the local depth. As a tsunami enters the nearshore region, roughly characterized by 

water depths of less than 100 m, the wave can undergo a major physical transformation. The 

properties of this transformation depend heavily on the characteristics of the beach profile and 

the wave itself. Complexities in bathymetry, such as islands, shoals, and canyons, can locally 

amplify or decrease the tsunami energy, sometimes in counter-intuitive ways (e.g., Stefanakis et 

al., 2014). In the simplest inundation case, where the beach profile is relatively steep1 and the 

tsunami wave height is small, the runup process closely resembles that of a wave hitting a 

vertical wall, and the runup height will be approximately twice the offshore tsunami height. In 

these special cases, a breaking bore front is not expected. Instead, horizontal fluid velocities 

near the shoreline would be very small. In this case, the tsunami inundation would closely 

resemble that of a quickly rising tide with only very minor turbulent, dynamic impacts. However, 

even in these cases, overland flow constrictions and other features can create localized 

energetic inundation.

If the beach profile slope is mild, as is typical of continental margins, or the tsunami wave height 

is large, then the shallow water evolution process becomes highly nonlinear. However, while 

the nonlinear effect becomes very important, in the majority of cases, frequency dispersion is 

still very small and can be neglected. Nearshore nonlinear evolution is characterized by a 

strong steepening and possible breaking of the wave front with associated large horizontal

velocities. In these cases, turbulent dissipation can play a major role. Dean and Dalrymple 

(1991) present the technical background and discussion of how beach slope impacts shoaling of 

water waves, as well as the theoretical tools needed to qualify such impacts.

While it may be intuitive to postulate that wave breaking dissipation at the tsunami front plays a 

significant role in tsunami inundation, this may not be entirely correct. This dissipation, while 

intense, is localized at the front and often contains only a small fraction of the total tsunami

energy and momentum. So, for tsunamis such as the 2004 Indian Ocean event, the related 

dissipation likely had only minor impact on leading-importance quantities such as the maximum 

runup and inland (off-beachfront) flow velocities. However, the properties of breaking are 

important to other aspects of tsunami inundation. The maximum forces on beachfront 

infrastructure, such as ports, terminals, piers, boardwalks, and houses, should include bore 

impact as well as the drag force associated with the ensuing quasi-steady flow (Ramsden, 

1996; Yeh, 2006). The bore turbulence may play an important role in understanding how 

bottom sediments are suspended, transported and deposited by a tsunami. Therefore, 

understanding the dynamics of a breaking tsunami front is not of particular importance for near 

real-time or operational tsunami forecast models. Engineers and planners can use this

information to design tsunami-resistant structures.

A second energy dissipation mechanism that plays a major role in determining maximum runup 

is bottom friction. On a fundamental level, this dissipation is caused by the flow interaction with 

the seafloor, where bottom irregularities lead to flow separations and the resulting turbulence. 

All natural bottoms result in some bottom friction; a smooth, sandy beach may generate only 

minor dissipation, while a coral reef or a mangrove forest can play a larger role in reducing 

tsunami energy (Fernando et al., 2005). Other means of energy dissipation will be largely local 

 

1Here "steep" is defined in terms of the tsunami wavelength. If the horizontal distance along the slope connecting 

deep water to the shoreline is small compared to the tsunami wavelength, the beach would be considered steep.
2-3

and may include enhanced mixing due to sediment or debris entrainment, large shallow-flow 

vortex generation by headlands or other natural or artificial obstacles, and flow through/around 

buildings and other infrastructure.

Up to this point, we have only discussed the "typical" nearshore tsunami evolution, which is 

portrayed as a wave without frequency dispersion, and may be called a linear or nonlinear 

tsunami depending on a number of physical properties. The rest of this section will be devoted 

to those situations where the above characterization may no longer be adequate.

Looking first to the tsunami source, waves generated by submarine landslides may not behave 

as non-dispersive waves in the open ocean (Lynett et al., 2003; Weiss et al., 2006). These 

source regions tend to be at least an order of magnitude smaller in spatial extent compared to 

their subduction zone sources. Physically, this implies that the generated waves will be of 

shorter wavelength. As a rule of thumb, if these generated waves have length scales of less 

than 10 times the local depth, then it should be anticipated that frequency dispersion will play a 

non-negligible role (Lynett and Liu, 2002). Under this constraint, individual component wave 

speeds near the dominant period become frequency dependent.

Understanding that an impulsively generated wave can be dispersive has serious implications.

For example, consider a hypothetical landslide located in the Atlantic Ocean that generates a 

dispersive tsunami (e.g., Ward and Day, 2001). As this tsunami travels across the Atlantic, to 

either the east coast of North America or the west coast of Europe, frequency dispersion effects 

will spread the wave energy in the direction of propagation. This will convert the initial shortperiod pulse into a long train of waves. By spreading this energy out, the inundation impact will 

be reduced in two ways. First, by taking a high-density energy pulse and stretching it into a 

longer, lower-density train, the maximum energy flux, and thus the intensity directed at the 

shoreline, will decrease. Second, by increasing the duration of the time series and creating 

many individual crests, energy dissipation can play a bigger role. Using simple energy 

arguments, it can be shown that a low-density, long-period train loses more energy through 

bottom friction and wave-breaking than a high-density, short-period pulse. Numerical studies 

have shown that for such cases, the individual wave crests are largely dissipated, and runup is 

dominated by the carrier wave, thus becoming a time-dependent, wave-setup problem 

(Korycansky and Lynett, 2007).

Some current research in the tsunami community argues that frequency dispersion may 

occasionally play a non-negligible role in even the long-wavelength, subduction-zone tsunamis. 

Two arguments are used to justify this assumption. First, that short-period energy generated at 

the source is significant and will lead to different patterns of runup if included (e.g., Kulikov, 

2006; Horrillo et al, 2006). And second, that shallow-water nonlinear interactions can generate 

short-period components which can become decoupled (or unlocked) from the primary wave, 

and change the incident tsunami properties (Matsuyama et al., 2007).

Under certain conditions, namely a nonlinear tsunami propagating across a wide shallow shelf, 

a process called fission may occur. Wave fission is a separation process where wave energy, 

initially part of a primary wave or pulse, attains certain properties, such as higher or lower phase 

speed, allowing it to disconnect from the primary wave and propagate as an independent wave. 

In the context of nearshore tsunami evolution, a standard mechanism is the cause of this 

fission. In order to explain wave fission, it is necessary to describe what a nonlinear, phaselocked wave is. To do this, we will examine the acceleration terms of the 1D inviscid flow 

conservation of momentum equation:
2-4

uutt + uuuuxx = −ppxx,

where uu(xx,tt) and pp(xx,tt) are the flow velocity and pressure at location xx and time tt, respectively; the 

subscripts indicate partial derivatives with respect to time tt and coordinate xx. Now assume that

there is a single wave component, under which the velocity oscillates as:

uu(xx,tt) = cos(kk − σσ )

where kk is the wavenumber, σσ the angular frequency, and σσ/kk the phase velocity of the wave.

If the wave is nonlinear, which is to say that the convective acceleration term (uuuuxx) in the above 

momentum equation is not negligible, the convective term will become:

uuuuxx = −kkcos(kk − σσtt) × sin(kk − σσtt)

= − kk

2 sin(2kk − 2σσtt)

Thus, this nonlinear term generates a new wave component with twice the wavenumber and 

frequency (or half the wavelength, and period). From linear wave theory, it is expected that this 

new wave with a shorter period will have a different wave speed than the original, primary wave. 

However, from the phase function of this new wave, there is the speed 2σσ⁄2kk = σσ⁄kk, which is 

identical to the speed of the primary wave. This new wave is locked to the phase of the original 

wave; this is the fundamental effect of nonlinearity. This locking connection can be delicate, 

and any disruptions to the primary wave, such as a varying seafloor, dissipation, or interactions 

with other free waves in the train or wave pulse, can cause the new waves to become unlocked. 

When this occurs, the now free waves retain their frequency,~2σσ, but take a wavenumber as

given by the linear theory dispersion relation. Since these freed waves will be of a shorter 

period than the primary wave, they will travel at a slower speed and generally trail the main 

wave front.

Long wave fission is most commonly discussed in the literature by describing a solitary wave 

propagating over an abrupt change in depth, such as a step (e.g., Madsen and Mei, 1969; 

Johnson, 1972; Seabra-Santos et al., 1987; Losada et al., 1989; Goring and Raichlen, 1992; Liu 

and Cheng, 2001). In these cases, there is a deep-water segment of the seafloor profile, where 

a solitary wave initially exists. At this depth, the solitary wave is of permanent form. As the 

solitary wave passes into shallower water, the leading wave energy will try to rediscover a 

balance between nonlinearity and dispersion. Since this new solitary wave will be of a different 

shape and contain a lower mass, by conservation law there must be some trailing disturbance 

to account for the mass deficiency. This trailing disturbance will take the form of a rank-ordered 

train of solitons. The solitons in the trailing train, while smaller in height than the leading solitary 

wave, tend to have a similar wavelength; this has been shown both analytically, numerically, 

and experimentally. However, the discussion of fission in this sense is not particularly relevant 

to ’real’ tsunami modeling, where the offshore wave approaching the shelf break rarely 

resembles a solitary wave solution (Tadepalli and Synolakis, 1996). However, the offshore 

wave does not specifically need to be a solitary wave for this process to occur; any long,

arbitrary pulse of energy with some nonlinearity will undergo this process.

Numerous eyewitness accounts and videos of the 2004 Indian Ocean and 2011 Japan tsunami

provide evidence of a tsunami approaching the coastline as a series of short period (on the 

order of 1 minute and less) breaking fronts, or bores (e.g., Ioualalen et al., 2007). These short 

period waves may be the result of fission processes of a steep tsunami front propagating across 

a wide shelf of shallow depth. Along the steep front of a very long period wave, nonlinearity is
2-5

important. There is a large amount of energy in high-frequency components with wavelengths 

similar to the horizontal length of the tsunami front (on the order of 1 km). As the wave 

continues to shoal, the high-frequency locked waves may eventually become free waves, and 

take the form of very short waves "riding" the main wave pulse. This situation is akin to an 

undular bore in a moving reference frame. This process, identical to that described in the above 

paragraph, simply takes place over a much longer distance. The newly freed waves, in the 

nonlinear and shallow environment, will attempt to reach an equilibrium state, where frequency 

dispersion and nonlinearity are balanced. Thus, the fission waves will appear as solitary waves, 

or more generally, cnoidal waves. This fact provides some guidance as to the wavelength of 

these fission waves; they can be approximately calculated via solitary wave theory using the 

tsunami height and depth of the shelf. A conclusion of this fission issue is that, if attempting to 

simulate tsunami propagation with dispersive equations, and if the grid is too coarse to resolve 

the short fission waves, the justification to use the dispersive model is greatly degraded.

It is well established that large-scale coastal features, such as small islands, large shoals, 

canyons and shelves play an important role in tsunami inundation due to conventional shallow 

water effects such as shoaling and refraction (e.g., Carrier, 1966; Yeh et al., 1994; Briggs et al,

1995; Liu et al., 1995; González et al., 1995). On the other hand, understanding the impact of 

small-scale features is in early development. Efforts to understand the effect of small-scale 

features was largely initiated by field observations. Synolakis et al (1995), while surveying the 

coast of Nicaragua for information about the 1992 tsunami in the region, noted that the highest 

level of damage along a particular stretch of beach was located directly landward of a reef 

opening used for boat traffic. It was postulated that the reef gap acted as a lower resistance 

conduit for tsunami energy, behaving like a funnel and focusing the tsunami. In contrast, along 

neighboring beaches with intact reefs, the tsunami did not have the intensity to remove beach 

umbrellas. Investigating impacts from the same tsunami, Borrero et al. (1997), discussed how 

small-scale bathymetry variations affected coastal inundation. One of the conclusions of this 

work was that bathymetry features with length scales of less than 50 m had leading-order 

impact on the runup.

Looking to the Indian Ocean tsunami, a survey team in Sri Lanka inferred from observations that 

reef and dune breaks lead to locally increased tsunami impact (Liu et al., 2005). Also in Sri 

Lanka, Fernando et al. (2005) performed a more thorough survey along the southeastern 

coastline, and concluded that there was a correlation between coral mining and locally severe 

tsunami damage. While additional research is needed to quantify the effects of small-scale 

features, the observations suggest that defense measures such as seawalls, once thought to be 

inconsequential to tsunami inundation, may actually provide some protection as it has been 

observed that small barriers can have large effects on local tsunami elevation and speed.

Onshore, tsunami propagation is affected by the general topography (ground slope), ground

roughness, and obstacles (e.g., Synolakis, 1987, Tadepalli and Synolakis, 1994, Lynett, 2007; 

Tomita and Honda, 2007). The composition of the ground, which may include sand, grass, 

mangroves or pavement, controls the roughness and the subsequent bottom friction damping. 

To predict tsunami inundation with high confidence, the ground type must be well mapped and 

the hydrodynamic interaction with that type must be well understood. If the tsunami approaches 

the shoreline as a bore, the process of "bore collapse", or the conversion of potential to kinetic 

energy, will cause the fluid to rapidly accelerate (Shen and Meyer, 1963; Yeh et al., 1989). This 

fast flow equates to high fluid forces on obstacles such as buildings. Tsunami interaction with 

these obstacles can lead to a highly variable local flow pattern (e.g., Cross, 1967; Tomita and 

Honda, 2007). As the flow accelerates around the corners of a building, for example, the scour 

potential of that flow increases greatly, and foundation undermining is a concern. As with any 
2-6

fluid flow past an obstacle, the backface of the obstacle is characterized by a low-pressure 

wake. Combined with the interior flooding of a building, this low-pressure wake may lead to an 

outward "pull" force on the back wall, causing it to fail by falling away from the center of the

building. Increasing the topographical complexity, structures built in coastal environments are 

located within close enough proximity to each other such that their disturbances to the flow may 

interact. This can lead to irregular and unexpected loadings (e.g., a second row building

experiences a larger force than beachfront buildings due to a funneling effect).
3-1 

3. TSUNAMI MODELS: THEORY

Linear Wave Theory

The simplest methods for estimating coastal wave properties utilize linear potential wave theory. 

While this theory is covered comprehensively in many textbooks (e.g., Dean and Dalrymple,

1984), the most relevant equations will be given here. Additionally, these equations will

establish the basic variable quantities most often used in coastal hydrodynamic wave modeling. 

The free surface elevation function, ηη, for a single frequency component is given by: 

ηη = HH

2 cos(kk − σσ ),

where HH is the wave height equal to twice the wave amplitude, aa; kk is the wave number; xx is 

the horizontal coordinate; σσ is the angular frequency; and tt is time. In general terms, kk and xx

can be vector quantities when directional waves are considered. Wavenumber and angular 

frequency are related to other important physical quantities by: 

kk = 2ππ

λλ , σσ = 2ππ

TT , ff = 1TT

with λλ is the wavelength, or horizontal distance between two successive crests or troughs, TT is 

the wave period, and ff is the wave frequency. The linear dispersion relation provides the 

relation between wave frequency and wavelength by: 

σσ2 = gg tanh(kkh)

where gg is the acceleration due to gravity and h is the local water depth. Expressions for fluid 

velocity and pressure under the wave can be found in textbooks, and will not be repeated here. 

Linear wave theory predicts closed, elliptical orbits of fluid particles under waves, and therefore 

no net mass transport.

In the limiting case of very shallow water which is often relevant to tsunami studies, here 

implying the ratio of water depth to wavelength goes to zero, simplifications can be made to the 

full linear potential theory discussed above. The dispersion relation reduces to

σσ2 = gghkk2

which yields a wave speed, cc, of

cc = σσ

kk = λλ

TT = �ggh

and a horizontal fluid velocity, uu, of

uu = HH

2 �gg

h cos(kk − σσ ).

It is remarked here that the horizontal fluid velocity, uu, also known as the particle velocity or 

orbital velocity, is typically much smaller than the speed at which the wave moves, cc. In shallow 

water, these two quantities can be related by 
3-2 

max(uu) = aa

h �ggh

In the open ocean, the aa/h parameter will often be less than 1/1000, indicating there is a large 

difference between the speed of the wave and the speed of the water under the wave. It is 

usually not until the tsunami reaches the coastal area that the aa/h parameter might approach

unity.

The expressions in the section should be applied only in the case of a very long wave, such as 

for tides and tsunamis. However, even outside of these cases, quick wave property calculations 

based on the above will still provide a useful preliminary reference point for coastal waves. As 

mentioned in the previous section, tsunamis are typically non-dispersive, but can contain 

dispersive frequency components depending on the properties of the source. Similarly, the 

wave may be linear or nonlinear, depending largely on the amplitude, or height, of the waves. 

To describe these characteristics for water waves, Figure 3-1 shows the various regimes of 

water waves. From this figure, we see that shallow water (non-dispersive) theory is valid to the 

left of the red-dashed line, for h/(ggTT2) less than 0.0015. The region where waves can be 

considered as linear, shown by the yellow shaded area, intersects the shallow water region only 

for very small HH/(ggTT2). It is in this seemingly small portion of the overall parameter space in 

which water waves exist that earthquake-generated tsunamis in the open ocean fall. As 

tsunamis enter the coastal area, the waves may become nonlinear. Many tsunami models, 

such as the MOST and ComMIT models discussed later in this report, are accurate throughout 

the shallow water region, for both linear and nonlinear waves.

Figure 3.1 Characteristics for water waves in various regimes
3-3 

Linear wave transformation includes the effects of shoaling, refraction and diffraction. Wave 

refraction is the process of wave crests bending due to propagation-normal variations in 

bathymetry (Figure 3.2). Wave crests will bend "into" shallow water regions, with the refraction 

coefficient KKrr (Dean and Dalrymple, 1991): 

KKrr = �cosαα0

cosαα

where αα0 is the angle of incidence relative to the shoreline, or some prescribed depth contour, 

and is related to the refraction angle αα through the Snell’s Law:

sinαα = cc

cc0

sinαα0,

where cc and cc0 are the wave speeds in the shallow and deep waters, respectively. Once a 

refraction coefficient is obtained, the transformation experienced by the wavelength can be 

computed from the relation below: 

KKrr = �nn0λλ0

nn

where the subscript “0” denotes deep water values and

nn = 1

2 �1 +

2kkh

sinh(2kkh)

�.

The "nn" value varies from 0.5 for a deep water wave to 1.0 for a shallow water wave. The basic 

linear theory presented here does not include bottom reflection, and thus can be applied 

between any two points regardless of the bathymetry between these two points. However, the 

Figure 3.2 Schematic of wave transformation through refraction, where the change in 

the blue shading indicates the depth transition
3-4 

theory implicitly assumes an inviscid propagation and a mildly sloping bottom, and this must be 

satisfied at a minimum for the shoaling estimation above to be physically appropriate.

These techniques have been displaced by numerical methods, some of which are discussed in 

later sections of this chapter. There is no general equation for the diffraction coefficient. The 

engineer should reference the numerous Diffraction Diagrams that can be found in the Coastal 

Engineering Manual (CEM) (2002), which plot the diffraction coefficient for various combinations 

of breakwater configurations and incident wave conditions.

Once the wave reaches the break point, the coastal transformation equation given above can no 

longer be used, and some other approach must be used to cap the wave height. A common 

engineering approach is the introduction of the "breaker index", which is the ratio of the breaking 

wave height to the local water depth. For a given wave period and bathymetry profile, a breaker 

index can be estimated. Once the shoaling and refraction analysis estimates a wave-height-todepth ratio that exceeds the breaker index, it is assumed that the wave starts to break; this is 

the break point. Breaking is then usually assumed to be "depth-limited", meaning that the 

breaker index remains constant as the water depth continues to decrease onshore. The 

difficulty in this approach is the reliable estimation of the breaker index. There are various

equations for this index; the CEM contains the most common equations. For waves shoaling 

over a very steep slope, the breaker index can be greater than 1.0, while for waves breaking 

over an essentially flat bottom, the breaker index may be lower than 0.4 (e.g., Raubenheimer et 

al., 1996). 

For studies that use straightforward analysis methods but also estimate the wave height 

envelope across a cross-shore beach profile, the linear Energy Flux method (e.g., Dally et al., 

1985) has shown to be a useful tool. Here, it is assumed that the spatial rate of change of wave 

energy flux, FF, is due to some dissipation function, φφ: 

dd

dd = −φφ(xx).

Breaking drives the dissipation function and is proportional to the difference between the local 

wave energy flux and some "stable" energy flux for the local depth. Determination of a stable 

energy flux, similar to the determination of a breaker index, is difficult for a wide range of 

conditions. Note that outside the breaker zone, where the left hand side of the above equation 

is zero, the Energy Flux method will reduce to a linear shoaling model. A number of large-scale 

morphology change simulators have adopted this wave transformation model, which require a 

good estimate of the breaker location and resulting longshore currents for large horizontal 

areas. However, this approach has been largely displaced by more recent numerical tools (e.g. 

Wei et al., 1995), which can better include directionality, dissipation, and nonlinear effects. 

Some of these tools are discussed in the following sections.

Nonlinear Shallow Water (Long Wave) Modeling

A long wave has a wavelength much larger than the local water depth. Waves that commonly 

meet this requirement include tides, storm surges, and tsunamis. When the wave disturbance is 

"long", it is reasonable to use the non-dispersive non- linear shallow water (NSW) wave 

equations. The NSW equations can be derived in a number of different ways, but all 

fundamentally arise from an integration of the Euler or Navier-Stokes equations with an 

assumption of vertically-invariant horizontal velocity and hydrostatic pressure. Due to the 
3-5 

simple and well-studied nature of the NSW equations, they can be solved by a wide variety of 

numerical schemes.

Models used to examine tides and storm surges are often similar, if not the same. Examples of 

this class of model include ADCIRC (Dawson et al., 2006) and DELFT3D (Deltares Systems), 

both of which have received widespread acceptance in the engineering community. These 

models need to include proper tidal forcing, wind stress for surge, and bottom friction in 

nearshore areas. Of particular importance for accurate surge prediction is high-resolution and 

accurate coastal bathymetry and topography, often down to 10 m or less. Therefore, a very 

wide range of scales must be simulated, a common challenge for simulation of long-wave 

generation and evolution. Numerical approaches that can accommodate this issue, such as 

finite element meshing or grid nesting, must be adopted.

With the large number of deadly tsunamis in the past decade, tsunami simulation capabilities

have increased relatively rapidly. The National Tsunami Hazard Mitigation Program, sponsored 

by NOAA, currently uses several tsunami computational models to produce tsunami inundation 

and evacuation maps for the states of Alaska, California, Hawaii, Oregon, and Washington. 

Tsunami computational models include: MOST (Method Of Splitting Tsunami), originally 

developed by researchers at the University of Southern California (Titov and Synolakis, 1998); 

COMCOT (Cornell Multi-grid Coupled Tsunami Model), developed at Cornell University (Liu et 

al., 1994); and TUNAMI (Tohoku University’s Numerical Analysis Model for Investigation), 

developed at Tohoku University in Japan (Imamura et al., 1988). All three models solve the 

same depth integrated and 2D horizontal (2DH) nonlinear shallow-water equations with different 

finite difference algorithms.

Successful simulation of tsunami propagation and accurate prediction of the arrival time and 

wave height at different locations rely on a correct estimate of the earthquake fault plane

geometry and fault type. Interplate faults in subduction zones are responsible for most of the 

large historical tsunamis. For interplate fault ruptures, the resulting seafloor displacement can 

be estimated using linear elastic dislocation theory (e.g., Manisha and Smylie, 1971; Okada, 

1985). More sophisticated fault models expect non-uniform stress fields (i.e., faults with various 

kinds of barriers, asperities, etc.; e.g., Kanamori, 1977) leading to a seafloor displacement field 

with large spatial variability. Regardless of the earthquake source and fault model, once the 

seafloor displacement is determined, the initial ocean free surface profile is assumed to take the 

same configuration, based on the assumptions that the upward seafloor movement is impulsive 

and seawater is incompressible.

For a given source-region condition specified by either the initial free surface elevations or a 

time history of seafloor displacement, hydrodynamic models can accurately simulate 

propagation of a tsunami over a long distance, provided that bathymetry data exists. The 

shallow-water equation models by definition lack the capability of simulating dispersive waves, 

which could well be the dominating features in landslide-generated tsunamis (Lynett and Liu, 

2002) and for tsunamis traveling a long distance (Grilli et al., 2007). To address this issue of 

dispersivity, a different set of governing equations must be employed, or some manipulation of 

the numerical truncation error must be made.

As a tsunami propagates into the nearshore region, the wave front undergoes a nonlinear 

transformation while it steepens through shoaling. If the tsunami is large enough, it can break at 

some offshore depth and approach land as a bore. Wave breaking in traditional NSW tsunami 

models has not been handled in a physically satisfactory manner. Numerical dissipation is 

commonly used to mimic breaking, thus producing grid-dependent results. An alternative is to 
3-6 

include ad-hoc breaking schemes to the governing equations, which been validated for a wide 

range of nearshore conditions (e.g., Chen et al., 2000).

High-Order Approaches: Boussinesq-type Model

A significant effort in the nearshore wave model community towards developing phaseresolving, time domain Boussinesq models has occurred in the past decade. A Boussinesq 

model is similar in mathematical form to the NSW equations, yet are able to include frequency 

dispersion. Assuming that both nonlinearity and frequency dispersion are weak and are in the 

same order of magnitude, Peregrine (1967) derived the "standard" Boussinesq equations for 

variable depth in terms of depth-averaged velocity and free surface displacement. Numerical 

results based on standard Boussinesq equations or equivalent formulations have given 

predictions that compared quite well with field data (Elgar and Guza, 1985) and laboratory data 

(Goring, 1978; Liu et al., 1985).

As mentioned above, the standard Boussinesq equations are derived based on an assumption 

between the nonlinearity of the wave εε = aa/h, and the frequency dispersion of the wave, μμ = kkh. 

The precise relation of these parameters comes from a non-dimensionalization of the full 

potential flow equations using a shallow water scaling, and is given by:

OO(εε) = OO(μμ2) ≪ 1,

which is the "true" Boussinesq assumption as it relates to nonlinear long surface waves. With

this assumption, nonlinear 2D-vertical potential flow can be reduced to the 1D-horizontal

equation set: 

continuity:; 

ηηtt + εε(ηηφφxx )xx + hφφxxxx = 0

momentum: φφtt + ηη +

εε

2φφ2

xx − μμ2

3 φφxxxx = OO(εεμμ2, μμ4)

where the subscripts represent partial derivatives and φφ is the depth-averaged velocity 

potential. Note that the horizontal velocity, uu, is calculated as the spatial gradient of potential, 

uu = φφxx. The non-dispersive, shallow-water wave equations are easily extracted from this 

equation set by neglecting the term on the left hand side of the momentum equation (φφtt); it is 

this single term that adds dispersion to the shallow water model. 

The mathematical effect of the additional dispersive term can be quantified by examining the 

linear form of the standard Boussinesq equations given above, and substituting in the linear 

wave solution form:

ηη = aaeeiiii φφ = iiΦeeiiii θθ = kk − σσ

Where Φ represents the amplitude of the velocity potential function. The dispersion relation of 

this approximate equation set can be expressed as:

σσ2 = gghkk2

1 + 1

3 (kkh)2
3-7 

which is also the [0, 2] Pade approximation of the full linear dispersion relation. In a practical 

sense, this additional φφxxxx term in the momentum equation allows for accurate linear 

propagation of waves up to kkh ∼ 0.5, which is approximately a ten-fold increase in applicability 

when compared to the shallow water model. This increase in dispersion accuracy is important 

for tsunamis that might contain energy in dispersive frequencies, such as landslide-generated 

tsunamis.

Because it is required that both frequency dispersion and nonlinear effects are weak and of the 

same order, the standard Boussinesq equations are not applicable to very shallow water depth, 

where the nonlinearity becomes more important than the frequency dispersion, and to the deep 

water depth, where the frequency dispersion is of order one. The standard Boussinesq 

equations break down when the depth is greater than one-fifth of the equivalent deep-water 

wavelength. A lesser depth restriction is desirable for many engineering applications where the 

incident wave energy spectrum consists of many frequency components. To extend the 

applications to shorter waves (or deeper water depth) many modified forms of Boussinesq-type 

equations have been introduced (e.g., Madsen et al., 1991; Nwogu, 1993; Chen and Liu, 1995). 

Although the methods of derivation are different, the resulting dispersion relations of the linear 

components of these modified Boussinesq equations are similar, and may be viewed as a slight 

modification of the [2, 2] Pade approximation of the full dispersion relation for linear water waves 

(Witting, 1984). It has been demonstrated that these modified Boussinesq equations are able to 

simulate wave propagation from intermediate water depth (water depth to wavelength ratio of

about 0.5) to shallow water depth including the wave-current interaction (Chen et al., 1998).

Despite the success of the modified Boussinesq equations in intermediate water depth, these 

equations are still restricted to weak nonlinearity. As waves approach shore, wave height 

increases due to shoaling until eventually breaking. The wave height to water depth ratios 

associated with this physical process violates the weakly nonlinear assumption. By eliminating 

the weak nonlinearity assumption this restriction is removed (e.g., Liu, 1994; Wei et al., 1995). 

Numerical implementations of the highly-nonlinear, Boussinesq-type equations include 

FUNWAVE (Fully Nonlinear Boussinesq Wave Model, e.g., Wei et al., 1995) and COULWAVE 

(Cornell University Long and Intermediate Wave Model, e.g., Lynett & Liu, 2002). These 

models have been applied to a wide variety of topics, including rip and longshore currents 

(Chen et al., 1999; Chen et al., 2003), wave runup (Chen et al., 2000), wave-current interaction 

(Ryu et al., 2003), and wave generation by submarine landslides (Lynett & Liu, 2002), among 

many others. Boussinesq models are steadily becoming a practical engineering tool. 

Directional, random spectrums can readily be generated by the models, which capture 

nearshore evolution processes such as shoaling, diffraction, refraction, and wave-wave 

interactions with very high accuracy. Limiting the applicability of the Boussinesq models is the 

fact that the models are fundamentally inviscid. The dissipation processes, such as breaking 

and bottom friction, must be parameterized in traditional Boussinesq models. While these 

processes are parameterized, the models are able to provide an accurate prediction of 

associated effects for a wide range of wave types, including tsunamis (e.g. Lynett, 2006).

Recently, a number of "non-traditional" Boussinesq approaches have been developed, with the 

goal of including horizontal vorticity explicitly in the flow field. Veeramony and Svendsen (2000)

attempt to include these dynamics under a breaking wave with further advances given in 

Musemeci et al. (2005). Integrated within a Boussinesq-type derivation, the stream function 

equation is used to determine the vertical variation of the velocity. This allows the inclusion of 

the vorticity generated by breaking. The breaking terms that appear as corrections to the 

momentum balance are a function of the amount of vorticity generated during the breaking 

process. This vorticity is obtained from the solution to the vorticity transport equation, and has 
3-8 

been shown to capture the flow field dynamics under a spilling breaker. Kim et al. (2009) made 

a similar attempt to include the viscous effects of a bottom shear and the associated rotationality 

directly in a Boussinesq-type derivation. While this leads to a far more complex equation model, 

it includes the physics necessary to simulate boundary shear and the complete coupling of 

these effects with a nonlinear, dispersive wave field. This model can predict the friction-induced 

changes to the vertical profile of velocity under weakly unsteady flow and thereby provide good 

estimates of internal kinematics. This model is also able to translate the bottom-created 

horizontal vorticity into a vertical vorticity field. Furthermore, such a model is able to couple the 

dissipative and nonlinear effects of a bottom shear with a dispersive and nonlinear wave field. 

While this and similar recent advances allow the depth-integrated equations to model a wider 

range of physical processes, the computing cost for this inclusion is high.

High-Order Approaches: Navier-Stokes Modeling

Phase- and depth-resolving surf zone hydrodynamic models, such as those that use the 

Reynolds-Averaged Navier-Stokes (RANS) equations along with a turbulence closure and a 

robust free surface-tracking scheme, are an ideal alternative for simulation of complicated 

nearshore processes that involve breaking waves. In general, the RANS-based models are 

capable of calculating turbulence energy and energy dissipation due to wave breaking and 

bottom friction. For example, in one of the RANS-based models, COBRAS (Cornell Breaking 

Wave and Structure; e.g., Lin and Liu, 1998 a, b; Hsu et al., 2002), the two-dimensional RANS 

equations are coupled with the kk − εε turbulence closure and volume of fluid (VOF) method for 

tracking free-surface position. This and similar models can adequately resolve the wave 

breaking processes and its interactions with the seabed.

This general class of model is still not practically applicable, even in the academic sense, for 

simulating waves in 3D coastal regions. 3D solvers, using various turbulent closure schemes, 

are however becoming a common basis for large-basin numerical wave tanks. One of the more 

promising programs is OpenFOAM, which is an open-source and freely available computational 

fluid dynamics package that includes numerical solvers for a large number of equation sets, 

including RANS and Large Eddy Simulation (LES). In the coming decade, it is anticipated that 

wind wave modeling will continually move towards the utilization of very high-resolution LES 

tools, as this class of model is the main option for predicting vertical and turbulent structure 

which drives mixing and transport.
4-1 

4. TSUNAMI MODELS: NUMERICAL SOLUTION METHODS

Numerical Solution Methods 

Solution methods are generally classified into grid/mesh methods and meshless/particle 

methods. Grid methods will include finite difference, finite volume, and finite element methods. 

Meshless methods cover Lagrangian approaches (e.g., Zelt, 1991) and smooth particle 

hydrodynamics (e.g., Dalrymple and Rogers, 2006). These meshless methods are not 

commonly used in tsunami hazard assessment, and are not discussed in detail here. 

At present, grid methods tend to be the more common approach found in operational models, 

with the finite difference option being the most widespread. Finite difference models solve the 

governing partial differential equations at a discrete point in space and/or time, where the local 

derivatives are approximated through Taylor-series-derived difference equations. The type of 

difference equation used controls the accuracy of the model; very high order difference 

equations can be employed to produce a solution with equally small numerical error. 

Established models such as MOST (Titov and Synolakis, 1998), COMCOT (Liu et al., 1994), 

TUNAMI (Imamura et al. 1988), and NEOWAVE (Yamazaki et al., 2010) are examples of fixedgrid, finite difference schemes. Upwind differencing of advection terms in the momentum 

equation is often employed to increase model stability (e.g., COMCOT and TUNAMI). However, 

this is accompanied by possibly large numerical diffusion, which can artificially smear solutions 

in the presence of sharp gradients in velocity. 

While finite volume schemes can reduce to an equivalent finite difference scheme for certain 

configurations, conceptually finite volume schemes consider discrete shapes or volumes, for 

which a grid point typically lies in the center. Fluxes are calculated along the sides of the 

discrete shape, and these fluxes are balanced according to the set of differential equations 

solved. As neighboring volumes share boundaries, a flux leaving one volume is necessarily 

equal to the flux entering its neighbor, and so this approach has the advantageous property of 

being conservative; this is a property that is not guaranteed with finite difference models. 

Various techniques have been developed to evaluate boundary fluxes, and flux limiters can be 

integrated to increase model stability in the presence of flow discontinuities (see Toro, 1999 and 

LeVeque, 2002, for a complete review of the finite volume method). Tsunami models that use 

the finite volume approach include TsunamiClaw (Berger and LeVeque, 1998), COULWAVE

(Kim et al., 2009), and FUNWAVE (Shi et al., 2012). Of these three models, TsunamiClaw can 

utilize adaptive mesh refinement, while FUNWAVE and COULWAVE solve weakly dispersive

Boussinesq-type equations.

Finite element approaches permit arbitrary location of grid points or nodes, where the 

connections, or elements, between the points create the numerical mesh. Fundamentally, finite 

element methods can be considered the generalized forms of both finite difference and finite 

volume schemes. Due to their ability to handle arbitrary spatial distributions of nodes, finite 

element methods have the ability to capture a wide range of spatial resolutions with a single 

mesh. Finite element models that used for tsunami studies include ADCIRC (Dawson et al., 

2006) and SELFE (Zhang and Baptista, 2008).

Coupled and Hybrid Techniques

Coupled and nested methods may be desirable to examine transient wave phenomena or wave 

hydrodynamics on a fine scale. Of the two typical phase-resolving, depth-integrated models 
4-2 

used in wave studies, the Boussinesq can be considered a more physically complete, or at least 

a higher-order, approximation compared to the NSW. In coastal regions, where the water depth 

is very shallow and thus amplitude and wavelength become high and short, nonlinear and 

bathymetric interactions across a wide range of frequencies occur. These interactions can

locally generate various shorter-crested, or dispersive wave components even if the offshore 

forcing is considered a long wave. A well-known example is the transformation of a tsunami 

front into an undular bore (e.g. Matsuyama et al., 2007). Thus, if the wave evolution is expected 

to be nonlinear and (possibly) dispersive, the Boussinesq model is appropriate. However, the 

additional physics included in the Boussinesq approximation come with a substantial 

computational cost, often making the model impractical for ocean basin-scale simulations. 

Furthermore, to use the physical advantages of the Boussinesq model for a local region in the 

nearshore zone, it becomes necessary to couple the Boussinesq model with some other source 

of wave information for its boundary conditions. The obvious coupling choice would be the 

NSW, proven for both efficient and accurate basin-scale tsunami prediction (Son et al., 2011). 

There are numerous challenges with this coupling, most notably that the two approximations 

(NSW and Boussinesq) are different, creating a physical mismatch across the coupling 

interface. Also, the NSW will typically have a low-order numerical solution approach (hence its 

computational efficiency), while high-order partial derivations in the Boussinesq model require a 

high-order scheme. Matching these two schemes can also create numerical stability issues. 

Details regarding a coupling NSW-Boussinesq approach can be found in Son et al. (2011), 

which indicates that a coupling approach might be useful for site-specific tsunami hazard 

assessment in coastal regions. 

Moving Shoreline Algorithms

In order to simulate the flooding of dry land by a tsunami, a numerical model must be capable of 

allowing the shoreline to move in time. Here, the shoreline is defined as the spatial location 

where the solid bottom transitions from submerged to dry, and is a function of the two horizontal 

spatial coordinates and time. Numerical models generally require some type of special 

consideration and treatment to accurately include these moving boundaries; the logic and 

implementation behind this treatment is called a moving shoreline, or runup, algorithm.

For typical tsunami propagation models, it is possible to divide runup algorithms into two main 

approaches: those on a fixed grid and those on a Lagrangian, or transformed, domain. Both 

approaches have their advantages and disadvantages; currently, fixed-grid methods are more 

commonly found in operational-level models (e.g., Titov and Synolakis, 1998), likely due in large 

part to their conceptual simplicity. This section provides a review of these two classes of 

models, followed by a review of the standard analytical, experimental, and field benchmarks 

used to validate the runup models. Pedersen (2004) provides additional information and a more

comprehensive review.

With a fixed-grid method, the spatial locations of the numerical grid points or control volumes 

are determined at the start of a simulation and do not change shape or location throughout the 

simulation duration. These methods can be classified as extrapolation, stair-step, auxiliary 

shoreline point, and permeable beach techniques. The extrapolation method has its roots in 

Sielecki and Wurtele (1970) with extensions by Hibberd and Peregrine (1979), Kowalik and 

Murty (1993), and Lynett et al. (2002). The theory behind this method is that the shoreline 

location can be extrapolated using the nearest wet points, such that its position is not required 

to be locked onto a fixed grid point; it can move freely to any location. Theoretically, the 

extrapolation can be of any order; however, from stability constraints a linear extrapolation is 
4-3

generally found. Hidden in the extrapolation, the method is roughly equivalent to the use of loworder diffusive directional differences taken from the last wet point into the fluid domain (Lynett 

et al., 2002). Additionally, there are no explicit conservation constraints or physical boundary 

conditions prescribed at the shoreline, indicating that large local errors may result if the flow in 

the extrapolated region cannot be approximated as a linear slope. The extrapolation approach 

can be found in both NSW and Boussinesq models with finite difference, finite volume, and finite 

element solution schemes, and has shown to be accurate for a wide range of non-breaking, 

breaking, two horizontal dimension, and irregular topography problems (Lynett et al, 2003; 

Pedrozo-Acuña et al., 2006; Cienfuegos et al., 2007)

Stair-step moving shoreline methods, one of the more common approaches found in tsunami 

models (e.g., Liu et al., 1994), reconstruct the naturally continuous beach profile into a series of 

constant elevation segments connected through vertical transitions. Essentially, the bottom 

elevation is taken as the average value across a single cell width. A cell transitions from a dry 

cell to a wet cell when the water elevation in a neighboring cell exceeds the bottom elevation, 

and transitions from wet to dry when the local total water depth falls below some small threshold 

value. These methods are particularly useful in finite volume and staggered grid approaches 

(e.g., Arakawa and Lamb, 1977; LeVeque and George, 2004), but can be difficult to implement 

in centered difference models, particularly high-order models or those sensitive to fluid 

discontinuities, where the "shock" of opening and closing entire cells can lead to numerical

noise. Auxiliary shoreline point methods require dynamic re-gridding very near the shoreline, 

such that the last wet point is always located immediately at the shoreline. Obviously, this 

method requires a numerical scheme that can readily accommodate non-uniform and changing 

node locations. The moving shoreline point must be assigned some velocity, which can be

extrapolated from the neighboring wet points using some extrapolation method discussed 

above. However, it is fundamentally different in that the shoreline point is explicitly included in 

the fluid domain. Therefore, the governing conservation equations near the shoreline are more 

precisely satisfied here, although still dependent on the appropriateness of the extrapolation. 

One such method, that of Titov and Synolakis (1995), has been successfully applied in the Nonlinear Shallow Water (NLSW) equation models.

Another fixed-grid treatment of moving boundary problems is employing a slot or permeableseabed technique (Tao, 1983 and 1984). Conceptually, this method creates porous slots, or 

conduits, through the dry beach, such that there is always some fluid in a "dry" beach cell, 

although it may exist below the dry beach surface. These porous, "dry" nodes use a modified 

form of the NLSW. Although in concept this approach is modeling a porous beach, it is not 

attempting to simulate the groundwater flow under a real, sandy beach, for example. The 

equations governing the "dry" domain contain a number of empirical parameters tuned to 

provide reasonable runup agreement with benchmark datasets. The advantage of this 

approach is that it allows the entire domain, including the fluid and "dry" nodes, to be

determined via a somewhat consistent set of governing equations without requiring a direct 

search routine to determine the shoreline location. The method has gained some popularity in 

wind wave models (e.g., Madsen et al., 1997; Kennedy et al., 2000) where a highly accurate 

estimate of the shoreline location is not the highest priority. However, the approach has also 

been successfully used in tsunami studies (e.g., Ioualalen et al., 2007) despite the fact that the 

empirical coefficients that govern the model accuracy cannot be universally determined for a 

wide range of problems (Chen et al., 2000).

An alternative to fixed-grid methods is the Lagrangian approach. The Lagrangian approach 

discretizes the fluid domain into particles, or columns of fluid in depth-integrated models, that 

are transported following the total fluid derivative. There are no fixed spatial grid locations; the 
4-4

columns move freely in space and time, requiring numerical flexibility for constantly changing 

space and time steps. The Lagrangian approach is both a more physically consistent and 

mathematically elegant method of describing shoreline motion. The shoreline "particle" is 

included in the physical formulation just as any other point in the domain (i.e., no extrapolations 

are necessary), and thus the shoreline position accuracy will be compromised only by the 

overarching physical approximation (e.g., long wave approximation) and the numerical solution 

scheme (e.g., second-order time integration). However, with the increased accuracy, this 

mathematical system can be more difficult and tedious to solve numerically, typically requiring 

domain transformation, mapping, and/or re-gridding. Lagrangian methods have been used 

successfully in finite difference and finite element NSW and Boussinesq equation models (e.g., 

Pedersen and Gjevik, 1983; Gopalakrishnan and Tung, 1983; Petera and Nasshei, 1996; Zelt, 

1991; Ozkan-Haller and Kirby, 1997; Birknes and Pedersen, 2006).
5-1 

5. APPLICATION OF NUMERICAL MODELS

Specification of Initial Condition

There are two general categories of initial tsunami conditions in numerical models: hot-start and 

cold-start. In a hot-start model, there is some profile of the free surface elevation (and/or 

velocity) forced at the starting time of the simulation. This profile must be taken from some 

other model. Alternatively, during a cold-start, all values of surface elevation and velocity are 

zero and any waves must be generated by either a lateral bottom (e.g., moving bottom due to 

landslide) or free surface (e.g., atmospheric pressure gradient) boundary condition.

Hot-start conditions are often found in earthquake-generated tsunami studies. Here, it is 

assumed that the seafloor deformation is impulsive (occurring over a time scale much shorter 

than the wave period) such that the initial free surface profile matches this deformation profile 

exactly, and the velocity everywhere is zero. For earthquake-generated tsunamis, the theory of 

Okada (1975) generates the initial condition with information about the earthquake location, 

geometry, displacement, focal depth, and dip, slip, and rake angles. Section 7.1 of this report 

presents the options for generating a hot-start initial condition for landslide-generated tsunamis. 

Care must be taken when using hot-start conditions with models that solve equations with highorder derivatives. For example, time integration of the highly-nonlinear Boussinesq-type model 

will require information at three previous time levels to solve the current time level. This 

requires that free surface (and velocity) profiles at three previous times are defined to hot-start 

the model. For an earthquake-generated tsunami, free surface profiles at these three previous 

time levels are not known, and a reasonable solution to this problem is to use the same initial 

surface profile at the three different time levels. However, this is problematic because such a 

specification implies, or enforces, a certain value of the time derivative of the free surface 

elevation for a situation in which these derivatives are not clearly defined. A practical 

workaround to this issue is to use a low-order time integration of the shallow water equations 

initially and then transition to the full, high-order solution after the third time step. 

Dynamic Bottom Boundary Conditions

Dynamic bottom boundary conditions are often used in conjunction with a cold-start simulation 

to model the waves generated by some type of bottom motion. In shallow water equation 

models, the moving bottom is accommodated by an additional forcing term in the continuity 

equation equal to the time rate of change of the bottom elevation. Some Boussinesq-type

models include additional terms in the conservation of momentum equations (e.g. Lynett and 

Liu, 2002) to account for the motion of relatively short segments of seafloor.

For earthquake-generated tsunamis, dynamic bottom boundary conditions are sometimes used 

for earthquakes with a long rupture duration; where “long” implies a duration greater than 1/10 

of the generated wave period. In these cases, the waves generated during the initial rupture 

may have traveled a non-negligible distance from the source region by the time the rupture has 

completed and this effect can play a role in the along-front tsunami energy distribution. 

Information regarding the rupture timing of individual subfaults can be extracted from finite fault 

solutions.

For landslide-generated tsunamis, specification of the seafloor time-history is a challenging task. 

Figure 5-1 shows the wide variety of forms that submarine mass movements can take a wide 

variety of forms, each with a different tsunamigenic efficiency. The tsunamigenic efficiency is a 
5-2 

measure of the transfer of the energy of the mass movement into free ocean surface waves. In 

general, mass movements that are "shallow," meaning the horizontal length scale of the slide is 

much greater than the local depth, "coherent," meaning that the mass fails as a single piece, not 

as a group of smaller, spatially and temporally separated segments, and "fast," meaning that the 

time scale of motion of the movement is on the order of the generated wave period, tend to be 

more efficient in transferring their energy to the water column. However, it is difficult to know 

beforehand how a landslide will fail, which implies that a deterministic analysis should use a 

conservative assumption about this efficiency. Section 7.1 provides further discussion on this 

topic. Note the above discussion details efficiency and does not address tsunami potential, 

which is a function of efficiency but also closely related to landslide volume. 

Once a time history of seafloor motion or rheological slide model has been chosen, this 

information is loaded into the hydrodynamic model and the waves are generated, as discussed 

earlier in this section. In depth-integrated models, the evolving seafloor must remain singlevalued in the vertical (i.e., a horizontal field). Furthermore, steep, moving segments are akin to 

short waves and are, therefore, not as accurate. Generally, Lynett and Liu (2002) note that in 

shallow water models, in order for the wave generation by a dynamic bottom boundary condition 

to be accurate, the "characteristic" horizontal length-scale of the slide divided by the local water 

depth should be no less than 25, but suggest a less restrictive lower limit of 10 for typical 

Boussinesq-type models. However, there can be numerous important length scales of an 

irregular and evolving mass; such slides should be examined spectrally to understand at which

length scales there is significant deformation. The primary conclusion for practical application is 

that it can be very difficult to satisfy the above accuracy restrictions for an arbitrarily deforming 

mass and/or a rheology that tends to develop steep fronts. In vertically resolving, fully nonhydrostatic models, these limitations do not exist and, theoretically, any arbitrary seafloor motion 

can be accommodated. Additional prescriptive details on choosing initial tsunami conditions 

can be found in Section 7.1.

Topographic Grid Creation

There are many sources of bathymetric and topographical data. Data at resolutions from 30 

arcseconds to 2 arcminutes are available across the globe (e.g., GEBCO) and are acceptable 

for use on open ocean propagation of long waves. For detailed, site-specific coastal simulation, 

higher grid resolution is often required. In many locations along the U.S. coast, gridded 1/3 

arcsecond (approximately 10 m) bathymetry and topography data are available through NOAA’s 

Tsunami Inundation Digital Elevation Models database. Additional publically available 

resources for coastal bathymetry and topography include NOAA’s Digital Coast database, which 

collects chart and recent LiDAR data. Furthermore, many states host their own online 

bathymetry databases, although most of these datasets are also found on the NOAA sites. 

Care must be taken to ensure datum consistency between different topographic datasets. In 

locales where no reasonable digital bathymetry data exist, navigation charts of the area can 

provide useful information, however these data are often outdated and coarse. In the absence 

of even nautical charts, caution must be taken when using GEBCO and SRTM data sources in 

coastal regions, as these often exhibit large errors for water depths less than 50 m. For sitespecific modeling of tsunami effects, topographic and protective features, such as seawalls and 

levees, can be included in the grid as they can have important effects on the local flow (e.g. 

Suppasri et al., 2013). However, such features must be shown to remain intact during the PMT 

event in order to be included in the grid.
5-3

Figure 5.1 Various types of movements, as a function of material class (taken from 

Cruden and Varnes, 1996.) 
5-4 

Grid and Time Steps

In most numerical models, the time step is controlled by a maximum allowable Courant number

(Fletcher, 1991), and is therefore connected to the grid size used. While topographic grid 

resolution is controlled by the available data sources, the numerical grid resolution is controlled 

by the wavelengths that need to be properly resolved. In general, for a reasonable numerical 

representation of the waveform, there should be at least 10 grid points covered in a single 

wavelength. While an apparently simple statement, in practice this can be challenging due to 

the changing length of the wave through shoaling and nonlinear generation of higher harmonics 

(shorter waves) in shallow water. Furthermore, certain "thin" coastal features, such as 

breakwaters, can be very important to the local tsunami effects and must be resolved. In 

practice, a horizontal grid spacing of 10 m is adequate for proper resolution of detailed tsunami 

evolution for a site-specific coastal/onshore modeling application. Numerical convergence tests 

(e.g., comparing 10 m and 15 m resolution modeling results) are often undertaken to 

demonstrate acceptable resolution of hydrodynamic processes. While 10 m resolution is 

adequate for coastal and onshore resolution, in the deep ocean grid lengths of a few kilometers 

may provide proper resolution, indicating the usefulness of variable gridding. Use of nested or 

telescoping grids is a common variable resolution option (e.g. Titov et al., 2011) in tsunami 

modeling

Bottom Friction and Turbulence Closure

In most shallow water models, the dominant physical turbulence dissipation mechanism is 

bottom friction. This dissipation appears as a drag-like momentum flux term, where the 

magnitude of the flux is controlled by an empirical coefficient. There are numerous formulations 

of this coefficient, such as Mannings, Chezy, and Moody. All require some information about 

the bottom type, including roughness height or character. The Manning’s formulation is found 

commonly in U.S. projects, and is controlled by bottom type, such as smooth, grassy, or rocky. 

Naturally, this bottom type changes throughout the domain, and physically so should the 

Manning’s coefficient. Accurate specification of the bottom friction coefficient is particularly 

important for overland flow and for propagation over wide, shallow continental shelves. Section 

7.1 provides guidance on the selection of bottom friction coefficients. 

Boussinesq-type models also attempt to capture breaking dissipation through ad-hoc breaking 

models. These models include numerous empirical coefficients, which are tuned with laboratory 

and field data; agreement for break point, surf zone wave height, and wave train 

asymmetry/skewness is quite good (e.g., Kennedy et al, 2000) despite the highly empirical 

nature of the breaking schemes. Additional dissipation mechanisms include sub-grid eddy 

viscosity (Chen et al., 2003) and vertical mixing models (Kim and Lynett, 2012).

Antecedent Water Level

NRC Japan Lessons Learned Project Directorate (JLD) Interim Staff Guidance (ISG) JLD-ISG2012-06 provides guidance for estimating the high antecedent water level to be used with sitespecific tsunami modeling. This water level consists of the 10% exceedance high tide, the initial 

rise (or sea level anomaly), and water level increase due to site-specific historical sea level rise 

over the design life. Regulatory Guide (RG) 1.59 and ANSI/ANS-2.8-1992 define the 10 percent 

exceedance high spring tide as the high tide level that is equaled or exceeded by 10 percent of 

the maximum monthly tides over a continuous 21-year period. It is remarked that this is a rare 

water level; if the crest of a high tide persists for one hour, then this tidal elevation is exceeded 

only once every 7200 hours (assuming 30 days in a month). Furthermore, if the peak flow 
5-5 

depth considered during a tsunami event lasts for one hour, the likelihood of the maximum 

tsunami elevation occurring at the same time as the 10% monthly tidal level defined in RG 1.59 

is likewise 1-in-7200 (occurs 0.014% of the time). It may be reasonable for an analysis to use a 

less conservative tidal level, if that analysis includes some type of sensitivity analysis or 

stochastic component wherein this type of probability might be consistently incorporated. 

Associated Effects

JLD-ISG-2012-06 discusses the associated effects, such as hydrostatic and hydrodynamic 

forces, debris and water-borne projectiles, and sediment erosion and deposition. Detailed 

discussion of these effects is out of scope of this document. Detailed simulation tools exist to 

quantify some of these processed (e.g. Son and Lynett, 2014) but such tools are in an earlier 

stage of development and are not yet widely used. Prescriptive approaches to estimate these 

effects are available in the Coastal Engineering Manual (CEM).

6-1 

6. INTERPRETATION AND ANALYSIS OF NUMERICAL RESULTS

Numerical Convergence

As mentioned in the Grid and Time Steps discussion in Section 5.4, a test to demonstrate that 

the hydrodynamic predictions of interest are not dependent on the grid resolution should be 

performed. While there are numerous methods to demonstrate this property, a reasonable 

approach would be to increase and/or decrease the target resolution by 25% and compare 

results. Here, the target resolution is the resolution of expected convergence, for example 10 m 

for a detailed, coastal/onshore simulation. If the output of interest is the maximum water surface 

elevation and current speed near a particular location, these values as predicted by the 

simulations with different resolution should be presented and compared. If the results are within 

a reasonable tolerance, here 5% of the target resolution result, it can be stated that the target 

resolution is convergent and numerically acceptable.

Types of Output

Numerical simulations provide a wealth of output that can be presented in different ways. The

focus in this report is on output types typically required for the presentation of tsunami hazard 

analysis.

6.2.1 Instantaneous Snapshots of Ocean Elevation and Speed

Snapshots of a physical solution parameter, such as ocean surface elevation or velocity, display 

the parameter in two or three spatial dimensions at a particular instant in time. A series of these 

snapshots is useful in graphically presenting the spatial evolution of the tsunami in a technical 

report, whereas an animation of the parameter might be used in an oral presentation. Figures 

6.1 and 6.2 are snapshot image examples. Figure 6.1 provides a snapshot of the tsunami on a 

transatlantic propagation grid, while Figure 6.2 gives the local, near-coastal ocean elevation. 

Note that in both figures, the time of the snapshot is given, and axes are labeled with 

appropriate spatial dimensions. 

Figure 6.1 Snapshot of the ocean surface elevation as the tsunami travels across the 

Atlantic
6-2

6.2.2 Time Series

Time series are used to present the temporal evolution of a simulation parameter at a single 

location. These images are useful for presenting ocean surface elevation time histories near 

the site or at some offshore location, for example the location at the offshore limit of the highresolution coastal inundation grid. This is often the product used to validate modeling results 

with real events, since usually a time series of water elevation during tsunami events can be 

obtained from coastal tide gauges.

6.2.3 Maximum and Minimum Surface Plots 

Maximum surface plots provide a summary graphic of an entire simulation. These images are 

useful to present when discussing the maximum levels reached near a particular location. 

Figure 6-3 shows the maximum ocean surface and fluid speed recorded during a simulation.

Figure 6.2 Snapshot of ocean surface elevation in the nearshore, high-resolution grid
6-3

6.2.4 Tabular Summary of Results

In addition to the graphical products described above, a hazard analysis should also include a 

summary table of the site-specific results. Such a table should include tsunami properties such 

as maximum flow elevation, minimum flow elevation, and maximum flow speed. These 

properties should be listed for all tsunami sources included in the site-specific analysis, for all 

relevant locations at the site (e.g. intake, outfall, reactor building, etc.).

Figure 6.3 Maximum elevation (left) and speed (right) predicted by a high resolution 

simulation

7-1 

7. EVALUATING TSUNAMI HAZARDS NEAR CRITICAL FACILITIES

This section describes the reasonable approaches for the determination of tsunami hazard at 

NPPs. These approaches are not requirements for analysis; any state-of-the-art method

consistent with the guidance provided in General Design Criterion 2 and discussed in Section 1 

would be sufficient. The information presented below builds closely off Section 8 in PMEL-136, 

"Template THA." Readers should be familiar with this document before performing a THA.

Deterministic Approach

The deterministic approach should employ a "hierarchal" method, where extremely high levels 

of conservatism are initially employed in the tsunami hazard analysis. If, after the initial 

analysis, it is shown that the PMT might impact the site, various levels of conservatism can be 

removed with careful justification. This method allows the analysis to show how the different 

conservative assumptions affect the final PMT assessment. Though not necessarily efficient, 

this approach is systematic and clear, which are highly important in justifying a PMT evaluation 

approach. A number of the conservative approaches recommended here may appear so 

unrealistic as to be unfair to the analysis. However, this is a necessary side effect of a 

deterministic approach for a poorly understood hazard that should have a low annual recurrence 

frequency.

For sites that require a PMT study, the first step is to provide a complete tsunami source 

characterization considering all local and distant sources that might affect the site. Information 

about historical sources can be found in tsunami catalogs, and supplementary materials in 

technical journals, particularly for landslide information. For the deterministic approach, every 

allowance should be made to describe the source in the most conservative, but physically 

reasonable, manner. For earthquakes, the maximum slip should be based on theoretical 

considerations of the maximum possible earthquake, not based on historical precedent of any 

particular fault or fault segment. Hypothetical ruptures should be allowed to pass through 

proposed or speculated rupture barriers, such as bends or transitions, to permit extremely long 

rupture lengths. For example, a tsunami source from the Caribbean Subduction Zone should 

include rupture along the entire subducting area (Figure 7-1), with slip based on the corner 

magnitude in a Gutenberg-Richter magnitude-frequency relationship for all global subduction 

zones. While such a distribution may not be valid for the local subduction zone, the 

conservative and statistically defensible option is that this distribution applies to all subduction 

zones. A deterministic approach, void of any measure of uncertainty quantification, requires 

such a conservative procedure.
7-2

Large earthquakes are considered the most common tsunami source because they have the 

potential for generating large amplitude waves of long enough wavelengths to travel thousands 

of miles from the source region to affect remote coastlines. For many coastal locations, even 

for those without a nearby tsunami source region, these seismically generated tsunamis can be 

the source for the PMT. However, at other locations, local landslides capable of generating very 

large wave heights over a small area can more commonly represent the PMT source. At any 

location, both types of tsunamis should be evaluated for their potential to represent the PMT. In 

the absence of any statistical information regarding the location of past and future landslides, it 

should be assumed that any local landslides that have occurred in similar geophysical 

conditions as found near the NPP may in the future occur immediately offshore of the NPP. For 

example, if there is geological evidence of a submarine landslide along the shelf break at a 

location distant to the NPP but in a similar geophysical configuration, then that landslide should 

be considered to possibly occur at the shelf break immediately offshore of the NPP. This 

Figure 7.1 Initial tsunami condition for an earthquake along the Caribbean Subduction 

Zone using the deterministic approach 
7-3

approach can be discarded only with a clear and convincing geological argument against 

permitting a change in the future landslide to a location closer to the NPP. Due to an overall 

lack in understanding of the evolution of large submarine landslides and the waves that they 

generate, a highly conservative approach must be employed to initially estimate the waves 

generated. An initial estimation method is to assume that the vertical change in seafloor 

elevation due to the postulated landslide is reflected exactly in the water surface elevation.

Figure 7-2 provides a graphical example of this procedure. Landslide dimensions, including 

width, length, and depth should first be estimated from available geophysical data. With this 

information, and the local bathymetry data, it is possible to construct a "before" and "after" 

landslide profile. The initial free surface profile is determined from the difference in these two 

profiles. Fundamentally, this method is assuming an impulsive vertical motion of the landslide,

which is conservatively defensible in light of the current state of knowledge for tsunami 

generation by landslide, and eliminates some of the significant uncertainty inherent in using 

landside-tsunami initial condition functions found in the literature (e.g. Watts et al., 2003). This 

"impulsive" landslide option can be discarded for alternative generation descriptions, such as 

modeling the time-history of the slide, only with clear and convincing arguments that the 

alternative approach is conservative as well as a physically reasonable estimation (e.g. Geist et 

al., 2009).

Figure 7.2 Initial condition for a tsunami generated by a submarine landslide. The 

"before" landslide profile is given by the solid yellow surface and the 

"after" landslide profile is given by the black dashed line. The difference 

in these two profiles results in the initial sea surface elevation profile as 

shown in solid blue.
7-4

If it is important or necessary to model the time history of the slide and the waves generated 

from it, the following general procedure provides a methodology. First, the geometry of the 

mass movement must be specified. Both volume and the two-dimensional slip surface must be 

prescribed; note that this information is also required in the above "impulsive" landslide option. 

Next, either a mass movement failure mechanism (e.g., rotational, translational, types shown in

Figure 5-1) or a slide material behavior model (e.g., solid body motion, Bingham plastic flow

(Bingham, 1916)) should be selected. If the failure mechanism is specified, then some 

analytical model to provide the slide motion must be provided. If the material behavior model is 

used, then the mathematical governing equations and the numerical solution technique, if 

applicable, must be given. Regardless of the method used to provide a slide time-history, all 

empirical coefficients that govern said the selected method should be clearly specified and 

explained. In particular, a discussion regarding the physical ranges of these coefficients and 

their effect on the resulting solution should be provided. Furthermore, justification that the 

chosen value for each coefficient is appropriate for the slide of study and conservative with 

respect to the resulting site-specific tsunami impacts should be provided. Such a demonstration 

might be accomplished through a sensitivity study wherein numerous simulations covering the 

potential variation of the coefficients controlling the slide motion are compared.

With a specification of the initial tsunami condition or a slide time-history completed, modeling of 

the wave propagation can commence. The choice of the model will depend on the properties of 

the generated waves as well as the expected transformation of the waves through shallow 

water. For example, if there are substantial dispersive frequencies in the source condition, then 

a dispersive model should be used (Lynett and Liu, 2002). This is most likely to be the case for 

landslide-generated tsunamis, although is not always the case (e.g. Kirby et al., 2013). If a 

dispersive model is chosen, care must be taken to ensure that the grid resolution utilized by the 

model is capable of resolving the dispersive components in the source. Reasonable arguments 

can often be made to show that dispersive effects should not be important along the 

propagation path, or that a non-dispersive model will provide a conservative result. However, in 

light of the complex wave transformations that might occur for the types of large waves 

associated with a PMT, as discussed in Section 2, a nonlinear-dispersive model (e.g. Wei et al., 

1995) is generally preferred. Regardless of the level of dispersion included in the modeling 

approach, nonlinearity in the wave evolution is likely to always be an important effect as the 

wave approaches very shallow water.

Site-specific simulation of tsunami impacts using one-horizontal-dimension (1HD) transects is 

the first level of analysis in the hierarchal approach. The most conservative bathymetry and 

topography configuration is a 1HD transect that is devoid of bottom friction. Use of a 1HD 

transect for a local landslide tsunami simulation is equivalent to assuming that the entire 

alongshore length of shelf break fails as a large, simultaneous landslide. There are no twodimensional spreading effects and removing bottom friction eliminates any damping effect on 

the wave. A 1HD transect with no friction provides an upper limit of the possible flow inundation 

extent. While removing bottom friction is unphysical (even perfectly smooth surfaces will add 

non-negligible bottom friction dissipation over long lengths of propagation), specification of an 

upper limit effect is useful for safety assessment. If, under these assumptions, the tsunami 

does not reach the NPP, the source under investigation can be ignored as it is no longer a 

justified PMT candidate. 

If the tsunami does reach the NPP site, the next step is to investigate physically reasonable 

values of bottom friction for the particular site. An approach would be to employ a conservative 

value for overland flow while maintaining no friction in offshore areas. For a site that is fronted 

by thick vegetation and brush, a conservative friction model might use a Mannings nn value of 
7-5 

0.025. This is a value more appropriate for smooth surfaces and implicitly assumes that the 

land use of the area may change in time from a natural bottom to an artificially smooth bottom.

While the 1HD simulations provide a conservative measure of the possible inundation extent, 

site-specific, refractive focusing may possibly lead to local amplifications in wave height All 

sources that reach the site with conservative friction values from the 1HD simulations should be 

performed in full 2HD, which is the second level of analysis in the hierarchal approach. If no 

source waves reach the site with realistic friction, then the two sources with the largest 1HDconservative-friction wave height near the NPP should be tested in full 2HD. Alternatively, a 

user may skip the 1HD screening of sources and examine all sources under consideration that

impact coastal areas near the site in full 2HD. Section 9 presents an efficient alternative 

approach to the 1HD computations for seismically generated tsunamis. 2HD simulations should 

be run with conservative friction to account for the lack of knowledge of how the area seaward 

of the NPP might transform in the future. Note that if the site-specific grid does not resolve the 

structures onsite, then, if the tsunami does reach the site, the maximum elevation predicted by 

the modeling should be increased to include runup on vertical walls. The resulting largest 

tsunami elevation at the NPP represents the PMT. A flow chart with the deterministic PMT 

evaluation logic is given in Figure 7-3. It would be reasonable to interpret this figure as a 

detailed expansion of the "Analyses" box in Figure 8-4 in PMEL-136, where the analysis utilizes 

a deterministic method.

Probabilistic Approaches

For tsunami hazard assessment working towards the goal of risk assessment, Probabilistic 

Tsunami Hazard Assessment (PTHA) will in general provide more useful and appropriate 

results as compared to a deterministic analysis. This section highlights the challenges in using 

PTHA for landslide sources and the goal of determining the PMT, and is not meant to be a 

primer on PTHA for earthquake sources. Additional technical details on PTHA for earthquake 

sources is provided in Section 10. The following excerpt, taken from PMEL-136, is still largely 

valid:

In principle, an effective PTHA methodology would provide a more realistic and 

scientifically rigorous framework for decision-making during NRC reviews of NPP 

applications, since such reviews would be based on quantitative hazard level 

estimates -i.e., the probability of occurrence for tsunami events with an estimated 

level of destructive potential. The PTHA methodology [for earthquake-generated 

tsunamis] has matured to the point of prototype application at Seaside, Oregon, 

and is being considered for adoption by Federal agencies and U.S. States 

charged with THA-related missions.

With the completion and the publication of the Seaside Study (González et al., 2009), the 

largest remaining hurdle in the application of PTHA for the evaluation of the PMT is 

incorporating landslide-generated tsunamis into the analysis, or more specifically, incorporating

local sources with significant heterogeneity in their spatial geometry and time evolution. This is

classified as epistemic uncertainty, driven by a lack of understanding of how small-scale spatial 

and temporal details generate and evolve. If this heterogeneity can be quantified, as it can be 

to a reasonable extent for earthquake sources, proceeding with a PTHA is justified. However, 

for landslides, this is a challenge. Furthermore, the average return periods of submarine 

landslides in a given area, a piece of information integral to a PTHA analysis, are often 

accompanied by leading order imprecision and uncertainty. Figure 7.4 provides an example 

methodology for integration of landslides in PTHA, focusing only on the modeling of the 
7-6

landslide time history. The first piece of information required is a distribution that relates slide 

volume and horizontal slide area to return period but slide area is equally important, as without 

this piece of information it is not possible to estimate the slide thickness, which largely controls 

the generated wave height. Alternatively, a distribution that relates slide volume to slide area 

(and/or maximum slide thickness), when used in conjunction with a slide volume to return period 

relationship, would close this problem. Currently, limited information exists for building these 

starting-point distributions.

These distributions should be sampled to arrive at a set of landslide volume plus thickness 

combinations with a unique time-history determined for each. The first step in this procedure is 

for the user to decide how he/she will describe the motion. There are two main categories of 

choices (which could be combined into a single category if desired): (1) specify the failure type 

or (2) specify the rheological model. In the "Specify Failure Type" route, following Figure 5-1, a 

set of different types of failure mechanisms, such as rotational slides, translational slides, etc., 

are included. Each of these different failure mechanisms will yield different tsunamigenic 

efficiencies as well as different generated wave properties. This step represents a branching 

logic tree and a weighting factor must be assigned for each different failure mechanism. The 

variety and number of different failure mechanisms can be reduced based on site-specific 

conditions, and the weighting factors must be location-specific as well. Once a failure 

mechanism is chosen, the parameters that govern its evolution must be selected. For example, 

if working with a translational landslide, the initial motion would likely follow solid body motion;

therefore, parameters such as drag coefficient, added mass, as well as material density must be 

given. These parameters should be selected from distributions created for each parameter; 

such information could be generated, but is not known to formally exist currently for submarine 

landslides. It is also worth noting that while the evolution model for a translational slide (solid 

body motion) is established in the literature (e.g. Watts et al., 2003), the same cannot be stated 

for any of the other motions associated with the failure mechanisms. With the parameters 

necessary for the failure model provided, the slide time history can be imported into the 

hydrodynamic model, and a single realization can be generated. In the Monte-Carlo analysis, 

this procedure is repeated until enough realizations are generated so that distributions can be 

constructed for a tsunami impact metric under analysis. Instead of choosing the failure 

mechanism route, it is possible to specify a set of different rheological models. Similar to the 

failure mechanism approach, the range and weightings of the employed rheological models is 

location specific. Each rheological model contains a handful of material and empirical 

parameters, and a distribution for each should be constructed.

The methodology outlined in Figure 7.4 and the previous paragraphs is one suggestion for a 

probabilistic analysis incorporating landslides - an analysis procedure that is currently not well 

defined. Alternative methods are viable, and significant simplifications may be made, with 

justification, to the generality given here. These simplifications will often exist in the description 

of the failure mechanism/rheological models and distributions of the parameters that govern 

them. The current knowledge of wave and landslide coupling may not justify a confident 

simplification for these models/parameters, hence the current preference for a conservative 

source estimate with the deterministic approach. Likewise, the computational cost for a PTHAwith-landslides, driven by the high computational cost of 3D/dispersive models often needed for 

accurate physical representation of the slide and waves, is significant and may only be justified 

for sites where the PMT is the design basis for one or more SSC.
7-7

Figure 7.3 Procedural flowchart for a deterministic evaluation of the PMT
7-8

Figure 7.4 Example procedure tree for the inclusion of slide variability and uncertainty 

into a PTHA
8-1

8. PRESENTATION OF TSUNAMI HAZARD ANALYSIS

The presentation of the PMT analysis is flexible, but should include at minimum:

• An overview image and accompanying discussion showing the locations of all the

sources to be discussed, as well as the location of the site to be studied. The image 

should include annotation for all shoreline and bathymetric features that are mentioned 

in text that have relevance to tsunami propagation, dispersal, focusing, etc. The image 

should also include an outline on basin scale figure that shows domain of site figures or 

other figures pertinent map. Figure 8-1 shows an example of what these figures may 

look like.

• An image, or series of images, with units and vertical datum noted, and accompanying

discussion showing the detailed local bathymetry near the site to be studied, as shown in 

sample Figure 8-2.

Figure 8.1 Example of source location map including shoreline and bathymetric 

features
8-2

• Physical detail about all the sources to be modeled, including a snapshot of the initial 

surface conditions. For 1HD transect analysis, provide an image showing the location of 

the transect as shown in sample Figure 8-3.

Figure 8.2 Sample figure showing bathymetry local to the site of interest
8-3

• For each source, provide a series of no less than four spatial snapshots of the evolving 

tsunami at different times. At least one of the snapshots should be at the time of 

maximum water elevation at the site.

• For each simulation, provide a surface of the maximum water level, the maximum 

current speed, and the minimum water level (if appropriate).

• Present a table of the maximum water level (or seaward runup elevation), maximum 

speed and minimum water level (if appropriate) at the site location.

• Time series at pre-defined locations of interest to the NPP.

Discussion and analysis centered around these pieces of data will adequately present the

modeling results and PMT specification.

Figure 8.3 Sample figure showing initial surface conditions and location of the 

transect AAAA′

9-1 

9. A PRACTICAL APPROACH TO DETERMINING PMT

In the spirit of the assessment approach put forth in Section 0, the following sections propose a 

comprehensive methodology conducive to the identification of a seismically-generated PMT for 

specific coastal locations by providing a systematic method for the selection of a credible worstcase scenario source. The approach introduced here advocates the use of a pre-computed 

tsunami database developed in recent years at the NOAA Center for Tsunami Research 

(NCTR), and recently used in tsunami modeling studies conducted at the NCTR for the NRC 

(Titov et al.,2016). This approach encourages the use of tsunami forecast tools developed by 

NOAA for real-time tsunami forecasting, currently in operation at both of NOAA’s Tsunami 

Warning Centers. The tsunami hazard assessment guidance provided in this is document is 

consistent with that outlined in JLD-ISG-2012-06 (U.S. NRC, 2-012). The JLD document 

recommends a general hierarchical approach, where sources are screened-out and levels of 

conservatism are sequentially removed if needed, identical to that provided in detail here. Other 

aspects of a tsunami assessment for a NPP, including calculation of antecedent water levels, 

bathymetry, and choice of numerical model are consistent between JLD and this document; 

effectively, this document provides a high level of detail to the guidance outlined in JLD.

Background

As stated previously in this report and in PMEL-136, one of the major difficulties in assessing 

the exposure of a coastal location to tsunami damage is the short historical record associated 

with events of this type. This is particularly true along U.S. coastlines where, at best, historical 

records date back a few hundred years forcing scientists to resort to paleo-tsunami and 

sedimentology studies (Atwater, 2005) to evaluate the exposure and recurrence period of 

tsunamis along our seaboard in an attempt to increase the length of the historical record. 

However, due to the long recurrence periods that devastating tsunamis can have, even 

additional information provided by paleo-tsunami science has proven to be insufficient in 

determining the level of exposure of a particular site as evidenced by the tsunami disasters of 

Sumatra in 2004 and Japan in 2011. Two major conclusions can be extracted from these 

events in connection to THA. First, the evaluation of a PMT should not be based solely on 

evidence of historical events, even though such information should undoubtedly be included in 

the analysis. Hazard assessment studies based only on historical evidence of past events to 

determine a credible worst-case scenario are likely to under-predict such cases, as 

demonstrated by recent events. Second, limitations on the magnitude of potential tsunami 

sources based on geophysical analysis of local plate tectonics, such as that proposed by the 

Japan Society of Civil Engineers (JSCE), the Nuclear Civil Engineering Committee, and the 

Tsunami Evaluation Subcommittee in Chapter 4 of their report Tsunami Assessment Method for 

Nuclear Power Plants in Japan (JSCE, 2002), have led to underestimations of the risk. This 

was seen during the 2011 Fukushima-Daiichi NPP disaster. A recent report by the American 

Nuclear Society Special Committee on the Fukushima Daiichi accident American Nuclear 

Society, 2012) states: 

"The tsunami design bases for the Fukushima NPPs were not consistent with the 

level of protection required for NPPs. If the return period for a tsunami of the 

magnitude experienced in Japan is as short as reported (once every 1000 years), 

a risk-informed regulatory approach would have identified the existing design 

bases as inadequate. In light of the March 2011 event, the tsunami design bases 

for the Fukushima NPPs were clearly inadequate. The Jogan tsunami in 869 AD 

and related tsunami stones in Iwake may seem to suggest a return period of 
9-2 

1000 years for such large tsunami waves. However, TEPCO’s analysis of the 

Jogan event prior to March 2011 predicted a resulting inundation height within 

the Fukushima NPP design basis. The discrepancy emphasizes a need for a 

coherent regulatory approach that is risk-informed, utilizing the most advanced 

evaluation methodologies, accounting for all relevant data available, and 

employing robust mitigation features for beyond-design-basis occurrences." 

Following the recommendations of conservatism in source selection advocated in the present 

report and in light of recent tsunami events, the methodology presented here suggests a 

relaxation of the structural tectonic constraints limiting the magnitude of PMT sources, inferred 

from seismic analysis of the local plate tectonics, unless compelling evidence to the contrary 

exists. This relaxation will result in hazard assessment studies including a large number of 

potential sources that might be considered non-credible under more stringent seismic analysis. 

One of the most efficient ways of conducting the type of hazard assessment study proposed in 

this section is via the combined use of NOAA’s deep-water tsunami propagation database and 

site-specific inundation forecast models integrated in the tsunami forecast software called Shortterm Inundation and Forecasting of Tsunamis (SIFT). The following sections describe NCTR’s 

tsunami forecasting approach as well as the tools developed by NCTR for the application of its 

forecasting methodology. The last section in this chapter explains how these tools can be used 

to conduct THA studies that not only include the reproduction of historical tsunami events, but 

also evaluate tsunamis generated by sources that would have not been considered by other 

methods.

In order to understand the use of these forecasting tools in hazard assessment studies, NCTR’s 

forecasting methodology will be presented next. 

NCTR’s Forecasting Methodology

NCTR’s tsunami forecasting methodology was developed to incorporate real-time tsunami 

simulations into the forecasting process. Real-time numerical simulations are produced during 

the propagation stage of a developing tsunami with the intent of obtaining an estimate of the 

tsunami’s impact on the coastline before it makes landfall. If generation, propagation, and 

inundation are considered the three main development stages of a tsunami, real-time 

simulations would be performed during the propagation stage, as early as possible after the 

generation phase. The following introduces NCTR’s forecast methodology and describes the 

tools designed for its implementation. The use of these tools in hazard assessment and an

explanation of how some of the products generated during a site-specific hazard assessment 

study could be used during a tsunami event to generate a forecast for the site are also 

presented. While the generation of a real-time tsunami forecast will not affect the hazard 

assessment for the site, the availability of event-specific arrival times, estimated inundation 

extent, or lack thereof and expected wave amplitudes and currents are data that can be of vital 

guidance in the determination of the adequate operational level of the NPP during a tsunami 

event.

Due to the operational need of obtaining simulated forecast results in a timely manner before

tsunami arrival at the coastline, NCTR’s methodology accelerates the generation of a forecast 

by taking advantage of the linear behavior of tsunami waves in deep water. This shows that the 

linearity of long waves is determined by what is known as the linearity parameter (Johnson, 

1997), that is, the ratio of the wave amplitude to the water depth, εε. For small values of εε, wave

amplitudes are small compared to the water depth, and current speed values
9-3 

(uu = HH/2�gg/h cos (kx − σt)), as described in Section 3.1, are small compared to the velocity of 

wave propagation, �ggh. Under these circumstances, the behavior of the wave can be 

considered linear. Typical tsunami wave amplitudes observed in the deep ocean range from a 

few centimeters to a couple of meters. These values are several orders of magnitude smaller 

than the average depth (h) of the Pacific Ocean (approximately 4500 meters). Accordingly, 

deep-water current speed values are estimated at a few centimeters per second, much smaller 

than the speed of wave propagation given by the expression �ggh. In light of this evidence, 

tsunami behavior in deep water can be considered linear. The linearity or non-linearity of the 

wave in shallow coastal waters is less straight forward. In this case the relative size of the 

amplitude of tsunami waves must be compared to the local depth, if wave amplitude is at least 

an order of magnitude smaller, wave behavior will be dominated by linear dynamics, if not, nonlinear effects may no longer be negligible. In very shallow waters of less than 10 meters, it is 

expected that a significant amount of non-linearity may be present in tsunamis large enough to 

cause damage.

Linearity of tsunami waves in deep water has major consequences for the development of 

NCTR’s forecasting methodology. Linearity of tsunami waves implies that basin-wide solutions

computed independently from two or more tsunami sources can be combined into a single 

solution associated with a source constructed as a linear combination of the individual sources. 

By appropriately combining propagation solutions from several unit sources in a way that fits 

observations by tsunami detection instruments deployed in deep water around the generation 

area during an actual event, a more realistic and complex source can be constructed within 

seconds simply by retrieving results of deep-water propagation stored in the database. 

Therefore, by using a pre-computed database of deep-water propagation from unit sources, 

computationally expensive and time consuming calculations of tsunami propagation across the 

ocean basin can be reproduced in a matter of seconds, limiting real-time simulations to the last 

phase of tsunami propagation into shallow coastal waters and inundation onto dry land.

This approach of reconstructing observed tsunami measurements via a linear combination of 

stored unit source simulations allows for the incorporation of tsunami wave elevation values 

reported by tsunami detection instruments as an integral part of the analysis and not merely for 

forecast verification purposes. The down side, however, is the time it takes the tsunami waves 

to travel from the source to the closest detection instrument. However, again, these 

considerations are only relevant to tsunami forecasting in real-time and do not apply to hazard 

assessment since this analysis is based on a previously identified PMT and not on the particular 

scenario generated during an event. 

NCTR’s Propagation Database

NCTR’s current propagation database contains tsunami propagation solutions from over 1,725 

unit sources in the Pacific, Atlantic, and Indian Oceans (Gica, 2008) and allow for the 

reconstruction of tectonic, fore-arc, and aft-arc events from most generation zones.

Each unit source in the propagation database takes into consideration local tectonic 

parameters. In particular, the estimated local dip angle of the subducting plate as it slides 

underneath the continental plate, the orientation of the fault plane with respect to the 

geographical North, or strike angle, and the probable rake angle, and the angle of motion of the 

subducting plate relative to the mounting plate along the fault plane are all estimated based on 

the best available science. Each unit source is designed to have a rupture length of 100 km and 

a width of 50 km. If the subducting and mounting plates are allowed to slip past each other by 
9-4 

1 m, and if one assumes a typical rigidity parameter of the earth crust of 4 × 1011 dynes/cm2

, the 

resulting magnitude of the associated earthquake is MMww = 7.5. This is traditionally considered a 

threshold magnitude for an earthquake that would generate a damaging tsunami away from the 

near field. Once these parameters are defined, the earth crust deformation associated with 

such a seismic event is computed via Okada’s deformation expressions for a semi-elastic half 

space (Okada, 1985). The associated ground deformation typically exhibits an area of uplift 

(usually coinciding with the offshore side of the fault) and an area of subsidence (on the onshore side) that may extend onto dry land, causing permanent inundation of the coastline. 

Figure 9.1 shows the computed ground deformation associated with a single source defined by 

the seismic parameters in Table 9.1. 

Seismic parameters

Longitude (deg) 180.0

Latitude (deg) 25.0

Length (km) 100

Width (km) 50

Dip Angle (deg) 20.0

Rake Angle (deg) 90.0

Strike angle (deg) 180.0

Slip Amount (m) 1

Depth (km) 15

Even though, in principle, an approximate deformation to any large real event can be obtained 

by an adequate combination of unit sources, the flexibility of the reconstruction process is 

Figure 9.1 Typical ground deformation as computed by Okada’s expressions. An area 

of uplift (hot colors) to the North contiguous to an area of subsidence (cold 

colors) to the South can be observed.

Table 9.1 List of seismic parameters corresponding to the sea floor deformation 

shown in Figure 9-1 
9-5 

limited when trying to reconstruct smaller events close in size to the design magnitude of a 

single unit source (MMww = 7.5). In this case, the number of unit sources used in the 

reconstruction must be reduced in order to keep the overall deformation area consistent with the 

magnitude of the event, consequently reducing the number of degrees of freedom in the 

reconstruction problem. The deformation associated with each single unit source is designed to 

produce a credible seismic deformation in itself, such that it will provide an appropriate 

representation of ground deformation for small events. Scaling the deformation associated with 

a single unit source up or down is, perhaps, the only possible degree of freedom.

Tsunami propagation from each unit source is computed by assuming that the ground 

deformation calculated in the manner described above is transferred identically and 

instantaneously to the ocean surface. This approach has been and continues to be standard 

practice in tsunami simulation. In recent years, the advent of deep-ocean tsunami detection 

systems (DARTs) and the deployment of dense networks of real-time Differential Global 

Positioning Systems (GPS), has provided mounting evidence of the relevance of ancillary 

processes associated with large seismic events that may, to some degree, question the 

similarity assumption between the sea floor and sea surface deformation.

For instance, initial conditions derived from direct measurements of recent tsunami events by 

DART systems and used to initialize hydrodynamic models, which in turn have produced 

accurate forecasts of tsunami wave heights on the coast, do not always correlate accurately 

with the estimated ground deformation field derived from seismic data.

Recent ground deformation observations by Japan’s network of differential GPS systems during 

the 2011 Japan event have provided evidence of the complexity of rupture kinematics during a 

large seismic event, with each point on the deformation field experiencing large time-dependent 

displacements of opposite sign before settling onto a final co-seismic deformation.

Processes typically associated with large earthquake events, such as rupture of splay faults and

the occurrence of submarine landslides, seemed to play a fundamental role in the amplification 

and early arrival time of tsunami waves onto the shores of Banda Aceh, during the 2004

Sumatra event, and in smaller events such as the Papua New Guinea event of 1998 (Heinrich 

et al., 2000).

The discrepancy between sea floor and sea surface deformation during an event can be 

compensated by the use of DART systems that directly measure sea surface deformation after 

tsunami generation. This is one of the strengths of NCTR’s forecasting methodology. 

Unfortunately, such an approach is not feasible when conducting hazard assessment studies 

and one has to resort to the traditional method of assuming identical sea floor and surface 

deformations to initialize hydrodynamic models.

 Real-time Simulation Computations

As explained Section 9.2 during a tsunami event, simulation solutions from the generation area 

to the impact coastline are extracted from a pre-computed database. When the tsunami 

approaches the coastline, the linearity parameter, εε, becomes increasingly larger and the linear 

combination of simulations corresponding to individual unit sources is no longer valid. At this 

time, real-time, non-linear calculations of tsunami evolution in shallow water and inundation onto 

dry land are necessary. In NCTR’s forecasting methodology, the development of tsunami 

inundation forecast models is accomplished using high-resolution, highly reliable digital 

elevation data of the area of interest. These forecast models are used to solve a boundary 
9-6

value problem by forcing the models through the boundaries of the grid with the local solution 

extracted from the propagation database. In the case of a site located in the near-field of the 

earthquake, and therefore subject to co-seismic deformation, this deformation is also used as 

an initial condition, effectively resolving an initial and a boundary value problem.

9.4.1 High-resolution Digital Elevation Models

The reason for the use of local high-resolution models is two-fold. First, High-resolution grids, 

with a grid spacing in the order of 30 m, are necessary due to tsunami waves shortening in 

wavelength as they propagate over increasingly shallow waters, as described in Section 0. 

However, finer grid resolution of approximately 10 m is strongly recommended. Grid resolution 

should be high enough so that a minimum number of grid nodes per wavelength is maintained 

throughout the simulation. As tsunami waves become shorter, this constraint requires the use 

of increasingly finer grids. It is also important to keep in mind that high frequency waves are 

likely to be generated in coastal areas as longer waves interact and reflect from local 

bathymetric features with smaller length scales. High-resolution grids will ensure accurate 

resolution of these high frequency waves as well.

Second, some small-scale local bathymetric and topographic features may have a significant

impact on wave propagation and local inundation at the coastal community. Due to their small 

scale in one or both of the horizontal dimensions, these features may easily be absent from 

lower-resolution digital elevation models (DEM). Examples of these features are breakwaters, 

piers, and docks. The presence of these structures will affect wave propagation by sheltering 

parts of the study area from tsunami impact, by redirecting wave energy away from the study 

area or focusing it on specific areas of the site, or even by completely blocking the wave from 

propagating into certain parts of the forecast domain.

9.4.2 Coverage of the Forecast Models

In addition to high-resolution topographic and bathymetric data, area coverage of the forecast 

model grids requires special attention. Two primary factors determine the geographical extent 

of the real-time computational grids, wave height and spatial extent. As wave height increases 

relative to the decreasing depth and current speed becomes faster compared to the slowing 

wave speed,�ggh, tsunami behavior transitions from linear to non-linear. Since tsunami 

behavior captured in NCTR’s database simulations is essentially linear, non-linear local 

calculations at sufficiently deep-water depths should be initiated to guarantee that non-linearity 

has not yet manifested itself in wave propagation. Any manifestation of non-linear behavior 

outside of the coverage area of the non-linear local forecast models will be lost, resulting in a 

deficient forecast prediction.

Another consideration determining the spatial extent of the non-linear model is the shortening

wavelength when the tsunami propagates onto shallow waters. Currently, NCTR’s propagation 

database is computed on a 4-arc-min-resolution grid on the equator with increasing resolution in 

latitudes away from the equator due to the spherical shape of the earth. The present 

configuration guarantees the presence of at least 10 grid nodes along a wave approximately 

75 km long. If tsunami wavelengths are expected to drop below 75 km as water depth 

diminishes, it is recommended that the high-resolution grid of the non-linear model be extended 

into deeper waters where longer wavelengths are found, before propagation database 

simulations degrade due to lack of resolution.
9-7

9.4.3 Development of NCTR’s Forecast Models

In order to be able to run these community-specific forecast models fast enough to comply with 

the stringent time constraints of real-time operational tsunami forecasting, a compromise 

between grid resolution and computational speed must to be found. A sufficient resolution of 

the forecast models will guarantee accurate reproduction of the general tsunami dynamics in the 

study region, but should not be exceedingly fine that it unnecessarily slows down calculations in 

the process of computing minor details that may not be relevant in the overall tsunami dynamics 

and inundation. In order to ensure that these objectives are met by all forecast models, 

development of the models begins with a very high-resolution model of the forecast area. This 

initial high-resolution model is called the Reference Model.

The high resolution and long computational time of the Reference Model does not make it a 

suitable candidate for real-time computations. However, its purpose is to ensure that the correct 

tsunami dynamics are captured and the influence of high frequency waves is correctly reflected 

in the simulation. In order to make sure this is the case, NCTR uses a suite of historical events 

in the Pacific Ocean for which DART and tide gauge data are available. These events are 

simulated using the Reference Model for a specific community and modeling results are then 

verified and validated with tide gauge observations at the site of interest. The set of standard 

historical events contains some of the largest tsunami events of the past century, such as the 

1964 Prince William Sound tsunami and the 1946 Alaska tsunami, but it mainly contains large 

events that occurred in the last 20 years, after the first deployment of DART systems in the 

Pacific Ocean. Table 9.2 shows a list of standard historical events used in forecast model 

development in the Pacific.

The reason for the inclusion of mostly post-DART deployment events is that an accurate and 

credible source definition of the tsunami can be computed using observations from DART

systems in deep water via an inversion process. The inversion solution will yield the optimized 

linear combination of unit sources to reproduce the tsunami as observed by the DART systems. 

Reproduction of pre-DART tsunami sources is more complicated and relies primarily on seismic 

inversions and tide gauge observations of the event. The risks of identifying ground dislocation 

with initial ocean surface displacements are addressed in Section 9.3. Inference of tsunami 

initial sources from tide gauge data as opposed to DART data has the added complexity of 

inversion of non-linear signals since the wave has presumably developed a highly non-linear 

behavior on arrival at the tide-gauge station. These difficulties usually result in more 

cumbersome inversions and less accurate tsunami source solutions.

Once the Reference Model has been developed, the second step is the gradual reduction of 

grid resolution and possibly coverage area to ensure that real-time simulations can be 

computed fast enough so as to provide four hours of tsunami simulation in approximately 10 

minutes of wall-clock time. During the process of reduction of grid resolution, numerical results 

are constantly monitored and compared with the results of the high-resolution Reference Model, 

ensuring that no major features in tsunami dynamics are missing from the coarser models. 

While some information is expected to be lost during this process, it is guaranteed that main 

wave propagation characteristics and reflections from local topographical and bathymetric 

features are present in both models. Once sufficient computational speed has been achieved, 

the resulting grids are tested for stability with a set of 18 large magnitude synthetic events 

(MMww = 9.3) from different locations in the Pacific Ocean, and with a set of five synthetic events in 
9-8

the Atlantic. Additionally, a synthetic MMww = 7.5 event and a micro event MMww = 6.2 are tested to 

ensure resolution of higher frequency waves associated with smaller events.

Table 9.2 List of historical events used in the validation of tsunami forecast models 

in the Pacific Ocean by the NOAA Center for Tsunami Research

Earthquake / Seismic

Event

USGS CMT

Date, Time (UTC) MMww

Epicenter

Date, Time (UTC)

Centroid

1946 Unimak 01 Apr 12:28:56.0

52.75°N, 163.50°W 8.5

1952 Kamchatka 04 Nov 16:58:26.0

52.75°N, 160.05°W 9.0

1957 Andreanov 09 Mar 14:22:31

51.5°N, 175.7°W 8.6

1964 Alaska 28 Mar 03:36

61.05°N, 147.48°W 9.2

1994 East Kuril 04 Oct 13:22:59.54

43.956°N, 147.412°E

04 Oct 13:22:59.54

43.956°N, 147.412°E 8.3

1996 Andreanov 10 Jun 04:03:35

51.564°N, 177.632°W

10 Jun 04:04:03

51.10°N, 177.410°W 7.8

2001 Peru 23 Jun 20:33:14

16.26°S, 73.64°W

23 Jun 20:34:23

17.28°S, 72.71°W 8.4

2003 Hokkaido 26 Sep 19:50:06

41.775°N, 143.904°E

26 Sep 19:50:38

42.21°N, 143.84°E 8.3

2003 Rat Island 17 Nov 06:43:07

51.13°N, 178.74°E

17 Nov 06:43:31

51.14°N, 177.86°E 7.8

2006 Tonga 03 May 15:26:39

20.130°S, 174.164°W

03 May 15:27:03

20.39°S, 173.47°W 8.0

2006 Kuril 15 Nov 11:14:16

46.607°N, 153.230°E

15 Nov 11:15:08

46.71°N, 154.33°E 8.3

2007 Kuril 13 Jan 04:23:20

46.272°N, 154.455°E

13 Jan 04:23:48

46.17°N, 154.80°E 8.1

2007 Solomon 01 Apr 20:39:56

8.481°S, 156.978°E

01 Apr 20:40:39

7.76°S, 156.34°E 8.1

2007 Peru 15 Aug 23:40:57

13.354°S, 76.509°W

15 Aug 23:41:58

13.73°S, 77.04°W 8.0

2007 Chile 14 Nov 15:40:50

22.204°S, 69.869°W

14 Nov 15:41:11

22.64°S, 70.62°W 7.7

2009 Samoa 29 Sep 17:48:10

15.509°S, 172.034°W

29 Sep 17:48:27

15.13°S, 171.97°W 8.1

2010 Chile 27 Feb 06:34:14.55

35.909°S, 72.733°W

27 Feb 06:35:15

35.95°S, 73.15°W 8.8

2011 Tohoku 11 Mar 05:46:24

38.297°N, 142.372°E

11 Mar 05:47:47

38.486°N, 142.597°E 9.0
9-9 

Application to Hazard Assessment

Standard deterministic hazard assessment studies traditionally focus on the study of the impact 

of historical events on the site of interest under the assumption that this set of historical events 

would include the PMT scenario. On occasion, identification of PMT may have been based on 

scientific interpretation of structural parameters of local tectonics (a combined approach is also 

possible). In those cases in which the PMT is shown to exceed PMF from other hazards, this 

PMT should be used as the design scenario. The main deficiency of the first approach is the 

short length of the historical record for this type of event, with return periods of several hundreds 

of years at best for the largest magnitude events. This can easily result in underestimation of 

the PMT in an assessment study. The risks associated with the second approach also come 

from an underestimation of PMT based on geophysical or paleo-seismic records. Examples of 

such an underestimation are the 2004 Sumatra and 2011 Japan events. In the case of the 2004 

Sumatra event, Geist et al. (2007) suggest that the slow rate of convergence of the northern 

section of the Sumatra-Andaman subduction zone associated with very long return periods of 

large tsunamigenic events resulted in under-sampling of the historic catalog. Geist et al. (2007) 

also recognize that, prior to the event of December 2004, it was generally unclear whether 

MMww > 9 earthquakes could occur in highly oblique subduction zones. In light of the previous 

findings, Bird and Kagan (2004) emphasize that all subduction zones, no matter their 

convergence rate, are capable of producing large earthquakes of tsunamigenic magnitude. In 

the case of the 2011 Japan event, the American Nuclear Society found that the PMT used in the 

THA study conducted by the Tokyo Electric Power Company (TEPCO) for nuclear power plants 

in Japan was largely underestimated (American Nuclear Society, 2012). It is likely that subsampling of the historical catalog was at least partially responsible for the misrepresentation of 

PMT in this study, which resulted in tsunami wave heights largely exceeding the design 

parameters for the Fukushima-Daiichi NPP.

In view of the underestimation of PMT by recent hazard assessment studies, NCTR has 

proposed a new methodology that makes use of the NCTR developed forecast tools described 

in previous sections to identify a more conservative PMT than those selected by existing 

methods. The conservatism of the new approach substantially reduces the possibility of 

underestimating the design scenario in the future, effectively limiting the shortcomings of 

previous approaches described above. NCTR’s hazard assessment process proceeds in four

main steps as follows:

1. Following the findings of Geist et al. (2007) and Bird and Kagan (2004), a mega-seismic 

event of magnitude MMww = 9.3 or higher is assumed to be credible from any subduction 

zone in the Pacific or Atlantic Oceans, unless there is sufficient evidence to the contrary. 

These events are reconstructed using a combination of sources in NCTR’s propagation 

database, so that the full transoceanic propagation of any of these events is available in 

seconds from the database. The location of the event is shifted by the length of a unit

source (100 km) along the subduction zone to account for any directionality effects on 

basin-wide tsunami energy distribution. This process yields, in a matter of minutes, 

deep-water propagation results for hundreds of simulations that would otherwise take 

multiple days of computation. 

2. The development of NCTR’s validated tsunami forecast models for specific coastal 

communities provides an excellent tool to compute the nonlinear stage of tsunami 

simulation as the waves approach the shallow waters of the coastal site of interest. 

Since the forecast models by design perform calculations under the time-constraints of a 

real tsunami event, they are capable of computing the last stage of propagation and 
9-10

inundation for the hundreds of potential scenarios extracted from the propagation 

database in a small amount of time. This process can be scripted and automated so 

that runup and inundation from hundreds of potential scenarios are computed in a matter 

of days.

3. Upon completion of step 2, a mapping of the estimated wave height at the impact site 

from each of the scenarios extracted from the propagation database can be generated, 

as represented in Figure 9.1. Since this process should also be automated, it is 

important to examine the results of Figure 9.2 for inconsistencies and outliers in order to 

quality-control the results and ensure that no instabilities have arisen during the process, 

resulting in spurious results. If any inconsistencies in trends or outliers are found in

Figure 9.2, the associated simulations should be carefully quality-controlled by 

examining animations of the results in the propagation and/or inundation grids and 

correcting any spurious values.

4. Based on the results of the previous step, a small set of potential PMTs is selected 

based on their potential for generating large tsunami waves at the site of interest. This 

last small set of events should then be run on the high-resolution Reference Model that 

has been generated as a by-product of the forecast model. Simulations using the 

Reference Model will necessarily be much more computationally expensive than those 

obtained with the forecast models, but the set of simulation scenarios is reduced from 

several hundred to perhaps 10 or 20 cases. Due to the high-resolution of the Reference 

Models, these new simulations will contain the level of detail and high accuracy 

necessary to select a final candidate as PMT for the site of interest.

NCTR developed forecast models for 75 communities along U.S. coastlines. If a NPP is located 

within the boundaries of one of NCTR’s existing inundation models (for instance in the case of 

the Diablo Canyon Power Plant in California), use of the NCTR developed grids could be made 

in a hazard assessment study. If the NPP is not included within the grid boundaries of any of 

the developed forecast models, a new set of grids could be developed following NCTR’s 

forecast model development process, which, in addition to the hazard assessment, would result 

in the creation of a tsunami forecast model for the NPP site. In this case, it is recommended 

that inclusion of the grids in the NCTR-developed SIFT system be requested so that a real-time 

tsunami forecast could be issued to the NPP in the event of a tsunami.
9-11

Additional Recommendations

In addition to the hazard assessment methodology described in this section, the following 

considerations should also be taken into account in a rigorous hazard assessment analysis.

• Due to the vulnerability of NPPs to both maximum tsunami inundation values and 

maximum draw down values for the reasons stated in PMEL-136, it is recommended 

that once a PMT scenario has been identified, two modelling studies be conducted. One 

should use an extremely high water level (see Section 5.6 for details on the 

recommended Antecedent Water Level), while the other should use MLLW or other low 

level consistent with the methodology proposed for high-water level. The 

maximum/minimum values obtained for runup and drawdown when combining results for 

both studies should be used as the expected extreme values.

• An assessment of local sea-level rise as determined by NOAA National Ocean Service, 

the U.S. Global Change Research Program and the Intergovernmental Panel for Climate 

Change (IPCC) during the period of validity of the assessment study should be

conducted and included. If expected sea-level rise has the potential to perceptibly affect

the results of the study, two cases should be modeled. One case should be based on

current sea-level at the time of the investigation to be used in the calculation of extreme 

Figure 9.2 Distribution of tsunami maximum amplitude at a random coastal location in 

the Pacific Ocean from each of the potential scenarios extracted from 

NCTR’s propagation data base. A subset of all the scenarios represented 

can be selected and simulated at high-resolution to arrive at a final PMT for 

the specific site under investigation
9-12

tsunami draw-down values and the other assuming sea-level values expected at the end

of the study’s period of validity to compute an extreme tsunami inundation line.

Corrections to the original DEM water level for extreme tides and long-term sea level 

trend should be included in the DEM prior to the numerical simulations, not added to the 

results after completion of computations at any other level.

• For many locations and events, the largest tsunami wave is not the first one to arrive, but

it could make landfall hours after arrival of the first wave. NCTR’s experience in the

development of tsunami forecast models has shown that reflection of tsunami waves 

from distant coastlines is a key factor driving this behavior. This is particularly true in the

Atlantic and Indian Oceans due to the smaller size of these basins when compared to

the Pacific Ocean. The duration of simulations in any basin should be long enough so 

that the effects of late-arriving reflections to the site of interest are well captured in the

simulation.

• To proceed with appropriate conservatism, it is recommended that bare-earth DEMs be

used in the generation of the modeling grids. Few, if any, tsunami models have the

capability of reproducing the complex hydrodynamic processes that develop as

inundation flows interact with small-scale buildings, harbor facilities, and other landbased structures. For the most part, the effect of these structures is that of reducing the

amount of runup and inundation. It is recommended that these structures be eliminated 

from the DEMs used for hazard assessment since their effect on flooding flows is 

mimicked generally by the use of a friction coefficient in the hydrodynamic models. Use

of the resulting bare-earth grids should result in more conservative values of inundation 

extent (See Section 5.3 for details on the type of structures to preserve in a bare-earth 

DEM).

• In cases in which the initiating event is expected to occur in the near-field of the NPP

location, the NPP site could experience either co-seismic subsidence or uplift associated

with the earthquake. This change in elevation with respect to mean sea level will result

in permanent inundation (or emergence of previously submerged land) that will persist

beyond the tsunami event. In hazard assessment studies, these areas of permanent

inundation should be identified and distinguished from those of temporary inundation 

due to the action of tsunami waves.

• It is also recommended that all DEMs be verified by a field survey of the area of interest 

to ensure that no significant changes in topography have occurred after the most recent

data used in the development of the DEM. Special attention should be given to new 

topographic features absent from the original DEM that can have a significant impact on 

local tsunami dynamics. For instance, the construction of a break-water can significantly 

alter tsunami dynamics in the harbor area by blocking the entrance of wave energy into 

the harbor and by modifying its natural frequencies. Sometimes man-made, land-fill 

areas can be found in ports and harbors that may be absent from the original DEM if 

they are of recent construction. In general, best judgement should be used to determine 

what new features, if any, should be included.

• As expressed earlier in this section and in PMEL-136, it is recommended that, if tsunami 

is found to pose considerable hazard to a NPP site, any new tsunami forecast models

developed in the process of a hazard assessment investigation should be included in the

SIFT tsunami forecast software currently in operation at the Tsunami Warning Centers. 

This would ensure NPP management receives timely and accurate tsunami forecast

during an event.
10-1 

10. ASSESSMENT OF EARTHQUAKE-GENERATED TSUNAMI 

INUNDATION BASED ON PROBABILISTIC OFFSHORE 

TSUNAMI HEIGHT

Background and Objectives

In addition to the deterministic approach, the Probabilistic Tsunami Hazard Analysis (PTHA) is 

an important method to assess the tsunami flooding risks at a NRC facility. Its applications, 

advances, and limitations were addressed as one of the main panel topics at the NRC 

Workshop on Probabilistic Flood Hazard Assessment (PFHA) in 2013, published as

NUREG/CP-0302. At this workshop, Panel 5 "Tsunami Flooding" addressed advanced 

methods for PTHA, with an overall emphasis on defining tsunami hazards at an annual 

exceedance frequency of 10−4 to 10−6. Much of the discussion relates to proper 

characterization of tsunami sources and development of robust tsunami models that can 

accurately simulate tsunami propagation, runup, and inundation. NUREG/CP-0302 also 

addresses the need for implementation and verification of the PTHA. This study describes how 

tsunami forecast tools recently developed at NCTR could be an important component of 

tsunami flooding assessment and of PTHA implementation in particular.

Geist and Parsons (2006) adapted the Probabilistic Seismic Hazard Assessment (PSHA) for 

tsunami hazard assessment. They showed that both the empirical and computational analysis 

could be incorporated into PTHA, independently or jointly. The Seaside study by González et 

al. (2009) has shown the PTHA to be a viable approach to tsunami hazard assessment. 

González et al. (2013) apply a similar approach to assessing the tsunami hazards for Crescent 

City, California. In contrast to their work in 2009, González et al. (2013) considered a few more 

earthquake source scenarios from Japan and implemented a method addressing tidal 

uncertainty. Thio et al. (2010) established a different PTHA approach consisting of a large 

amount (several thousands) of tsunami scenarios that include both epistemic uncertainty

through the use of logic trees as well as aleatory variability. Thio et al.’s (2010) study provides 

probabilistic offshore tsunami heights along California’s coastline. Thio et al. (2010) also 

extended their tsunami offshore heights to estimate the tsunami inundation for a few coastal 

communities in California. PTHA methods are also widely used in Japan (Annaka et al., 2007), 

Australia (Burbridge et al., 2008) and New Zealand (Power et al., 2012). Cadeno (2012) carried 

out detailed inundation modeling and risk assessment for several coastal communities in 

Australia based on Geoscience Australia’s tsunami risk assessment for New South Wales.

As stated in Section 10.2, NCTR’s modeling tools can be used for the long-term tsunami hazard 

assessment, in terms of deterministic (Tang et al., 2009; Uslu et al., 2010) or probabilistic 

(González et al. 2009) approach. These models provide unprecedented quality and scope of 

tsunami hazard assessment for a particular community, or a critical NRC facility, along U.S. 

Pacific and Atlantic coastline. Together with PMEL’s model database of tsunami propagation, 

these models are able to relate the PTHA offshore wave height to onshore flooding zones for 

tsunami hazards associated with certain design return period. PMEL is currently collaborating 

with URS Corporation and American Society of Civil Engineers (ASCE) to explore 

methodologies to develop 2,500-year (2% probability of exceedance in 50 years) tsunami

flooding zones based on max tsunami amplitude at 100 m depth obtained through PTHA and 

their disaggregated tsunami sources. This method can be implemented for probabilistic tsunami 

inundation assessment of NRC facilities.
10-2 

Methodology and Procedure

10.2.1 PTHA Offshore Maximum Tsunami Amplitude

Thio et al. (2010) detail the process of obtaining the PTHA offshore maximum tsunami 

amplitude. The methodology used by Thio et al. (2010) is adopted, for the most part, from the 

PSHA proposed by McGuire (2004), except the PTHA is interested in the exceedance of 

maximum tsunami amplitude. The method of Thio et al. (2010) computes the probabilities in 

terms of the annual exceedance frequency of a Poissonian distribution:

PP = 1 − ee−γγ

where P is the probability of exceedance in a time period t. is the annual rate of exceedance 

that can be obtained from the above equation as:

γγ = (−ll (1 − PP))/tt

This section shows a case study of obtaining the probabilistic tsunami flooding at an annual 

exceedance frequency of 4 x 10-3 (corresponding to an average return period of 2,500 years). 

Of note, this average return period is tied to the performance level of a coastal structure in the 

tsunami flooding zone based on the ASCE standards (ASCE, 2016). Thio and Li (2015) stated 

that the use of the Gutenberg-Richter distribution has its limit to characterize the recurrence of 

large earthquakes due to the short history of the earthquake catalog even when paleo-seismic 

data is included. Alternatively, Thio et al. (2010) applied a rigorous PTHA model using a 

subjective approach of the logic trees to weight the likelihood of the recurrence of large

earthquakes. 

Thio et al. (2010) first identify the subduction zones and setup the subfault partition for the 

earthquake sources, and then employed the shallow water wave models to establish a database 

of Green’s function for each of a set of subfaults that adequately describe the earthquake 

rupture. With a defined earthquake recurrence model, a large set of synthetic scenarios is 

generated to represent the full integration over earthquake magnitudes, locations, and sources, 

for every logic tree branch. The synthetic tsunami waveforms for any slip distribution are then 

summed from individual subfault tsunami waveforms to obtain the maximum tsunami wave 

amplitude along the 100 m water depth offshore. Thio et al. (2010) methods include 

consideration of both aleatory and epistemic uncertainties, which account for uncertainties 

resulting from the random nature of modeling, as well as uncertainties due to incomplete 

understanding of natural processes of the earthquake sources. This approach also provides 

source disaggregation, identifying the source regions and magnitudes that contribute the most 

to those offshore tsunami amplitudes.

As the first step of the probabilistic inundation study, we adopt the PTHA offshore amplitudes 

obtained by Thio et al. (2010). Based on the PTHA source disaggregation, we reconstruct 

tsunami sources to the detail of source parameters so that the reconstructed tsunami scenarios 

provide good approximation of the PTHA offshore amplitudes at a site of interest. As a result, 

we are able to extend the PTHA offshore amplitudes to obtain the 2500-year tsunami disaster 

zone (TDZ ) using tsunami inundation models. 
10-3

10.2.2 Tsunami Inundation Model and ComMIT

MOST can be used to compute the tsunami inundation. Sections 9.3 and 9.4 of this NUREG-CR

describe the development of tsunami propagation database and inundation models using MOST. 

It is worth pointing out that the existing propagation database should not be directly used as 

solutions for water depth shallower than 1,000 m due to the coarse resolution of the grids. 

However, all PTHA offshore amplitudes are available at the 100 m water depth. To solve this 

issue, the existing database was extended to include an additional database of tsunami 

waveforms computed using a grid resolution of 24 arc sec (~ 720 m). For coastlines of interest, 

we develop model grids of 24-arc-sec resolution to compute waveforms at the PTHA offshore 

points, adopting boundary conditions provided by the existing propagation database.

In the present study, the existing unit tsunami scenarios and the extended propagation database 

are used to reconstruct the disaggregated PTHA sources through an inversion method. This 

inversion process searches for a best match between the PTHA offshore wave amplitudes and 

the MOST-computed results, the details of which are provided in the next section.

To develop inundation maps for coastal communities of the Indian Ocean region a community 

model was identified as the primary tool at ICG/IOTWS-II. Subsequently, USAID funded 

PMEL/NOAA to develop such a tool which has been named: ComMIT: COMmunity Model 

Interface for Tsunami (Titov et al., 2011). ComMIT enables government agencies and others in 

the region to run tsunami models, using data from local or remote databases. This approach 

has several advantages. First, it allows communities without a significant cadre of trained 

modelers to build tsunami modeling capability for forecast and hazard assessment. Second, it 

allows communities with restrictions on sharing geo-spatial data to input that data locally and 

not share it with other web-based model users, at the same time but share the model results 

regionally or globally. Finally, and most significantly, the internet-based approach creates a 

virtual regional and global community of modelers using the same tools and approaches to 

understand tsunami threats, all able to share information and insights among themselves.

Tsunami models (both deep-water propagation and inundation) require information on: (1) 

bottom and coastal topography; (2) initial and boundary conditions; and, (3) model run specific 

information such as time-step, spatial resolution and length of model run. The aim of ComMIT 

is to provide an interface which allows for the selection of model input data (initial condition, 

bathymetry grids, etc.) as well as a platform to display model output through a graphical user 

interface (GUI). The interface also allows for internet sharing of the model results and use of 

shared databases. ComMIT has been written in the Java programming language (requiring 

version 1.5) and uses NetCDF format for model input and output thus making ComMIT platform 

independent (i.e. it may be run on different platforms such as MS WINDOWS, MAC OS or 

UNIX). ComMIT may be used with different computational model (or a combination of models) 

with the requirement that the models are able to input and output data in specified format 

(mostly NetCDF format – network Common Data Format), and that input parameters for the 

model are read from a simple text (ASCII) file. Use of universally excepted and standardized 

NetCDF format provides access to a number of open-source software for model data 

presentation and analysis. At present, the MOST model is implemented to work with the 

interface.
10-4

10.2.3 Tsunami Source Inversion and Reconstruction

As previously mentioned, the PTHA approach of obtaining 2,500-year offshore maximum tsunami 

amplitude consists of tens of thousands of numerical results of synthetic scenarios (Thio et al., 

2010). For an inundation study, it is unnecessary and time consuming to compute all synthetic

scenarios used to derive offshore amplitudes. Alternatively, a small subset of these synthetic 

scenarios is used to compute the inundation zone. For example, Thio et al. (2010) applied only 

scenarios that were selected based on their source disaggregation study, and Power et al. (2012) 

chose the largest 100 tsunamis for their probabilistic inundation study. Similarly, only the source 

regions and magnitudes contributing the most to those offshore tsunami amplitudes at sites of 

interest were chosen. Figure 10.1 shows an example of the source disaggregation obtained by 

the PTHA method of Thio et al. (2010). For a site (118.34°W, 34.146°N) in California, this source 

disaggregation map indicates that the most dominating source regions are the Aleutian Trench 

and the Alaska subduction zone. Therefore, for this California site, tsunami sources were 

reconstructed only in these two rupture areas, to produce model results matching the PTHA 

offshore amplitudes.

We use a nonlinear least squares method to realize the reconstruction of the tsunami sources, 

which in turn will be used for inundation computation. Based on the PTHA source disaggregation, 

we first select a group of unit sources in the dominating rupture zones. By the use of the precomputed propagation database (described in Section 2.3), the inversion method then adjusts the 

combination of the slip amount of each unit source until the model results match the PTHA 

offshore amplitudes. The nonlinear least squares method, expressed in the equation below, starts 

Figure 10.1 2500-year PTHA source disaggregation for a site (118.36°W, 34.136°N) in 

California, where the blue bars denote the source contributions (%) to the 

site indicated by the red circle
10-5

with an initial guess of slips for selected unit sources. This provides an initial tsunami source. 

The maximum tsunami amplitudes at every PTHA offshore point can be quickly obtained through 

a linear combination of the pre-computed propagation waveforms weighted by the slip amount. 

These model results are then compared with the PTHA values. The inversion method iteratively 

modifies the slip combination for those selected unit sources until a least squares error is reached 

between the model results and the PTHA offshore amplitudes. This solution of the slip 

combination at the source region is further refined until two conditions are satisfied: (1) the 

absolute error between the model results and the PTHA is less than 20%; and (2) all individual 

model results are greater than 80% of the PTHA values. As a result, the final solution of slip 

combination for the selected unit sources gives a workable tsunami source to compute the 

tsunami inundation.

( )

( ) ( ) j

m

i

j ij i

n

j

j x x

f x t x A

f x f x

−



 



 = ⋅















 =

∑

∑

=

=

1

1

2 2

2

max

min || ( )|| min

η

where ηηiiii(tt) is the wave amplitude time series at point jj due to iith unit source; xi is the slip 

coefficient on the ii

th unit source; and AA is PTHA offshore amplitude at jjth point.

10.2.4 Procedure of Obtaining Probabilistic Tsunami Inundation

Using ComMIT, a computation tool integrating NCTR’s tsunami propagation database and 

inundation models (Titov et al., 2011), the tsunami inundation associated with a design 

probability can be produced based on the offshore tsunami heights of Thio et al. (2010). 

Figure 10.2 shows a flow chart describing the steps to carry out this methodology, which are 

summarized as follows:

1. For a selected site or facility, obtain the tsunami heights associated with a design 

probability at 100-m (or other) water depth offshore from Thio et al. (2010).

2. Relate the tsunami height obtained in (1) to a listing of the disaggregated governing 

attributive seismic scenarios.

3. Use a combination of NCTR’s "unit tsunami sources" to reconstruct each disaggregated 

scenarios in (2) with comparable source location and magnitude.

4. Conduct ComMIT model runs for all sources obtained in (3) and compare the computed 

offshore tsunami height with Thio et al. (2010) with results collected in (1).

5. Adjust by scaling the sources in (4) until ComMIT results agree with Thio et al (2010) at 

all selected offshore points (within an acceptable error range).

6. Rerun ComMIT with the adjusted source in (5) to obtain the inundation limit at the 

selected site or facility.

7. Use the envelope of inundation lines from all model scenarios used in (6) to represent 

the probabilistic tsunami inundation at the selected site or facility.

The propagation database built in ComMIT allows this methodology to be carried out for model 

experiments. The next section gives an example of how this methodology is applied to 

determine the 2,500-year flooding zone at Monterey Bay, California.
10-6

10.2.5 Relevance for NPPs

Similar PTHA methodology and procedures to those described above can be implemented for 

longer return periods of relevance for NPPs. The challenges of PTHA for NPPs mainly arise 

from integration over a much broader range of tsunami sources with varying sizes and 

recurrence rates. The inclusion of both aleatory and epistemic variability in the PTHA deals with 

the natural (physical) uncertainties of the earthquake processes. The aleatory uncertainty

includes consideration of variability in magnitude, slip amount, tidal level and tsunami modeling. 

The epistemic uncertainty includes human’s subjective understanding of the rupture process 

and faulting mechanism. Another challenge worth noting is that, especially along the U.S. East 

Coast of Atlantic, the longer return periods of the THAs for NPPs (an annual exceedance 

frequency of 10-6 to 10-4

) is usually dominated by the landslide-generated tsunamis. As 

discussed in Section 7.2, the challenges of PTHA for landslide tsunamis exist due to lack of 

information to relate the slide volume and horizontal slide area to return period.

Figure 10.2 Flow chart of the methodology to assess tsunami inundation based on 

probabilistic offshore tsunami height
10-7 

Case Study

Monterey Bay, California, is used as an example to demonstrate the feasibility of the 

methodology proposed in this study.

10.3.1 Study Area

The City of Monterey is situated on the Monterey Bay National Marine Sanctuary (Figure 10-3), 

a federally protected ocean area extending 45 km along the coast, where the San Andreas Fault 

System traverses in a northwest-southeast direction and controls much of the overall geologic 

character of the region. This series of sub-parallel faults forms the boundary between the 

Pacific and North American tectonic plates, the former of which is sliding northwest several 

centimeters per year relative to the latter. In the vicinity of the sanctuary, the San Andreas Fault 

System is composed of four fault zones: the San Gregorio Fault, extending predominantly 

offshore from Monterey north to Half Moon Bay; the Monterey-Tularcitos Fault zone, covering a 

wide area from Monterey to Santa Cruz within Monterey Bay; the San Simeon Fault; and the 

San Andreas Fault that is almost entirely onshore in this region. An earthquake probability 

study by the USGS (Working Group on California Earthquake Probabilities, 2003) determined 

that there is a 62% chance of a magnitude 6.7 or greater earthquake occurring on one of the 

faults in the greater San Francisco Bay Area between 2003 and 2032. In this time period, there 

is a 10% chance of a magnitude 6.7 or greater earthquake on the San Gregorio Fault and a 

21% chance of a similar earthquake on the San Andreas Fault.

Monterey is subject to both distant and local tsunami threats. The last major tsunami to hit the 

coast of this area occurred in 1964, which affected the entire California coastline and the 

tsunami waves were particularly high from Crescent City to Monterey (~650 km of coastline, 

Figure 10-4) with heights on the open coast ranging from 2.1 to 6.4 m. The recorded wave 

amplitude at Monterey Bay tide gauge, located on the south side of the Monterey Bay, was

about 1 m, but reached as high as 3.4 m at Santa Cruz Harbor situated on the north side of the 

Monterey Bay. Similarly, the 1946 Unimak tsunami barely produced any noticeable waves at 

Monterey Harbor but reached over 3 m at Santa Cruz. Other recorded tsunami waves in the 

past 20 years are mostly smaller than 0.2 m in amplitude, causing no damage to the coastline.

The submarine canyon offshore of Monterey Bay is identified as a region of mass movement 

features. Slumps, debris flows, and other submarine landslides are concentrated along canyon 

walls and the lower continental slope, with many additional distinct and youthful slumps at the 

base of the headward walls of Monterey Canyon (Greene et al., 2002). Land mass movement 

features in the Monterey Bay region suggest that a potential for tsunami generation exists 

(Greene and Ward, 2003). A small landslide occurred at the head of Monterey Canyon during 

the 1989 Loma Prieta MMww = 6.9 earthquake with a small tsunami ~0.5 m high reported to have 

entered the Moss Landing Harbor and a turbidity current reported to have traveled down the 

canyon axis (Greene and Hicks, 1990; Schwing et al., 1990; Garfield et al., 1994). Ward (2005) 

showed that a 0.1 km3

 of material failure in Monterey Canyon could induce more than 7 m of 

runup along 25 km of the coast or more, posing severe tsunami hazard to the City of Monterey.

The City of Monterey is one of the forecast sites developed at NCTR and included in the 

tsunami forecast software SIFT. In the following, the current PTHA methodology is applied to 

this site and forecasting tools generated during the development of the Monterey forecast model

are used. Wei et al. (2013) provide a detailed explanation of the forecast model development 

for Monterey.
10-8 

10.3.2 Development of the 2500-year Tsunami Inundation Zone from the Probabilistic 

Offshore Tsunami Height

10.3.2.1 2500-year Offshore Tsunami Height at Monterey Bay

Thio et al. (2010) provide 2,500-year tsunami heights for 16 offshore locations (Figure 10-4. 

shows that these points are located at water depths ranging from 21 to 69 m, and the 

exceedance tsunami heights are between 4.41 and 6.07 m. Thio et al. (2010) uses a grid 

resolution of 120 arc sec (~3700 m) sampled from ETOPO2 (NGDC) to compute the tsunami

heights offshore. All these points are covered by the intermediate grid of NCTR’s Monterey 

model, which employs a grid resolution of 18 arcsec (40 m) sampled from NGDC’s tsunami 

DEM. Different data sources and grid resolutions usually lead to discrepancies in water depth. 

Table 10.1 compares the depth differences at all 16 locations between Thio et al. (2010) and the 

NCTR model. The large depth differences occur at locations where the ocean bottom is 

channelized due to the deep canyons offshore, such as location P1 (Carmel Bay), P5, and P7 

(east end of Monterey Canyon). These differences pose challenges for model comparison 

between NCTR models and Thio et al. (2010). Therefore, it is important to examine the model 

results in terms of error range rather than a point-by-point match.

Figure 10.3 Aerial photo overlooking Monterey Harbor
10-9 

Table 10.1 The 2,500-year tsunami heights at 16 offshore locations obtained by Thio 

et al. (2010). The water depths at these locations are compared with those 

extracted from NCTR’s tsunami model. 
10-10

10.3.2.2 Source Disaggregation Using ComMIT

An advantage of Thio et al.’s (2010) probabilistic approach is that it allows the source 

disaggregation, which identifies the most contributive source locations and magnitudes to the 

probabilistic offshore tsunami heights. For Monterey Bay, Thio et al. (2010) show that the most 

contributive sources for Monterey 2500-year tsunami heights are sources in Alaska-Aleutian 

and Kuril-Kamchatka source regions, as seen in Figure 10.5. In this study, two sources, one in 

Alaska-Aleutian and the other in Kuril-Kamchatka, are used to showcase the methodology.

For the Alaska-Aleutian scenario, 30 unit sources, a rupture area of 1,500 km by 100 km, were 

selected to construct sources of different magnitudes and varying slip amounts (uniformly 

applied to each unit source) (Figure 10.6). Similarly, as suggested by Thio et al. (2010), 14 unit 

sources, a rupture area of 700 km by 100 km, were used to construct sources in the KurilKamchatka source region (Figure 10.7). Both source selections can be quickly done using 

ComMIT.

Figure 10.4 Locations at Monerey Bay, California, where the offshore tsunami heights 

were obtained (from Thio et al., 2010). The circles are areas of large water 

depth discrepancies between Thio et al. (2010) and NCTR’s tsunami model.
10-11

Figure 10.5 Source disaggregation of 2,500-year offshore tsunami heights at Monterey, 

California courtesy of Thio et al. (2010). The vertical axis indicates how 

much each of the sources contributes in percentage, probabilistically, to 

the tsunami impact at Monterey.
10-12

Figure 10.6 Source selection in Alaska-Aleutian using ComMIT

Figure 10.7 Source selection in Kuril-Kamchatka using ComMIT
10-13

10.3.2.3 Source Inversion to Match Probabilistic Offshore Tsunami Heights

In the process of source inversion described in Section10.2.3, the slip amount was inverted at 

the source region to reach an agreement of the 2500-year offshore tsunami heights between 

ComMIT and the PTHA offshore amplitudes in Thio et al. (2010). Figure 10.8 shows the 

comparison of the tsunami heights for a M9.5 Alaska-Aleutian scenario with a slip amount tuned 

to 33.25 m. The computational results obtained from ComMIT indicate that this scenario 

produces very comparable results with those of Thio et al. (2010). Most of the errors are within 

10%, and the large errors (20-30%) coincide with large water depth discrepancies (Table 10.2). 

For all 16 locations, the average error is about 12.1% and the root-mean-square error is about 

0.89 m. This bolsters the feasibility of the methodology proposed in this study, but also 

highlights the need of matching bathymetric and topographic grids.

The methods of Thio et al. (2010) also provide a map of probabilistic offshore wave periods. By 

matching the wave periods, the tsunami scenario obtained through the above procedures can 

further be constrained and tuned for the inundation mapping effort. Although the off-shore wave 

periods of Thio et al. (2010) were not available for this report, the spectrum of the tsunami wave 

computed at all 16 locations for the M9.5 Alaska-Aleutian scenario was analyzed. f shows that 

the dominant wave periods have a range of 30 to 60 minutes. P1 has the longest wave period 

of 61 minutes. The waves in the north of Monterey Canyon (P5-12) have a common wave 

period of 51 minutes, while the waves near Monterey Harbor (P2-4) and northwest of Monterey 

Bay (P13-16) have a wave period of half an hour. It is worth noting that the longer wave period 

in the north of Monterey Canyon is likely attributed to a well-known bay resonance triggered by 

tsunami waves (Tolkova and Power, 2011).

Figure 10.8 Comparison of the 2,500-year offshore tsunami heights for the M9.5

 Alaska- Aleutian scenario
10-14

Site Thio et al. (2010) ComMIT Difference

(%)

P1 5.68 4.51 -20.7

P2 4.64 5.01 9.6

P3 4.64 4.70 1.3

P4 4.77 4.55 4.6

P5 5.92 4.15 -30.0

P6 4.60 4.75 3.2

P7 4.60 5.66 23.0

P8 4.61 4.91 6.6

P9 7.39 6.25 -15.5

P10 4.79 4.67 -2.5

P11 7.39 7.24 -2.1

P12 7.23 8.13 12.4

P13 4.41 4.36 -1.1

P14 5.43 3.85 -29.2

P15 4.41 4.78 8.5

P16 6.04 4.66 -22.9

Table 10.2 Comparison and errors of the 2500-year offshore tsunami heights between 

ComMIT and Thio et al. (2010) for the M9.5 Alaska-Aleutian sources
10-15

Figure 10-10 further investigates the same procedure for a Kuril-Kamchatka scenario (). For all 

16 locations, the average error is about 14.4% and the root-mean-square error is about 1.1 m 

(Table 10.3). To reach similar accuracy for the probabilistic tsunami heights, the slip amount 

needs to be tuned up to 142 m with a resulting magnitude of 9.7. Such a scenario is unrealistic 

for the Kuril-Kamchatka source region.

Figure 10.9 Spectrum analysis of the tsunami waves at all 16 locations where Thio et al. 

(2010) results are available
10-16

Figure 10.10 Comparison of the 2500-year offshore tsunami heights for the M9.7 KurilKamchatka scenario
10-17

Site Thio et al. (2010) ComMIT

Difference

(%)

P1 5.68 3.84 -32.4

P2 4.64 4.68 0.9

P3 4.64 4.30 -7.4

P4 4.77 4.32 -9.4

P5 5.92 3.96 -33.1

P6 4.60 4.49 -2.5

P7 4.60 5.26 14.3

P8 4.61 4.53 -1.7

P9 7.39 6.76 -8.6

P10 4.79 7.06 47.2

P11 7.39 7.01 -5.1

P12 7.23 7.56 4.5

P13 4.41 4.23 -4.1

P14 5.43 3.84 -29.3

P15 4.41 4.55 3.0

P16 6.04 4.42 -26.9

10.3.2.4 Tsunami Inundation at Monterey Bay

From the aforementioned steps, the M9.5 Alaska-Aleutian scenario was identified as a most 

probable source for the 2500-year tsunami heights offshore of Monterey Bay, California. Figure 

10-11 shows the maximum tsunami wave amplitude in the Pacific due to this scenario, which 

directs the main energy towards the U.S. west coast with far reach to the Pacific coast of South 

America.

As explained in Chapter 0, assembling unit sources in the propagation database provides the 

initial and boundary conditions for NCTR’s tsunami models. In this study, a 10-m reference 

model was used to compute the tsunami inundation at Monterey, California. Figure 10-12

shows the computed maximum water level, up to 7 m above mean high water, for the AlaskaAleutian scenario Figure 10-13 indicates the flow speed on land due to the same scenario could 

reach up to 7-8 m/s with complex flow pattern in the area of Monterey Harbor.

Table 10.3 Comparison and errors of the 2500-year offshore tsunami heights between 

ComMIT and Thio et al. (2010) for the M9.7 Kuril-Kamvhatka source
10-18

Figure 10-14 shows the ultimate product of the method proposed in this study: the inundation 

limit associated with the 2500-year design return period. Also shown in this figure are the runup 

heights along the inundation limit.

Figure 10.11 Computed maximum tsunami wave amplitude in the Pacific due to the 

M9.5 Alaska-Aleutian source
10-19

Figure 10.12 Computed maximum tsunami water level along Monterey Bay’s coastline 

due to the M9.5 Alaska-Aleutian source
10-20

Figure 10.13 Computed maximum tsunami flow speed along Monterey Bay’s coastline 

due to the M9.5 Alaska-Aleutian source
10-21

Figure 10.14 The 2500-year tsunami inundation zone for Monterey Bay, California

11-1

11. CONCLUSIONS

An introduction to the current state of the science with regards to tsunami hazard has been 

presented in this report with a special emphasis on its application to the mitigation of the hazard 

to NPP. Any hazard assessment study for an NPP site should take into account the 

peculiarities of this type of facilities and incorporate them into the analysis. The main 

peculiarities that need to be taken into account are:

• Operations at NPPs can be seriously disrupted not only by maximum values of tsunami 

wave elevation and currents, but also by the minimum values reached by both of these 

quantities during a tsunami event. Minimum values are particularly relevant to the 

cooling system of an NPP, which can be easily affected by a cooling flow deficit.

• If any significant sea-level change is expected in the area of the NPP site during the 

period of validity of the study, this correction should be incorporated into the analysis. 

Long term changes to mean sea level at a particular location could be due to different 

phenomena with periods ranging from several months to several years, such as long 

term tidal fluctuations. Other phenomena affecting long term sea level trends are 

changes in offshore currents and winds, change in ocean circulation patterns and 

climate change.

The rest of the recommendations presented in this report apply equally to tsunami hazard 

assessment of other type of facilities, however, the assessment of the tsunami impact on an 

NPP should be particularly observant and scrupulous due to the added potential for catastrophic 

consequences.

Whether the PMT generating source is of seismic or landslide origin, an approach in which the 

very worst credible source and most conservative set of model parameters are used in the 

estimation of the impact is recommended. Results obtained via this approach can then be 

scaled down or re-computed if sufficient scientific evidence can be found to justify the reduction 

of the initial PMT case.

In the same spirit, for landslide generated PMTs, initial hazard assessments based on onedimensional modeling approaches, in which the entire shelf is assumed to fail simultaneously, 

and no lateral spreading of the wave energy is allowed should be acceptable as an initial PMT. 

Only once sufficient information is available, so that one can accurately limit the lateral extent of 

the slide, should a two-dimensional modeling approach be adopted.

In the identification of a PMT for a specific location, the tsunami forecasting tools developed at 

the NCTR can be of great assistance by allowing the modeler to examine the effect of hundreds 

of tsunami scenarios on any location with a minimal development and computational cost.

The use of non-dispersive NSW models for tsunami simulation is recommended for seismically 

generated tsunamis. These tsunamis are not expected to exhibit a strong dispersive behavior in 

deep water and NSW models are expected to err on the side of caution i.e., overestimate 

extrema if any significant dispersive behavior is manifested in shallow waters.

When it is determined that the PMT is likely to be generated by a landslide, the choice of a 

physically dispersive model such as those based on Boussinesq approximations of the 

equations of fluid motion, is recommended. The potential for larger amplitude and shorter 

frequency waves that in the case of seismic tsunamis will result in the generation of strong 

dispersive effects that should be accounted for in the model.
11-2

If a probabilistic rather than deterministic assessment is to be conducted, we propose a 

methodology using NCTR’s tsunami forecast models and tools to develop the tsunami 

inundation zone based on PTHA offshore tsunami heights for the case of seismically-generated 

tsunamis. This approach utilizes ComMIT, a tool that integrates NCTR’s tsunami propagation 

database and inundation models, to tune the disaggregated earthquake sources according to 

the PTHA offshore tsunami heights. The envelope of computed tsunami inundation from 

selected sources then defines the probabilistic tsunami inundation zone. The study using 

Monterey Bay, California strongly encourages the feasibility of the proposed methodology. The

future implementation of this method to assess the PTHA tsunami inundation at NRC-regulated

facilities is recommended.
12-1

12. REFERENCES

American Nuclear Society, (2012). Fukushima Daiichi: ANS Committee Report., 

http:/fukushima.ans.org/report/Fukushima_report.pdf

Annaka, T., Satake, K., Sakakiyama, T., Yanagisawa, K., and Shuto, N., (2007). Logic-tree 

approach for probabilistic tsunami hazard analysis and its applications to the Japanese Coasts.

Pure and Applied Geophysics, 164, 577-592.

Arakawa, A., and V. R. Lamb, (1977). Computational design of the basic dynamical processes 

of the UCLA general circulation model. Methods in Computational Physics, J. Chang, Ed., 

Academic Press, 174-267.

Atwater, B.F., Satoko, M.-R., Satake, K., Yoshinobu, T., Kazue, U. and Yamaguchi, D.K., 

(2005). The Orphan Tsunami of 1700-Japanese Clues to a Parent Earthquake in North 

America. U.S. Geological Survey Professional Paper 1707, 144 p.

Berger, M. J., and R. J. LeVeque (1998). Adaptive mesh refinement using wave-propagation 

algorithms for hyperbolic systems. SIAM J. Num. Anal., 35:2298-2316.

Bird, P., and Y. Y. Kagan, (2004). Plate-tectonic analysis of shallow seismicity: apparent 

boundary width, beta-value, corner magnitude, coupled lithosphere thickness, and coupling in 7 

tectonic settings, Bull. Seism. Soc. Am., 94, 2380- 2399.

Birknes, J., G. Pedersen, (2006). A particle finite element method applied to long wave runup,

International Journal for Numerical Methods in Fluids, Vol. 52(3), p 237-261.

Borrero, J.C., Bourgeois, J., Harkins, G., and Synolakis, C.E. (1997). How Small-Scale 

Bathymetry Affected Coastal Inundation in the 1992 Nicaraguan Tsunami, AGU Fall Meeting,

San Francisco.

Briggs, M. J., Synolakis, C. E., Harkins, G. S. and Green, D. (1995). Laboratory experiments of 

tsunami runup on a circular island, Pure and Applied Geophysics, 144 (3/4), 569-593.

Burbridge, D., Cummings, P.R., Mleczko, R. and Thio, H.K., (2008). A probabilistic tsunami 

hazard assessment for Western Australia, Pure and Applied Geophysics, 165, 2059-2088.

Cadeno, (2012). NSW tsunami inundation modeling and risk assessment, NSW Office of 

Environmental and Heritage and NSW State Emergency Service, LJ2874/Rep2703, p56. 

Carrier, G. F., (1966). Gravity waves on water of variable depth. J. Fluid Mech, 641-659.

Chen, Q., Kirby, J. T., Dalrymple, R. A., Shi, F. and Thornton, E. B., (2003). Boussinesq 

modeling of longshore currents, J. Geophys. Res, 108 (C11), 3362, 

doi: 10.1029/2002JC001308

Chen, Q., Kirby, J. T., Dalrymple, R. A., Kennedy, A. B., & Chawla, A. (2000). Boussinesq 

modeling of wave transformation, breaking, and runup. II: 2D. Journal of Waterway, Port, 

Coastal, and Ocean Engineering, 126(1), 48-56.
12-2

Chen, Q., Dalrymple, R. A., Kirby, J. T., Kennedy, A. and Haller, M. C., (1999). Boussinesq 

modeling of a rip current system, J. Geophys. Res, 104, 20,617 -20, 637.

Chen, Q., Madsen, P. A., Schaffer, H. A., and Basco, D. R. (1998). Wave-current interaction 

based on an enhanced Boussinesq approach, Coast. Eng., 33, 11-39.

Chen, Y. and Liu, P. L.-F. (1995). The unified Kadomtsev-Petviashvili equation for interfacial 

waves, J. Fluid Mech., 288, 383-408, doi: 10.1017/S0022112095001182

Cheung, K., Phadke, A., Wei, Y., Rojas, R., Douyere, Y., Martino, C., Houston, S., Liu, P. L.-F., 

Lynett, P., Dodd, N., Liao, S., and Nakazaki, E. (2003). Modeling of Storm-induced Coastal 

Flooding for Emergency Management, Ocean Engineering, 30, p. 1353-1386.

Cienfuegos R., Barthelemy, E., Bonneton P. (2007). A fourth-order compact finite volume 

scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: Boundary 

conditions and model validation, Int. J. Num. Meth. Fluids, 53 (9), 1423-1455.

Coastal Engineering Manual (2002). Engineer Manual 1110-2-1100, U.S. Army Corps of

Engineers, Washington, DC. (6 volumes), http://www.usace.army.mil/inet/usace-docs/engmanuals/cecw.htm

Cross, R. H., (1967). Tsunami surge forces. J. Waterways and Harbor Division, ASCE, 93 (4), 

201-231

Cruden, D.M, and Varnes, D.J. (1996). Landslide types and processes, Special Report 247: 

Landslides: Investigation and Mitigation, Transportation Research Board, Washington, D.C.

Dally, W.R., Dean, R.G., and Dalrymple, R.A. (1985). Wave height variation across beach of 

arbitrary profile, J. Geophys. Res., 90 (C6), 11917-11927.

Dalrymple, R.A. and B.D. Rogers (2006). Numerical Modeling of Water Waves with the SPH 

Method. J. Coast. Eng., 53/2-3, 141-147, DOI: 10.1016.

Danielsen F, Sorensen MK, Olwig MF, Selvam V, et al. (2005) The Asian tsunami: a protective 

role for coastal vegetation, Science, 310, 643 pp.

Dawson, C., J.J. Westerink, J.C. Feyen, and D. Pothina (2006). Continuous, Discontinuous and 

Coupled Discontinuous-Continuous Galerkin Finite Element Methods for the Shallow Water 

Equations, Int. J. Numer. Meth. Fluids, 52, 63-88.

Dean, R. G. and R.A. Dalrymple (1991). Water Wave Mechanics for Engineers and Scientists, 

Englewood Cliffs: Prentice-Hall, Inc., Reprinted Singapore: World Scientific Publishing Co.,

ISBN981-02-0420-5.

Elgar, S., and R.T. Guza (1985). Shoaling gravity waves: a comparison between data, linear 

finite depth theory and a nonlinear model, J. Fluid Mech., 158, 47-70.

Fernando, H.J.S., J.L. McCulley, S.G. Mendis, and K. Perera (2005). Coral poaching worsens 

tsunami destruction in Sri Lanka. Eos Trans. AGU, Trans. AGU, 86(33) 30, doi:

10.1029/2005eo330002.
12-3

Fujima, K., Kenji Masamura and Chiaki Goto (2002). Development of the 2d/3d Hybrid Model 

for Tsunami Numerical Simulation, Coast. Eng., 44 (4) 373-397.

Garfield, N., T. A. Rago, K. J. Schnebele, and C. A. Collins (1994). Evidence of a turbidity 

current in Monterey Submarine Canyon associated with the 1989 Loma Prieta earthquake, 

Cont. Shelf Res, 14, 673-686.

Geist, E., Lynett, P., and Chaytor, J. (2009). Hydrodynamic Modeling of Tsunamis from the 

Currituck Landslide, Marine Geology, 264 (1,2), 41-52.

Gica, E., M. Spillane, V.V. Titov, C. Chamberlin, and J.C. Newman (2008). Development of the 

forecast propagation database for NOAA’s short-term Inundation Forecast for Tsunamis (SIFT). 

NOAA Tech. Memo. OAR PMEL-139, 89 pp.

Geist, E.L., V.V. Titov, D. Arcas, F.F. Pollitz, and S.L. Bilek, (2007). Implications of the 26 

December 2004 Sumatra Andaman Earthquake on tsunami forecast and assessment models 

for great subduction-zone earthquakes, Bull. Seism. Soc. Am., 97(1A), doi:

10.1785/0120050619, S249-S270.

Geist, E.L. and T. Parsons (2006). Probabilistic analysis of tsunami hazards, Natural Hazards, 

7, 272-314.

González, F.I., R.J. LeVeque, and L.M. Adams (2013) Probabilistic Tsunami Hazard 

Assessment (PTHA) for Crescent City, CA (Final Report on Phase I). Pilot Study funded by 

BakerAECOM, p70.

González, F.I., E.L. Geist, B. Jaffe, U. Kanoglu, H. Mofjeld, C.E. Synolakis, V.V. Titov, D. Arcas, 

D. Bellomo, D. Carlton, T. Horning, J. Johnson, J. Newman, T. Parsons, R. Peters, C. 

Peterson, G. Priest, A. Venturato, J. Weber, F. Wong, and A. Yalciner, (2009). Probabilistic 

tsunami hazard assessment at Seaside, Oregon, for near- and far-field seismic sources, J. 

Geophys. Res., 114 (C11023), doi: 10.1029/2008JC005132.

González, F. I., Satake, K.; Boss, E. F.; Mofjeld, H. O., (1995). Edge wave and non-trapped 

modes of the 25 April 1992 Cape Mendocino tsunami, PAGEOPH, 144 (3-4)409-426.

Gopalakrishnan, T.C. and C.C. Tung (1983). Numerical analysis of a moving boundary problem 

in coastal hydrodynamics. Int. J. Num. Meth. Fluids., 3, 179-200.

Goring, D.G. and Raichlen, F. (1992). Propagation of long waves onto shelf, J. Waterway, 

Port, Coast, and Ocean Eng., 118 (1), 43-61.

Goring D.G. (1978). The propagation of long waves onto a shelf, Ph.D. thesis, California 

Institute of Technology, Pasadena, California.

Greene, H.G. and S.N. Ward (2003). Mass movement features along the central California 

margin and their modeled consequences for tsunami generation. Submarine Mass 

Movements and Their Consequences, Kluwer Academic Publishers, Netherlands.

Greene, H.G., N.M. Maher and C.K. Paull (2002). Physiography of the Monterey Bay National 

Marine Sanctuary and implications about continental margin development. Marine Geology,

181, 55-82.
12-4

Greene, H.G. and K.R. Hicks (1990). Ascension-Monterey canyon system: History and 

development, in Geology and Tectonics of the Central California Coastal Region, San Francisco 

to Monterey, Am. Assoc. Petro. Geol. (AAPG), GB 67: 229-250

Grilli, S. T., Ioualalen, M., Asavanant, J., Shi, F., Kirby, J. T. and Watts, P. (2007). Source

constraints and model simulation of the December 26, 2004 Indian Ocean tsunami, J. 

Waterway, Port, Coast., Ocean Eng.,ASCE, 414-428, DOI: 10.1061/(ASCE)0733-

950X(2007)133:6(414).

Guignard, S., Grilli, S., Marcer, R., and Rey, V. (1999). Computation of shoaling and breaking 

waves in nearshore areas by the coupling of BEM and VOF methods, Proc. 9th Offshore and 

Polar Eng. Conf., 3, 304-309.

Hammack, J. and Segur, H. (1978). Modelling criteria for long water waves, J. Fluid Mech.,

84(2), 359-373.

Heinrich, P., A. Piatanessi, E. Okal, and H. Hebert (2000). Near-field modeling of the 17 July 

1998 tsunami in Papua New Guinea, Geophys. Res. Lett., 27(19), 3037-3040, 

doi:10.1029/2000GL011497.

Hibberd, S. and Peregrine, D. H. (1979). Surf and run-up on a beach, J. Fluid Mech., 95, 323-

345.

Horrillo, J., Kowalik, Z, and Shigihara, Y. (2006). Wave Dispersion Study, in the Indian OceanTsunami of December 26, 2004, Marine Geodesy, 29, 149-166, DOI: 

10.1080/01490410600939140

Hsu, T-J, Sakakiyama, T., and Liu, P. L.-F. (2002). A numerical model for wave motions and 

turbulence flows in front of a composite breakwater, Coast. Eng. J., 46(1), 25-50, doi: 

10.1016/S0378-3839(02)00045-5.

Imamura, F. (1995). Review of tsunami simulation with a finite difference method, long-wave 

runup models, World Scientific, 25-42.

Imamura, F., Shuto, N. and Goto, C. (1988). Numerical simulations of the transoceanic

propagation of tsunamis, Proc. 6th Congress Asian and Pacific Regional Division, IAHR, Japan, 

265-272.

Ioualalen, M., Asavanant, J., Kaewbanjak, N., Grilli, S.T., Kirby, J.T. and P. Watt (2007). 

Modeling the 26 December 2004 Indian Ocean tsunami: Case study of impact in Thailand, J. 

Geophys. Res., 112 (C07024), 1-21, doi: 10.1029/2006JC003850.

Japan Society of Civil Engineers (2002). Tsunami Assessment Method for Nuclear Power

Plants in Japan, The Tsunami Evaluation Subcommittee, The Nuclear Civil Engineering 

Committee, JSCE (Japan Society .of Civil Engineers), 

http://committees.jsce.or.jp/ceofnp/system/files/JSCE_Tsunami_060519.pdf

Johnson, R. S., (1997). A Modern Introduction to the Mathematical Theory of Water Waves, 

Cambridge Texts in Applied Mathematics, ISBN:9780521598323
12-5

Johnson, D.B., Raad, P.E., and Chen, S. (1994). Simulation of impacts of fluid free surfaces 

with solid boundaries, Int. J. Num. Meth. Fluids., 19, 153-176.

Johnson, R. S. (1972). Some numerical solutions of a variable-coefficient Korteweg-de Vries 

equation (with application to solitary wave development on a shelf), J. Fluid Mech., 54 (1), 81-

91.

Kanamori, H. (1977). The energy release in great earthquakes, J. Geophys. Res., 82, 2981-

2987, doi:10.1029/JB082i020p02981.

Kennedy, A. B., Chen, Q., Kirby, J.T. and Dalrymple, R. A., (2000). Boussinesq modeling of 

wave transformation, breaking, and runup. I: 1D., J. Waterway, Port, Coastal, and Ocean Eng., 

ASCE, 126, 39-47.

Kim, D.-H. and P. Lynett (2012). A sigma-Coordinate Transport Model Coupled with Rotational 

Boussinesq-type Equations. Environmental Fluid Mechanics, doi: 10.1007/s10652-012-9256-1.

Kim, D.-H., P. J. Lynett, S. A. Socolofsky (2009). A depth-integrated model for weakly 

dispersive, turbulent, and rotational fluid flows, Ocean Modelling, 27 (3-4), 198-214.

Kirby, J. T. and Dalrymple, R. A. (1994). Combined Refraction/Diffraction Model REF/DIF 1, 

Version 2.5, Documentation and User’s Manual, Research Report No. CACR-94-22, Center for 

Applied Coastal Research, Dept. Civil Eng., University of Delaware, Newark.

Kirby, J. T. (1986). Higher-order approximations in the parabolic equation method for water 

waves, J. Geophys. Res., 91, 933-952.

Komen, G., L. Cavaleri, M. Donelan, K. Hasselmannn, S. Hasselmannn and P.A. E. M. 

Janssen, (1994). Dynamics and Modeling of Ocean Waves, Cambridge University Press,199p.

Korycansky, D. G. and Lynett, P., (2007). Runup from Impact Tsunami, Geophys. J. Int., 

170 (3): 1076-1088, doi: 10.1111/j.1365-246X.2007.03531.x. 

Kowalik, Z., and Murty, T.S. (1993). Numerical simulation of two-dimensional tsunami runup, 

Marine Geodesy, 16, 87-100.

Kulikov, E. (2006). Dispersion of the Sumatra tsunami waves in the Indian Ocean detected by 

satellite altimetry, Russ. J. Earth Sci., 8, ES4004, 1-5. doi: 10.2205/2006ES000214

Lay, T., Kanamori, H., Ammon, C., Nettles, M., Ward, S., Aster, R., Beck, S., Bilek, S., 

Brudzinski, M., Butler, R., DeShon, H., Ekstro ́m, G., Satake, K., Sipkin, S., (2005). The 

Great Sumatra-Andaman Earthquake of 26 December 2004, Science, 1127-1133, 

doi:10.1126/science.1112250.

LeVeque, R.J. and D. L. George, (2004). High-resolution finite volume methods for the shallow 

water equations with bathymetry and dry states, Advanced Numerical Models for Simulating 

Tsunami Waves and Runup, v . 10, Advances in Coastal and Ocean Engineering, 43-73.

LeVeque, R. J. (2002). Finite Volume Methods for Hyperbolic Problems, Cambridge University 

Press.
12-6

Lin, P., and Liu, P.L.-F., (1998a). A numerical study of breaking waves in the surf zone, J. Fluid

Mech., 359, 239-264.

Lin, P., and Liu, P. L.-F., (1998b). Turbulent transport, vorticity dynamics, and solute mixing 

under plunging breaking waves in surf zone, J. Geophys. Res., 103 (C8), 15677-15694.

Liu, P. L.-F., P. Lynett, H. Fernando, B. E. Jaffe, H. Fritz, B. Higman, R. Morton, J. Goff, and C. 

Synolakis (2005). Observations by the International Tsunami Survey Team in Sri Lanka, 

Science, 308, 1595.

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

rundown generated by three-dimensional sliding masses, J. Fluid Mech., 536: 107-144.

Liu, P.L.-F. and Cheng, Y. (2001). A numerical study of the evolution of a solitary wave over a 

shelf, Phys. Fluids, 13 (6), pp. 1660.

Liu, P.L.-F. (1995). Model equations for wave propagation from deep to shallow water, Adv. 

Coast. Ocean Eng., 1, 125-157.

Liu, P. L.-F., Cho, Y.-S., Briggs, M. J., Kanoglu, U. and Synolakis, C. E. (1995). Runup of

solitary waves on a circular island. In J. Fluid Mech., 320, 259-285.

Liu, P. L.-F., Cho, Y.-S., Yoon, S. B. and Seo, S. N., (1994). Numerical simulations of the 1960 

Chilean tsunami propagation and inundation at Hilo, Hawaii, Progress in Prediction. Disaster 

Prevention and Warning, Kluwer Academic Publishers, 99-115.

Liu, P. L.-F. (1994). “Model equations for wave propagation from deep to shallow water.” 

Advances in Coastal Engineering, P. L.-F. Liu, ed., Vol. 1, World Scientific, Singapore, 125–157.

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

long-wave run-up, J. Fluid Mech., 229, 675-688.

Liu, P.L.-F., Yoon, S.B. and Kirby, J.T. (1985). Nonlinear refraction-diffraction of waves in 

shallow water. J. Fluid Mech., 153, 184-201.

Lynett, P., (2007). The Effect of a Shallow Water Obstruction on Long Wave Runup and 

Overland Flow Velocity J. Waterway, Port, Coast., Ocean Eng., (ASCE), 133(6), 455-462.

Lynett, P. (2006). Modeling Using High-Order Boussinesq Equations, J. Waterway, Port, 

Coast., Ocean Eng., (ASCE), 132(5), 348-357.

Lynett, P., Borrero, J., Liu, P. L.-F., and Synolakis, C.E. (2003). Field Survey and Numerical 

Simulations: A Review of the 1998 Papua New Guinea Tsunami, PAGEOPH, 160, 2119-2146.

Lynett, P. and Liu, P. L.-F. (2002). A Numerical Study of Submarine Landslide Generated 

Waves and Runup, Proc. Royal Society of London, 458, 2885-2910.

Losada, M.A., V. Vidal, and R. Medina (1989). Modeling Experimental study of the evolution of 

a solitary wave at an abrupt junction, Coastal Eng., 46(2), 89-107.
12-7

Madsen, P.A., Sørensen, O.R., Schaffer, H.A. (1997). Surf zone dynamics simulated by a 

Boussinesq-type model: Part I. Model description and cross-shore motion of regular waves, 

Coastal Eng., 32, 255-287.

Madsen, P. A., Murray, R., and Sørensen, O. R. (1991). A new form of the Boussinesq 

equations with improved linear dispersion characteristics (Part 1), Coastal Eng., 15, 371-388.

Madsen, O.S., and C. C. Mei (1969). The transformation of a solitary wave over an uneven 

bottom, J. Fluid Mech., 39, 781-791.

Mansinha, L. and Smylie, D. E. (1971). The displacement fields of inclined faults, Bull. Seism.

Soc. Am., 61, 1433-1440.

Matsuyama, M., Ikeno, M., Sakakiyama, T., Takeda, T. (2007). A study on tsunami wave fission 

in an undistorted experiment, PAGEOPH, 164, 617-631.

Musumeci, R., I. A. Svendsen, J. Veeramony. (2005). The flow in the surf zone: a fully 

nonlinear Boussinesq-type of approach, Coastal Eng., 52 (7), 565-598.

NGDC. www.ngdc.noaa.gov/mgg/filters/01mgg04.html (Accessed: 04/16/2015)

NRC NUREG/CP-0302 (2013). Proceedings of the Workshop on Probabilistic Flood Hazard

Assessment (PFHA): Held at the U.S. Nuclear Regulatory Commission Headquarters,

Rockville, MD, January 29-31, 2013.

Nwogu, O. (1993). Alternative form of Boussinesq equations for nearshore wave propagation, J. 

Waterway, Port, Coastal and Ocean Eng., 119(6), 618-638.

Okada, Y., (1985). Deformation due to shear and tensile faults in a half space, Bull. Seism. Soc. 

Am., 75, 1135-1154.

Özkan-Haller, H.T., Kirby J.T. (1997). A Fourier-Chebyshev collocation method for the shallow 

water equations including shoreline run-up, Applied Ocean Research, 1997; 19: 21-34.

Pedersen, G. (2004). On Long Wave Runup Models, Proceedings of the 3rd International 

Workshop on Long-Wave Runup Models, June 17-18 2004, Catalina Island, California.

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

Pedrozo-Acuña, A., D. J. Simmonds, A. K. Otta, A. J. Chadwick (2006). On the cross-shore 

profile change of gravel beaches, Coastal Eng., 53(4), 335-347.

Peregrine, D.H. (1967). Long waves on a beach, J. Fluid Mech., 27 (4): 815-827. doi:

10.1017/S0022112067002605.

Petera J, Nassehi V. (1996). A new two-dimensional finite element model for the shallow water 

equations using a Lagrangian framework constructed along fluid particle trajectories, Int. J. 

Num. Meth. Eng., 39 (24), 4159-4182.

Power, W., X. Wang, E. Lane, P. Gillbreand (2012). A probabilistic tsunami hazard study of the 

Auckland region, Part I: propagation modeling and tsunami hazard assessment at the shoreline, 

PAGEOPH, 170, 1621-1634.
12-8

Priest, G.R., et al. (1997). Cascadia Subduction Zone Tsunamis: Hazard Mapping at Yaquina 

Bay, Oregon. Final Technical Report to the National Earthquake Hazard Reduction Program, 

DOGAMI Open File Report, 0-97-34, 144pp.

Ramsden, J.D. (1996). Forces on a Vertical Wall due to Long Waves, Bores, and Dry-Bed 

Surges, J. Waterway, Port, Coast. Ocean Eng., 122 (3), 134-141.

Raubenheimer, R.T. Guza, S. Elgar (1996). Wave transformation across the inner surf zone, J. 

Geophys. Res., 101 (C10), 25589-25597.

Ryu, S., Kim, M.H., and Lynett, P. (2003). Fully Nonlinear Wave-Current Interactions and 

Kinematics by a BEM-based Numerical Wave Tank, Computational Mechanics, 32, 336-346.

Science Review Working Group. (2007). Scientific and Technical Issues in Tsunami Hazard 

Assessment of Nuclear Power Plant Sites, NOAA Technical Memorandum OAR PMEL-136.

Seabra-Santos, F. J., D. P. Renouard, A. M. Temperville (1987). Numerical and experimental 

study of the transformation of a solitary wave over a shelf or isolated obstacles, J. Fluid Mech.,

176, 117-134.

Shen, M. and Meyer, R. (1963). Climb of a bore on a beach Part 3. Run-up, J. Fluid Mech., 16 

(1), 113-125.

Shi, F., Kirby, J. T., Harris, J. C., Geiman, J. D. and Grilli, S. T. (2012). A high-order adaptive 

time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation, 

Ocean Modelling, 43-44, 36-51.

Sielecki, A. and Wurtele, M.G., (1970). The numerical integration of the nonlinear shallow-water 

equations with sloping boundaries, J. Comput. Phys., 6 (2), 219-236.

Sitanggang, K., Lynett, P. (2009). Multi-scale Simulation with a Hybrid Boussinesq-RANS 

Hydrodynamic Model, Int. J. Num. Meth. Fluids, 62 (9), 1013-1046., doi: 10.1002/fld.2056.

Sittanggang, K., Lynett, P., and Liu, P., (2006). Development of a Boussinesq-RANS VOF 

Hybrid Wave Model, Proceedings of 30th ICCE, San Diego, 24-25.

Sitanggang, K. and Lynett, P. (2005). Parallel Computation of a Highly Nonlinear Boussinesq 

Equation Model through Domain Decomposition. Int. J. Num. Meth. Fluids, 49(1), 57-74.

Smith, J.M. (2007). Full-plane STWAVE with bottom friction: II. Model overview, CHETN-I-75. 

Vicksburg, MS: U.S. Army Engineer Research and Development Center,

http://chl.erdc.usace.army.mil/chetn.

Son, S., Lynett, P., and Kim, D.-H. (2011). Nested and Multi-Physics Modeling of Tsunami 

Evolution from Generation to Inundation, Ocean Modelling, 38 (1), 96-113, 

doi:10.1016/j.ocemod.2011.02.007.

Schwing, F.B., Norton, J.G., Pilskaln, C.H. (1990). Earthquake and bay: Response of

Monterey Bay to the Loma Prieta Earthquake, Eos, Transactions, AGU, 71. p. 402, 

doi:10.1029/90EO00035.
12-9

Stefanakis, T. S., Contal, E., Vayatis, N., Dias, F., & Synolakis, C. E. (2014). Can small islands 

protect nearby coasts from tsunamis? An active experimental design approach. In Proc. R. 

Soc A, 470, 20140575. 

Suppasri, A., Shuto, N., Imamura, F., Koshimura, S., Mas, E., & Yalciner, A. C. (2013). 

Lessons learned from the 2011 Great East Japan tsunami: performance of tsunami 

countermeasures, coastal buildings, and tsunami evacuation in Japan. Pure and Applied 

Geophysics, 170(6-8), 993-1018.

Synolakis, C.E., Imamura, F., Tsuji, Y., Matsutomi, S., Tinti, B., Cook, B., and Ushman, M. 

(1995). Damage, Conditions of East Java tsunami of 1994 analyzed, EOS, Transactions, AGU,

76, (26), 257 and 261-262.

Synolakis, C.E. (1987). The Runup of Solitary Waves, J. Fluid Mech,, 185, 523-545.

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

Lett., 77, 2141-2144

Tadepalli, S. and Synolakis, C. (1994). Runup of N-Waves on sloping beaches, Proc. Royal

Soc. London., 445, 99-112.

Tang, L., V. V. Titov, C. D. Chamberlin (2009). Development, testing, and applications of sitespecific tsunami inundation models for real-time forecasting, J. Geophys. Res., 114, C12025,

doi:10.1029/2009JC005476.

Tao, J. (1984). Numerical modeling of wave runup and wave breaking, Acta Oceanol., Sin. 6 

(5), 692-700 (in Chinese).

Tao, J. (1983). Computation of wave runup and wave breaking, Int. Report, Danish Hydraulics 

Institute, Denmark.

Thio, H.K. and Li, W. (2015). Probabilistic tsunami hazard analysis of the Cascadia Subduction 

Zone and the role of epistemic uncertainties and aleatory variability, In 11th Canadian 

Conference on Earthquake Engineering, Victoria, BC, 21–24 July 2015.

Thio, H.K., P. Somerville, and J. Polet (2010). Probabilistic tsunami hazard in California, PEER 

report 2010.108. Pacific Earthquake Engineering Research Center, College of Engineering, 

University of California, Berkeley.

Thornton, E.B., and Guza, R.T. (1983). Transformation of wave height distribution, J. 

Geophys. Res., 88, 5925-5938.

Titov, V. V., Moore, C.W., Spillane, M., Wei, Y., Gica, E., and Zhou, H. (2016). Tsunami hazard 

assessment based on wave generation, propagation, and inundation modeling for the U.S. East 

Coast, NUREG/CR-7222, NOAA, PMEL.

Titov, V.V., C. Moore, D.J.M. Greenslade, C. Pattiaratchi, R. Badal, C.E. Synolakis, and 

U.Kanoglu (2011). A new tool for inundation modeling: Community Modeling Interface

forTsunamis (ComMIT), PAGEOPH, 168(11), doi: 10.1007/s00024-011-0292-4, 2121-2131.
12-10

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

Port, Coastal and Ocean Eng., ASCE, 124(4), 157-171.

Titov, V. V. and Synolakis, C. E. (1995). Modeling of Breaking and Nonbreaking Long Wave 

Evolution and Runup using VTCS-2, J. Harbors, Waterways, Port, Coastal and Ocean Eng., 121 

(6), 308-316.

Tolkova, E., and W. Power (2011). Obtaining natural oscillatory modes of bays and harbors via 

Empirical Orthogonal Function analysis of tsunami wave fields, Ocean Dynam., 61(6),

doi:10.1007/s10236-011-0388-5, 731-751.

Tomita, T. and Honda, K., (2007). Tsunami estimation including effect of coastal structures and 

buildings by 3D model, Coastal Structures, 07, Venice.

Toro, E. F. (1999). Riemann Solvers and Numerical Methods for Fluid Dynamics, SpringerVerlag.

U.S. Nuclear Regulatory Commission (2013). Guidance for performing a tsunami, surge, or seiche 

hazard assessment (JLD-ISG-2012-06), ADAMS Accession No. ML12314A412.

Uslu, B., V.V. Titov, M. Eble, and C. Chamberlin (2010). Tsunami hazard assessment for

Guam. NOAA OAR Special Report, Tsunami Hazard Assessment Special Series, 1, 186 pp.

Veeramony, J., & Svendsen, I. A. (2000). The flow in surf-zone waves. Coastal Engineering, 

39(2), 93-122.

Ward, S.N. and S. Simon, (2005). Tsunami thoughts, Canadian Society of Exploration

Geophysicists (CSEG) Recorder, 30(10), 38-44.

Ward, S. N., and S. Day, (2001). Cumbre Vieja Volcano-- Potential collapse and tsunami at 

La Palma, Canary Islands, Geophys. Res. Lett., 28(17), 3397-3400.

Watts, P., Grilli, S.T., Kirby, J.T., Fryer, G.J., Tappin, D.R. (2003). Landslide tsunami case 

studies using a Boussinesq model and a fully nonlinear tsunami generation model, Nat. Hazards 

Earth Syst. Sci.,3, 391-402.

Wei, G., Kirby, J. T., Grilli, S. T. and Subramanya, R. (1995). A fully nonlinear Boussinesq 

model for surface waves, Part I. Highly nonlinear, unsteady waves, J. Fluid Mech., 294, 71-92.

Weiss, R., Wunnemann, K., and Bahlburg, H., (2006). Numerical modelling of generation, 

propagation, and run-up of tsunamis caused by ocean impacts: model strategy and technical 

solutions, Geophys. J. Int., 67, 77-88.

Witting, J. M., (1984). A unified model for evolution of nonlinear water waves, J. Comput.

Phys., 56, 203-236.

Woo, S-B., and Liu, P. L.-F. (2004). A finite element model for modified Boussinesq equations. 

Part I: Model development, J. Waterway, Port, Coastal and Ocean Eng., 130(1), 1-16.

Working Group on California Earthquake Probabilities (2003). Earthquake Probabilities in the 

San Francisco Bay Region: 2002-2031. U.S. Geological Survey Open File Report 03-214, 235p.
12-11

Wornom„ S., David J.S. Welsh, Keith W. Bedford., (2001). On Coupling the SWAN and WAM 

Wave Models for Accurate Nearshore Wave Predictions, Coas. Eng. J., 43, 161. doi: 

10.1142/S0578563401000335

Yamazaki, Y., Cheung,K.F., and Kowalik, Z. (2010). Depth-integrated, non-hydrostatic model 

with grid nesting for tsunami generation, propagation and runup, Int. J. Num. Meth. Fluids, 

doi:10.1002/ftd.2485.

Yeh, H. (2006). Maximum Fluid Forces in the Tsunami Runup Zone, J. Waterway, Port, Coastal 

and Ocean Eng., 132(6), 496-500.

Yeh, H., P. Liu, M. Briggs, and C. Synolakis (1994). Propagation and Amplification of Tsunamis 

at Coastal Boundaries, Nature, 372, 353-355.

Yeh, H., Ghazali, A., and Marton, I. (1989). Experimental study of bore runup. J. Fluid Mech., 

206, 563-578.

Zhang, Y., and Baptista, A.M. (2008). A semi-implicit Eulerian-Lagrangian finite element model 

for cross-scale ocean circulation, Ocean Modeling, 21(3-4), 71-96.

Zelt, J.A. (1991). The run-up of nonbreaking and breaking solitary waves, Coastal Engineering, 

15(3), 205-246

13-1

13. GLOSSARY

Bore – a wave with a steep vertical front. A tsunami wave may form a bore as it approaches 

shore.

Cnoidal Wave – is a nonlinear and periodic surface gravity wave with wavelengths large 

compared to the water depth. The sharp crests and very flat troughs are characteristic of these 

waves.

Courant Number – a measure of how much information traverses a computational grid cell in a 

given time-step.

Finite Difference Method (FD) – a numerical solution approach where a differential equation is 

solved at discrete points, typically arranged in a rectangular grid. The derivative operators are 

expressed with finite differences, which are Taylor Series based approximations, and can exist 

at any accuracy order and taken in any direction. FD methods are often the easiest to 

construct, but can lack flexibility in gridding.

Finite Volume Method (FV) – a numerical solution approach where a differential equation is 

solved within a discrete area or volume. Various FV methods exist, with their variations 

dependent on how the physical variables are averaged in or interpolated across the discrete 

volumes. FV methods can employ irregular grids and have robust solution techniques. 

Finite Element Method (FE) – a numerical solution approach where a differential equation is 

solved at arbitrarily located nodes, which are connected to form elements. Various techniques 

are available to approximate the governing equations across elements, and rely on 

interpolations between nodes and within elements. The strength of FE methods is that they are 

applicable for problems requiring complex gridding.

Manning Formula – is an empirical formula estimating the average velocity of a liquid flowing in 

an open channel, is driven by gravity.

Pade Approximation – approximation to a function by expanding it as a ratio of two power 

series and determining both the numerator and denominator coefficients.

Runup or run-up — Vertical difference between the elevation of tsunami inundation and the 

sea level at the time of a tsunami. Runup is the elevation of the highest point of land inundated 

by a tsunami as measured relative to a stated datum, such as mean sea level.

Short-term Inundation Forecasting for Tsunamis (SIFT) — A tsunami forecast system that 

integrates tsunami observations in the deep ocean with numerical models to provide an 

estimate of tsunami wave arrival and amplitude at specific coastal locations while a tsunami 

propagates across an ocean basin.

Shoaling – The effect by which surface waves entering shallower water increase in wave height

to maintain a constant energy flux.

Soliton – a wave in which the effects of dispersion and nonlinearity cancel each other, resulting 

in the wave propagating as a single pulse, maintaining its shape.

Technical 

1 2 1 2 2



NUREG/CR-7223 Tsunami Hazard Assessment: Best Modeling Practices and State-of-the-Art Technology December 2016