UNITED STATES

DEPARTMENT OF COMMERCE

Otto J. Wolf

Acting Secretary

NATIONAL OCEANIC AND

ATMOSPHERIC ADMINISTRATION

Mary M. Glackin

Acting Undersecretary for Oceans

and Atmosphere/Administrator

Office of Oceanic and 

Atmospheric Research

Richard W. Spinrad

Assistant Administrator

NOAA Technical Memorandum OAR PMEL-144

EXTRACTING TSUNAMI SOURCE PARAMETERS

VIA INVERSION OF DART® BUOY DATA

Donald D. Percival1

Diego Arcas2,3

Donald W. Denbo2,3

Marie C. Eble3

Edison Gica2,3

Harold O. Mofjeld3,4

Mick C. Spillane2,3

Liujuan Tang2,3

Vasily V. Titov2,3

1

Applied Physics Laboratory

University of Washington, Seattle, WA

2

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

University of Washington, Seattle, WA

3

NOAA/Pacific Marine Environmental Laboratory

Seattle, WA

4

Retired from NOAA

Pacific Marine Environmental Laboratory

Seattle, WA

February 2009
NOTICE from NOAA

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

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

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

opinions, findings, and conclusions or recommendations expressed in this material are those

of the authors and do not necessarily reflect the views of the National Oceanic and Atmospheric Administration.

Contribution No. 3134 from NOAA/Pacific Marine Environmental Laboratory

Also available from the National Technical Information Service (NTIS)

(http://www.ntis.gov)

ii
Contents iii

Contents

Abstract 1

1 Introduction 1

2 DART® Buoy Data for 2006 Kuril Islands Tsunami Event 2

3 Models for DART® Buoy Data 3

4 Fitting a Single Model to Data from a Single Buoy 5

5 Assessing Sampling Variability in Fitted Model 8

6 Selecting Amount of Data to Use for Estimating α 10

7 Handling Multiple Unit Sources and Buoys 12

7.1 Multiple Unit Sources and a Single Buoy . . . . . . . . . . . . 13

7.2 Multiple Unit Sources and Multiple Buoys . . . . . . . . . . . 17

8 Summary and Discussion 20

9 Acknowledgments 22

10 References 22

Extracting tsunami source parameters via inversion of DART®

buoy data

D.B. Percival,1, D. Arcas2,3, D.W. Denbo2,3, M.C. Eble3, E. Gica2,3, H.O. Mofjeld2,3,4,

M.C. Spillane2,3, L. Tang2,3, and V.V Titov2,3

Abstract. The ability to accurately forecast the hazards to coastal communities from earthquakegenerated tsunamis requires the knowledge of certain tsunami source parameters. The Short-term

Inundation Forecast for Tsunamis (SIFT) tool developed by the National Oceanic and Atmospheric

Administration estimates these source parameters based upon a combination of observed data from

DART® buoys and a set of precomputed models. We give a detailed description of the procedures

within SIFT Version 2.0 for estimating the source parameters and for assessing the sampling variability in the estimated parameters. We illustrate these procedures using detided 1-min data collected

by two DART® buoys (21414 and 46413) in the northwestern Pacific during the 15 November 2006

Kuril Islands tsunami event, generated by an offshore moment magnitude 8.3 earthquake. We discuss the advantages of using a least squares (L2) method for fitting model tsunami time series to

the data. The model series originate from two model source segments (50 km by 100 km) near the

earthquake epicenter. The fitting procedure estimates the source parameters of each segment, both

individually and together. Realistic confidence limits on the source parameters are obtained by fitting the residual series (from the data-fitted model) to a first-order autoregressive model to account

for the obvious correlation between adjacent residual values. Using data from both buoys improves

the source parameter estimates, at an operational cost of some delay in time for the tsunami to

reach the farther station. An appropriate fit uses the first full tsunami wave at each station with

no improvement found by extending the fitted data segment further. Formulae are derived for estimating the source parameters and computing the confidence limits for multiple DART® buoys and

source segments.

1. Introduction

The National Oceanic and Atmospheric Administration (NOAA) is charged

with primary responsibility for providing tsunami warnings to the United

States and has assumed a leadership role both in research on tsunamis and

in gathering data related to tsunamis (see http://www.tsunami.noaa.gov/

for details on NOAA’s tsunami program). The ability to accurately forecast

potential dangers to coastal communities due to an earthquake-generated

tsunami requires knowledge of tsunami source parameters. The Short-term

Inundation Forecast for Tsunamis (SIFT) tool that has been developed at

the NOAA Center for Tsunami Research extracts estimates of these parameters based upon two key components. The first is observed data from

DART® buoys collected during a tsunami event, and the second is a set

of precomputed models for each DART® location. The model series are

extracted from the numerical solution to the propagation of tsunami waves

throughout the ocean basin from “unit sources” and, because of the dynam1Applied Physics Laboratory, University of Washington, Box 355640, Seattle, WA

98195-5640 2Joint Institute for the Study of the Atmosphere and Ocean (JISAO), University of

Washington, Seattle, WA 98195 3NOAA/Pacific Marine Environmental Laboratory, Seattle, WA 98115

4Retired from NOAA
2 Percival et al.

ics of these waves in the open ocean, can be linearly combined to mimic the

observed data. An inversion procedure combines these key components by

using the DART® buoy data to extract parameters that adjust the amplitudes of the precomputed models. These parameters, once determined and

applied to the basin-wide set of “unit source” solutions, provide the boundary conditions for inundation models that generate high-resolution forecasts

of incoming tsunami waves at threatened coastal communities (Gica et al.,

2009). The intent of this memorandum is to describe the key steps in inversion procedures in general and to explain the rationale behind the particular

procedure implemented in SIFT, which is being installed at the U.S. Tsunami

Warning Centers.

The remainder of the memorandum is organized as follows. To motivate

our discussion, we use data recorded by two DART® buoys in November

2006 that include a tsunami event originating from an earthquake near the

Kuril Islands (see Horrillo et al. (2008); Kowalik et al. (2008) for details

about this event). We describe the data in Section 2, after which we consider precomputed models for the observed data (Section 3). For simplicity,

we first focus on a single buoy and a single model and discuss several criteria

for fitting the model to the data (Section 4). Use of a least squares criterion has several advantages, not the least of which is the ability to assess

the sampling variability in the estimated parameters (Section 5). We then

consider the important question of how much data should be used to carry

out the inversion (Section 6), after which we expand our methodology to

allow for more than one buoy and one model (Section 7). We make some

concluding remarks in Section 8.

2. DART® Buoy Data for 2006 Kuril

Islands Tsunami Event

The map in Fig. 1 shows the location of the source of the Kuril Islands

tsunami event of 15 November 2006, along with the locations of two DART®

buoys that recorded the event, namely, 21414 and 46413 (referred to hereafter as buoys 1 and 2). Figure 2 shows plots of the data taken during the

event by these two buoys. The data in these plots have been subjected to

a detiding procedure. The time associated with each detided data point is

in hours since the beginning of the event, which was at 11:14:16 UTC. Note

that the event was observed at buoy 2 about half an hour after it occurred at

buoy 1, which is consistent with their locations relative to the Kuril Islands

shown in Fig. 1. The spacing in time between adjacent values is 1 min, and

each data point represents a 1-min average. There are 105 data points for

buoy 1 and 114 for buoy 2. In Sections 3 to 6 we concentrate on the data

and models for buoy 1 alone, after which we consider both buoys. In what

follows, let T

be the set of times associated with the data from buoy 1, and,

for a given t ∈ T

, let xt represent the corresponding detided data.
Extracting tsunami source parameters via inversion of DART® buoy data 3

Figure 1: Map showing the approximate origin of the 15 November 2006 Kuril Islands tsunami

event, which was generated by a moment magnitude 8.3 earthquake at 11:14:16 UTC. The locations

of DART® buoys 21414 and 46413 (henceforth referred to as buoys 1 and 2) are also shown, along

with a grid of 50 km by 100 km “unit sources.” Two of these unit sources are marked as A12 and B13

(henceforth referred to as sources a and b) and are close to the epicenter of the earthquake (denoted

by a star on the map). For both unit sources, models of what would be observed at the two buoys

were precomputed under the assumption that the generating event is equivalent to a magnitude 7.5

reverse thrust earthquake located within the unit source.

3. Models for DART® Buoy Data

Here we consider models for the detided DART® data from buoy 1 (see Gica

et al., 2008, for details about how these models were formulated). The earthquake that generated the Kuril Islands tsunami event is presumed to have

originated from one or more unit sources, which are defined by nominal 50 km

by 100 km rectangles located along a fault line. We will consider initially just

one of these unit sources marked on Fig. 1 as “A12” and henceforth referred

to as unit source a. The circles in Fig. 3 show a precomputed model for what

we should expect to observe at buoy 1 under the assumption that there was

a magnitude 7.5 reverse thrust earthquake of appropriate strike, dip, and

depth originating from unit source a. The model was run with a 15-sec time

step, but, to save space on the operational database, was subsampled down
4 Percival et al. −0.06 0.00 0.06 meters

buoy 1 (21414)

1.5 2.0 2.5 3.0 3.5

hours since earthquake

−0.06 0.00 0.06

meters

buoy 2 (46413)

Figure 2: Detided data from 15 November 2006 Kuril Islands tsunami event as derived from data

recorded by two DART® buoys.

1.5 2.0 2.5 3.0 3.5

hours since earthquake

−0.005 0.005

meters

Figure 3: Model for data observed at buoy 1 generated assuming a magnitude 7.5 reverse thrust

earthquake originating from unit source a (i.e., A12).

to a discrete grid of times with a spacing of 1 min. In general, the times

used in a precomputed model might or might not correspond to the times at

which the DART® buoy data were actually collected. To facilitate matching

up the model and the observed data, we use a cubic spline to interpolate the

model, the result of which is shown by the solid curve passing through the

circles in Fig. 3. We let g(t) represent the spline-interpolated model in what

follows.
Extracting tsunami source parameters via inversion of DART® buoy data 5

1.5 2.0 2.5 3.0 3.5

hours since earthquake

−0.05 0.00 0.05

meters

Figure 4: Comparison of detided data from buoy 1 (dots) and its model from unit source a (solid

curve, as also shown in Fig. 3). In the inversion process the model will be rescaled to match the

observed amplitude of the detided data.

4. Fitting a Single Model to Data from a

Single Buoy

Figure 4 compares the precomputed model g(t) from unit source a (the solid

curve) and the detided data xt for buoy 1 (the dots). The model and the data

look different because the model was computed for a magnitude 7.5 reverse

thrust earthquake. We assume that an earthquake with an arbitrary magnitude can be handled by multiplying g(t) by a parameter α. Accordingly

we entertain the model

xt = αg(t) + et, t ∈ T

,

where et is a residual term that represents any remaining mismatch between

the data and the amplitude-adjusted model.

What criterion should we use to set the parameter α? For any proposed

α, we can compute a corresponding set of residuals et = xt −αg(t) for t ∈ T

.

As examples, let’s compute residuals corresponding to the choices α = 2,

4. and 8. The adjusted models αg(t) and data are shown in the left-hand

column of Fig. 5, while the corresponding residuals et are displayed in the

right-hand column (the residual with the largest absolute value is marked

with a circle). Arguably the “best” choice for α is such that the residuals et

are “small” by some measure. There are many possible measures, but three

common ones are to

1. make the sum of the squared residuals as small as possible:



t∈T

e2

t = 

t∈T

[xt − αg(t)]2 ≡ f2(α);
6 Percival et al. −0.05 0.00 0.05 meters

α = 2 data residuals

−0.05 0.00 0.05

meters

α = 4

1.5 2.0 2.5 3.0 3.5

hours since earthquake

−0.05 0.00 0.05

meters

α = 8

1.5 2.0 2.5 3.0 3.5

hours since earthquake

Figure 5: Adjusted model αg(t) with α = 2, 4, and 8 (solid curves in left-hand column of plots),

detided data from buoy 1 (dots in left-hand column) and corresponding residuals et (dots in righthand column). The circle in each right-hand plot indicates the residual with the largest magnitude.

The vertical dashed and dotted lines delineate the beginning and end of the first complete wave from

unit source a.

2. make the sum of the magnitudes of the residuals as small as possible:



t∈T

|et| = 

t∈T

|xt − αg(t)| ≡ f1(α);

3. and make the largest magnitude of the residuals as small as possible:

max

t∈T|et| = maxt∈T|xt − αg(t)| ≡ f∞(α).

These measures correspond to what are known in the literature as the L2,

L1, and L∞ norms. A specialized measure of some interest is to

4. make the sum of the squared residuals at the model peak and trough

as small as possible:

e2

t0 + e2t1 = [xt0 − αg(t0)]2 + [xt1 − αg(t1)]2 ≡ fpt(α),
Extracting tsunami source parameters via inversion of DART® buoy data 7 residual measure 0.030 0.040

f2(α)

1.2 1.3 1.4 1.5

f1(α)

0 2 4 6 8 10

α

0.045 0.055

residual measure

f∞(α)

0 2 4 6 8 10

α

0.001 0.003 0.005

fpt(α)

Figure 6: Residual measures f2(α), f1(α), f∞(α) and fpt(α) versus parameter α. The vertical

dashed lines show the location of the minima for each measure.

where t0 ∈ T

and t1 ∈ Tare such that g(t0) = maxt∈T{g(t)} and g(t1) =

mint∈T{g(t)}. The above can be expressed more compactly as

fpt(α) = 

t∈{t0,t1}

e2

t = 

t∈{t0,t1}

[xt − αg(t)]2 ,

which points out the similarity between this measure and f2(α), both of

which are based on squared residuals and are commonly referred to as least

squares criteria. In what follows, it will be convenient to let fls(α) stand for

either f2(α) or fpt(α) and to write

fls(α) = 

t

e2

t = 

t

[xt − αg(t)]2 ,

where it is understood that t takes on values in the set T

for f2(α) and in

{t0, t1} for fpt(α).

Figure 6 shows plots of f(α) versus α for the four measures, where we

vary α from 1 to 10 in steps of 0.2, a total of 46 different settings in all.

For any given measure, the estimated parameter is the value of α among

the 46 settings for which f(α) is minimized. The estimated parameters

corresponding to f2(α), f1(α), f∞(α) and fpt(α) are ˆα2 = 4.2, ˆα1 = 5.0,

αˆ∞ = 3.8 and ˆαpt = 4.2. These estimates are indicated on Fig. 6 by vertical

dashed lines.

While all four measures give estimated parameters that are roughly the

same, the two measures based on least squares criteria have certain advantages. First, the value of α minimizing fls(α) is given by a simple formula,
8 Percival et al.

namely,

αˆls =



t xtg(t)



t g2(t) (1)

(this follows by setting the derivative of fls(α) with respect to α to zero

and solving the resulting equation for α). This formula yields the more

refined estimates ˆα2 .= 4.17 and ˆαpt .= 4.26, without the need to locate the

minimizing α using an exhaustive grid search. Second, as we show in the next

section, we can quantify the statistical variation in ˆαls readily once we make

some assumptions about the nature of the residuals et. Similar assessments

for either ˆα1 or ˆα∞ are not as easy to come by. Third, the f1(α) measure is

known to be less influenced by large data values than least squares measures.

De-emphasizing large data values is a disadvantage in the current context

because they are arguably the best indicators in the measured data of the

magnitude of the tsunami event. Fourth, while f∞(α) is more influenced by

large data values than least squares measures, it can be unduly influenced

by a single rogue data point. Given these considerations, we will focus on

the least squares measures f2(α) and fpt(α) henceforth.

5. Assessing Sampling Variability in

Fitted Model

Here we consider how we can quantify the sampling variability in a least

squares estimator of the parameter α. To facilitate our discussion, let us

reformulate the model xt = αg(t) + et in vector notation as

x ̃ = αg ̃ + e ̃,

where x ̃, g ̃, and e ̃ are column vectors whose tth elements are, respectively,

xt, g(t), and et. The dimension Nof these vectors is the number of elements

in T

for f2(α) and is two for fpt(α). The least squares estimator of αls in

(1) can be rewritten as ˆαls = g ̃T x ̃/g ̃T g ̃, where g ̃T denotes the transpose of

g ̃.

In order to quantify the sampling variability in ˆαls, we need to consider

the statistical properties of the residuals et. A common assumption is that

the residuals are a sample from a collection of independent and identically

distributed (IID) random variables from a Gaussian (normal) distribution

with mean zero and variance σ2

e . With this assumption, a standard statistical

result says that ˆαls is Gaussian distributed with mean α and with a variance

given by σ2

e/g ̃T g ̃ (see, e.g., Draper and Smith, 1998). Let eˆ ≡ x ̃ − αˆlsg ̃ be a

vector containing the observed residuals. An unbiased estimator of σ2

e is

σˆ2

e = eˆT eˆ

N− 1

.

The estimated standard deviation for ˆαls is the square root of ˆσ2

e/g ̃T g ̃. Using

this estimate, we can express the uncertainty in the least squares estimates
Extracting tsunami source parameters via inversion of DART® buoy data 9

as ˆα2 = 4.17 ± 0.59 and ˆαpt = 4.26 ± 2.69, with the interpretation that the

corresponding intervals [3.58, 4.76] and [1.57, 6.95] represent 68% confidence

intervals for the unknown parameter α. Note that the sampling variability

in ˆαpt is considerably larger than that of ˆα2, presumably a reflection of the

fact that ˆα2 is based on significantly more data than is ˆαpt (105 data points

versus 2).

Unfortunately the IID assumption is highly questionable for the data

used to estimate ˆα2. A plot of the observed residuals eˆ associated with this

estimate looks virtually the same as the α = 4 case shown in Fig. 5 (middle

plot, right-hand column). The IID assumption dictates that the residuals

have the same variance, but note that those associated with the data prior to

the tsunami arrival clearly have less variability than elsewhere. Data prior

to the wave arrival are associated with portions of model for which g(t) is

practically zero. As can be seen from (1), these data values do not influence

the estimate ˆα2 since the products xtg(t) are also close to zero. For values of

t corresponding to the actual tsunami, the observed residuals appear to have

the same variability to a viable approximation, but they cannot reasonably

be assumed to be independent: residuals that are adjacent to one another

are quite similar, exhibiting a significant degree of autocorrelation.

We can adjust our model to account for the obvious violations of the IID

assumption by considering the reduced model

x = αg + e, (2)

where the vector x of dimension N = 70 is the portion of x ̃ that contains just

the data starting from the earliest time at which g(t) differs substantively

from zero, with g and e being similarly extracted from g ̃ and e ̃. The least

squares estimator of α is now expressed as ˆα2 = gT x/gT g. To account for

the apparent autocorrelation in the observed residuals, we assume a simple

model, namely, that the residuals obey a first-order stationary autoregressive

process et = φet−1+wt, where |φ| < 1, and the wt variables are IID Gaussian

deviates with mean zero and a common variance. This model implies that

the correlation between residuals et and et+τ separated by τ minutes is φ|τ|

.

Let us redefine the vector of observed residuals to be eˆ ≡ x − αˆ2g. Since φ

is the correlation between adjacent residuals et and et+1, we can estimate it

via

φˆ ≡

eˆT

(f)eˆ(l)

eˆT eˆ

.

= 0.86, (3)

where eˆ(f) consists of all of eˆ except for its first element, and eˆ(l) has everything but the last element. Statistical theory says ˆα2 is Gaussian distributed

with a mean of α and a variance given by σ2

e gT V g/(gT g)2, where V is a matrix whose (j, k)th element is φ|j−k|

. An approximately unbiased estimator

of σ2

e is given by

σˆ2

e = (1 − φ)2eˆT eˆ

N(1 − φ)2 − (1 − φ2)+2φ(1 − φN )N −1 . (4)

We can thus estimate the standard deviation of ˆα2 by taking the square

root of ˆσ2

e gT V g/(gT g)2 after replacing φ by φˆ. We now have ˆα2 = 4.17 ±
10 Percival et al.

1.50, which reports a larger uncertainty than under the questionable IID

assumption, namely, ˆα2 = 4.17 ± 0.59. (The procedure for getting a more

realistic estimate of the uncertainty in ˆα2 breaks down if we attempt to

adjust it to work with ˆαpt. This estimator is based on just two data points,

and hence getting a good estimate of φ is problematic.)

6. Selecting Amount of Data to Use for

Estimating α

The inversion algorithm in Version 2.0 of SIFT uses the least squares estimator ˆα2 to estimate the parameter α based upon available DART® buoy

data. In the example in the previous section, we used all of the data relevant

to the Kuril Islands tsunami event observed by buoy 1. Using the estimated

parameter to forecast the potential dangers of a tsunami to coastal communities leads naturally to the question of how this estimate would be affected

if we were to use less data than what eventually became available. In theory,

using more data should yield a better estimate of the parameter, but this

must be balanced off with the necessity of obtaining an estimate as quickly

as possible. Here we use the example presented in the last section to illustrate how estimates of the parameter change as we use smaller amounts of

data.

The top plot of Fig. 7 shows the data (the dots) for the Kuril Islands

tsunami event recorded at buoy 1, along with the model (shown as a curve)

from unit source a after multiplication by the parameter ˆα2 = 4.17, which

is based on all the available data. There are two vertical dashed lines on the

plot, which mark the times t0 = 2.046 and t1 = 2.196. The first of these is

the earliest time at which the model g(t) is measurably different from zero;

the second, the time just following the crest of the first full tsunami wave

as measured at the buoy (i.e., just after the so-called quarter-wave point).

There are 10 data points in all between (and including) t0 and t1; similarly,

there are 70 data points between t0 and the last observed value (te = 3.196).

Various results related to estimating α using just the first 10 points and all 70

points are shown by, respectively, the left and right ends of the curves in the

bottom four plots. The rest of the curves shows results using 11, 12,... , 69

points. From top to bottom, the second plot in the figure shows the least

squares estimates ˆα2; the third, the estimated standard deviations for ˆα2;

the fourth, the root mean square (RMS) errors of the residuals over just the

data used to determine ˆα2 (solid curve) and over all 70 data points (dashed

curve, which necessarily corresponds exactly to the solid curve at the very

end); and the fifth, two versions of the so-called R2 statistic, which is the

percentage of the sample variance of the data explained by the model (the

solid curve is the squared correlation—expressed as a percentage—between

the observed data and the fitted model, while the dashed curve shows the

correlations under the assumptions that the true means of both the data

and the model are zero).
Extracting tsunami source parameters via inversion of DART® buoy data 11 −0.05 0.00 0.05 meters

α = 4.17

4.5 5.5 6.5

estimated 

α

123

sd(estimated 

α)

0.010 0.020

RMS error

1.5 2.0 2.5 3.0 3.5

hours since earthquake

0 50 100

R2

Figure 7: Detided data from buoy 1 (dots, top plot) and model αg(t) from unit source a, with α

set to the least squares estimate ˆα2 = 4.17 using all the data (solid curve, top plot). The bottom

four plots summarize the effect of using increasing amounts of data to estimate α, with the smallest

amount used being indicated by the vertical dashed lines in the top plot. From top to bottom, these

plots show ˆα, the estimated standard deviation for ˆα, the RMS error over the selected data (solid

curve) and all the data (dashed) and two variations on the R2 statistic. The left-most dashed line

and the dotted line in the top plot delineate the 24 data points recorded during the first complete

wave from unit source a; the dotted lines in the other plots mark values corresponding to an estimate

of α based on these 24 values.
12 Percival et al.

Some interesting patterns emerge as the amount of data used to estimate

α increases from 10 to 70. The second plot shows that the estimated parameters ˆα2 roughly decrease with increasing amounts of data, so the estimated

parameter does depend on the amount of data used. The third indicates

that the associated standard deviation of the estimator increases initially,

but then decreases as more data are used, as intuition would suggest. The

fourth shows that the root mean square errors over the data used in the

fit and over all the data are relatively constant after about t = 2.4. The

last plot shows that the percent variance explained is initially fairly high

(around 90%) when just 10 data points are used, but then decreases monotonically after about t = 2.5 with the use of more data, eventually ending up

at around 32%. This indicates that the model has little explanatory power

beyond the first complete wave observed at the buoy.

This Kuril Islands tsunami event is one in which the very first wave is the

largest when observed at buoy 1. This is because there is an unobstructed

path for the tsunami to propagate from the source of the event to this

buoy (see Fig. 1, where buoy 1 is marked as 21414). For this case, we

are thus better off just using the data up to the first complete wave (as

defined by model a) to estimate α since the data and the model disagree

substantially beyond that point. Doing so yields ˆα2 = 5.15 ± 2.36 based

upon N = 24 data points. There are other situations in which later waves

can be larger, in which case it would be desirable to use a longer stretch of the

data for estimating the parameter α. How best to select the amount of data

to use in the inversion procedure is an open question requiring additional

research. Version 2.0 of SIFT allows the user to manually select the portions

of DART® data to use from each buoy, but there are provisional plans for

future versions to offer guidance based upon statistical tests.

7. Handling Multiple Unit Sources and

Buoys

So far we have considered an inversion procedure that is based upon a model

from a single unit source for the data received at a single DART® buoy. In

practice, certain tsunami events are potentially better described by considering multiple source models. This possibility is illustrated in Fig. 1, where

sea floor disturbances are highlighted in several unit source rectangles that

might have contributed to the tsunami. In addition the same tsunami event

can be recorded by multiple buoys. Figure 2 reflects the fact that the Kuril

Islands event was clearly detected at at least two DART® buoys, thereby

providing additional observations for use in source determination. In the

subsections that follow, we describe extensions to our inversion procedure

that handle multiple source models and/or multiple buoys. We use unit

sources a and b and buoys 1 and 2 to illustrate how these extensions work.

For j = 1 or 2, we let xj,t, t ∈ Tj , represent the relevant data from buoy

j, and we let xj be the corresponding vector of length Nj , where Nj is the
Extracting tsunami source parameters via inversion of DART® buoy data 13

1.5 2.0 2.5 3.0 3.5

hours since earthquake

−0.005 0.005

meters

Figure 8: Model series for buoy 1 from unit sources a (thin curve, which is the same as the curve

in Fig. 3) and b (thick curve with dots) generated under the assumption that a magnitude 7.5

reverse thrust earthquake originated from each source alone. The dashed and dotted lines mark the

beginning and end of the first complete wave from source a.

number of elements in the set of indices Tj . With this new notation, we can

re-express the reduced model of (2) for buoy 1 as

x1 = αag1;a + e1, (5)

where the subscript on αa says that this parameter is associated with the

model from unit source a; the subscripts on g1;a indicate that this vector

contains the model for buoy 1 from source a; and the subscript on e1 says

that these are the residuals associated with buoy 1. We can now express the

least squares estimator of αa as

αˆ1;a = gT

1;ax1/gT1;ag1;a = 

gT

1;ag1;a

−1

gT

1;ax1 (6)

(cf. Equation (1)), where the subscripts on ˆα1;a emphasize that this estimate of the parameter for unit source a is based on data from buoy 1 (the

right-most expression above is of interest to compare with generalizations to

come).

7.1 Multiple Unit Sources and a Single Buoy

The event that generates a tsunami is sometimes better described as originating from more than just one unit source. When this is the case, we

assume that the overall effect of the generating event is a linear combination

of the individual models from each source. As an example, the thin curve

in Fig. 8 shows the model from unit source a for buoy 1 (shown previously

in Fig. 3), while the thick curve with superimposed dots shows the model

from source b. Note that the first quarter wave from b arrives later than

that from a, but that the duration of the first complete wave is shorter for b

than for a. If we compare the data from buoy 1 with the fitted model from

a in the top plot of Fig. 7, we see that the peak in the data occurs after the

peak in the model, and hence the peak location might be better represented
14 Percival et al.

by the model coming from b; on the other hand, the duration of the first

wave in the data is more in keeping with that exhibited by a than by b. We

might get a better overall fit if we postulate an overall model for the data

from buoy 1 as being a blend of the models coming from a and b.

Accordingly we now augment the model of Equation (5) to include an

additional term:

x1 = αag1;a + αbg1;b + e1,

where αb is the parameter for unit source b, and g1;b is a vector of dimension

N1 containing the unadjusted model for buoy 1 from source b. As before,

we assume that the residual term e1 is multivariate Gaussian with a mean

vector of zeros and a covariance matrix given by σ2

1V1, where the (j, k)th

element of the N1 × N1 matrix V1 is given by φ|j−k|

1 for some |φ1| < 1. We

can express this model more compactly as

x1 = G1α + e1, (7)

where G1 is an N1 × 2 dimensional array whose first and second columns

are g1;a and g1;b, while α = [αa, αb]

T . The least squares estimator αˆ1 of α

satisfies the equation

GT

1 G1αˆ 1 = GT1 x1, (8)

where the subscript on αˆ1 emphasizes that the parameters are being estimated based upon the data from buoy 1. Under the assumption that the

2 × 2 matrix GT

1 G1 is invertible, the least squares estimator is given by

αˆ1 = 

GT

1 G1

−1

GT

1 x1,

which is similar in form to the estimator in Equation (6) when using unit

source a alone. Standard least squares theory says that αˆ1 is an unbiased

estimator of α and obeys a Gaussian distribution with a covariance matrix

given by

Σ1 ≡ σ2

1



GT

1 G1

−1

GT

1 V1G1



GT

1 G1

−1 (9)

(see, e.g., Draper and Smith, 1998 for details). Once we have determined

αˆ1, we can compute the observed residuals

eˆ1 = x1 − G1αˆ1

and use these to form an estimator of φ1 via the obvious analog of Equation (3). Conditional upon this estimate of φ1, we can obtain an approximately unbiased estimator of σ2

1 using the analog of Equation (4), following

which we have all the pieces needed to determine the elements of the covariance matrix Σ1.

Using the N1 = 24 data values collected by buoy 1 during the first

complete wave from source a, the least squares estimate of α is

αˆ 1 =

 αˆ1;a

αˆ1;b

 .

=

 4.24

3.00 

. (10)
Extracting tsunami source parameters via inversion of DART® buoy data 15 −0.05 0.00 0.05 meters

a only data residuals

−0.05 0.00 0.05

meters

b only

1.5 2.0 2.5 3.0 3.5

hours since earthquake

−0.05 0.00 0.05

meters

a and b

1.5 2.0 2.5 3.0 3.5

hours since earthquake

Figure 9: Comparison of fitted models involving unit sources a and b with data from buoy 1. The

left-hand plots show the data (dots) along with the models after adjustment using the estimated

source parameters (solid curves). The right-hand plots show the corresponding observed residuals.

The vertical dashed and dotted lines mark the stretch of 24 data points used to estimate the parameters. In the first row, only unit source a is used; in the second, only b; and the bottom row uses

both sources together.

The corresponding estimates of φ1 and σ2

1 are 0.89 and 0.00061, based upon

which we calculate

Σˆ 1 .

=

 3.43 −0.11

−0.11 2.19 

. (11)

The estimated parameters, along with their standard errors, for the two

sources are thus ˆα1;a .

= 4.24 ± 1.85 and ˆα1;b

.

= 3.00 ± 1.48.

The left-hand plot in the bottom row of Fig. 9 shows the fitted model

G1αˆ1 over T1 and its extension outside of T1 (solid curve), where the vertical

dashed and dotted lines indicate the earliest and latest times in T1. The data

from buoy 1 that were used to determine αˆ1 are shown as dots between the

vertical lines, and the remaining data, as the dots before and after these

lines. The right-hand plot shows the observed residuals eˆ1 (dots between

the vertical lines) and the differences between the observed data and the

extended fitted model (dots outside the vertical lines). The top row in
16 Percival et al.

1.5 2.0 2.5 3.0 3.5

hours since earthquake

−0.005 0.005

meters

Figure 10: As in Fig. 8, but now for buoy 2.

Fig. 9 shows corresponding plots when the data are modeled using just unit

source a, and the middle row, just source b. We see that the use of both

unit sources gives somewhat smaller observed residuals over the first quarter

wave of the observed data.

We can readily adapt the above procedure to handle more than two

unit sources. If we are interested in Ns sources in all, then G1 becomes

an N1 × Ns dimensional array whose kth column contains the unadjusted

model for buoy 1 from the kth source, while α becomes an Ns dimensional

column vector. Two problems, however, can crop up when dealing with

multiple unit sources. The first is that, due to dependencies amongst the

models from different unit sources, the matrix GT

1 G1 might not be of full

rank or might be close enough to being rank deficient so that computing the

inverse of GT

1 G1 becomes numerically unstable. We deal with this potential

problem by using a singular value decomposition in solving Equation (8)

for α to determine αˆ 1. The second problem is that some of the estimated

parameters might turn out to be positive, while others are negative. A

mixture of positive and negative parameters is difficult to reconcile with

the physics of earthquake generation. While use of such a mixture might

achieve a closer fit between model and observations, the application of these

parameters throughout the ocean basin can lead to unrealistic predictions

for inundations of coastal regions. We can deal with this problem using a

constrained least squares procedure, whereby, given x1 and G1, we minimize

x1 − G1α2 over α subject to either the condition α ≥ 0 or α ≤ 0,

i.e., all of the elements of α are either nonnegative or nonpositive (in the

above · denotes the Euclidean norm, and 0 is an Ns dimensional column

vector of zeros). This constrained least squares procedure can result in some

parameters being set to zero, which in effect removes the corresponding unit

sources from our model for the data. The above minimization problem is a

special case of Problem 10.1.1 of Fletcher (1987), and the method we use to

solve it is a variation of Algorithm 10.3.4 of Fletcher (1987).
Extracting tsunami source parameters via inversion of DART® buoy data 17

7.2 Multiple Unit Sources and Multiple Buoys

We now consider the data collected by buoy 2 (i.e., 46413), which is shown

in the bottom plot of Fig. 2. Figure 10 shows the unadjusted models for

these data from unit sources a and b. We can make use of these additional

data and associated models in three ways. First, we can use the least squares

estimates of the parameters from Equation (10) based upon the data from

buoy 1 to see how well it matches the data observed by buoy 2 (i.e., “crossvalidation”). The first row of Fig. 11 shows the data for buoy 2, its adjusted

model using unit sources a and b and the associated residuals (the layout

of the plots follows that of Fig. 9). The agreement between the model and

the data is relatively good, thus lending some additional credence to the

estimated parameters. Second, we can formulate a model that is analogous

to what we used for buoy 1 and estimate α based upon just the data from

buoy 2. Let x2 be a vector of dimension N2 = 19 containing the data from

buoy 2 that was sampled over the first complete wave in the model from

unit source a (this region is marked by the dashed and dotted vertical lines

in Fig. 10 and in the middle two rows of Fig. 11). The model parallels that

of Equation (7) and states that

x2 = G2α + e2, (12)

where G2 is an N2 ×2 dimensional array whose first and second columns are

g2;a and g2;b (vectors of dimension N2 containing the unadjusted models for

buoy 2 from sources a and b), and e2 is an N2 dimensional vector of error

terms that is multivariate Gaussian with zero mean and a covariance matrix

given by σ2

2V2 (the (j, k)th element of the N2 × N2 matrix V2 is given by

φ|j−k|

2 for some |φ2| < 1). The least squares estimate of α is

αˆ 2 =

 αˆ2;a

αˆ2;b

 .

=

 3.43

3.63 

.

The estimates of φ2 and σ2

2 are 0.95 and 0.00098, based upon which we

calculate the analog of Σˆ 1 in Equation (11), namely,

Σˆ 2 .

=

 2.49 0.19

0.19 1.94 

.

The estimated parameters are ˆα2;a .

= 3.43 ± 1.58 and ˆα2;b

.

= 3.63 ± 1.39,

which, when sampling variability is taken into account, are in reasonably

good agreement with what we obtained using just the data from buoy 1.

The second row of Fig. 11 has plots for this model that parallel those in the

top row.

A third use of the additional data from buoy 2 is to estimate α using

a model involving data x1 and x2 from both buoys. This joint model can

be formulated by “stacking” together the models given by Equations (7)

and (12):

x1:2 = G1:2α + e1:2, (13)
18 Percival et al. −0.05 0.00 0.05 meters

4.24, 3.00 buoy 2 data residuals

−0.05 0.00 0.05

meters

3.43, 3.63

−0.05 0.00 0.05

meters

3.80, 3.35

1.5 2.0 2.5 3.0 3.5

hours since earthquake

−0.05 0.00 0.05

meters

buoy 1 data

1.5 2.0 2.5 3.0 3.5

hours since earthquake

Figure 11: As in Fig. 9, with upper three rows comparing the data from buoy 2 to models based

upon unit sources a and b with parameters determined by data from buoy 1 alone (upper row), buoy

2 alone (second row) and both buoys (third). The two values given in the upper left-hand corners

are the estimated parameters involving unit sources a and b. The bottom plot compares the data

from buoy 1 to the model using the parameters determined by data from both buoys (i.e., the same

as in the third row). In the left-hand plots the vertical dashed and dotted lines indicate the starts

and ends of the data used to estimate the parameters (these lines are replicated in the right-hand

plots).
Extracting tsunami source parameters via inversion of DART® buoy data 19

where x1:2 is an N1 + N2 dimensional vector obtained by stacking x1 on top

of x2; G1:2 is an (N1 + N2) × 2 dimensional array obtained by stacking G1

on top of G2; and e1:2 is the residual vector; thus

x1:2 =

 x1

x2



and G1:2 =

 G1

G2



.

We assume that e1:2 is multivariate Gaussian with a mean vector of zeros

and an (N1 + N2) × (N1 + N2) covariance matrix given by

V1:2 =

 σ2

1V1 0

0T σ2

2V2



,

where the (j, k)th element of the N1 × N1 dimensional submatrix V1 is given

by φ|j−k|

1 while that for the N2 × N2 dimensional submatrix V2 is φ|j−k|2 ,

where |φ1| < 1 and |φ2| < 1 (in the above 0 and 0T are taken to represent an

N1×N2 matrix of zeros and its transpose). This model assumes that any two

error terms are uncorrelated if they are associated with different buoys (this

assumption is a reasonable starting point, but is subject to further research

pending examination of residuals from operational use of this model). The

least squares estimator αˆ1:2 of α satisfies the equation

GT

1:2G1:2αˆ1:2 = GT1:2x1:2, (14)

which is analogous to Equation (8). As with the models involving data from

a single buoy, αˆ1:2 is an unbiased estimator of α and follows a Gaussian

distribution with a covariance matrix given by

Σ1:2 ≡ 

GT

1:2G1:2−1

GT

1:2V1:2G1:2 

GT

1:2G1:2−1

(this is analogous to Equation (9)). The least squares estimate of α is

αˆ1:2 =

 αˆ1:2;a

αˆ1:2;b

 .

=

 3.80

3.35 

.

The corresponding estimates of φ1, φ2, σ2

1, and σ22 are 0.89, 0.96, 0.00059

and 0.00113, based upon which we calculate

Σˆ 1:2 .

=

 1.48 0.04

0.04 1.04 

.

The estimated parameters using data from both buoys are thus ˆα1:2;a .

=

3.80 ± 1.22 and ˆα1:2;b

.

= 3.35 ± 1.02. When we take into account sampling

variability, these estimates are consistent with what we got using the data

from each buoy by itself. The standard errors associated with these estimated parameters are smaller than those obtained when we used data from

just one buoy. The last two rows of Fig. 11 compare the observations, model

fits, and residuals at buoy 2 and buoy 1 when sources a and b and data from

the first complete wave at each buoy are considered.
20 Percival et al.

8. Summary and Discussion

We have described in some detail the SIFT inversion procedure by which

data from DART® buoys are used to estimate parameters that adjust precomputed models for the data arising from one or more unit sources. The

models are computed under the assumption that a magnitude 7.5 reverse

thrust earthquake occurred within each unit source. The parameters α essentially adjust the models to take into account earthquakes whose magnitude differs from 7.5, and does not coincide with, or is of greater spatial

extent than, a single unit source. We advocate a least square criterion for

estimating the parameters. This criterion has appeal for picking out appropriate estimates for the problem at hand and, in addition, leads to quick

computation and realistic assessment of the effects of sampling variability.

Version 2.0 of SIFT allows the user to select which unit sources are

to provide models for the observed data. For upcoming versions we plan to

implement ideas from statistical model selection to help the selection process.

Given the now completed U.S. array of DART® buoys in the Pacific and

Atlantic Oceans and given the extensive database of precomputed models

from over 1200 unit sources providing coverage of potential tsunami sources

in these oceans, a statistical approach to selecting appropriate combinations

of models and data can relieve some of the burden on users, particularly when

using SIFT during an ongoing tsunami event. Two simple approaches involve

step-up or step-down strategies. Suppose, for example, that we entertain

eight unit sources in all. A step-up approach would start with the unit

source whose model gives the best individual fit to the data (as measured

by residual variance) and would then use statistical tests to add models

from other unit sources one at a time, stopping when using an additional

model does not significantly improve the match with the data. A step-down

approach would start with a model that involves all eight unit sources and

then would use statistical tests to remove sources that do not appear to be

important contributors to the tsunami event.

Other areas of on-going research that might impact future versions of

SIFT are (1) the use of statistical criteria to select the amount of data to be

used to estimate the parameters; (2) allowing the precomputed models to be

shifted in time a bounded amount to compensate for uncertainties in propagation times between the unit sources and the buoys; (3) real-time detiding

of the DART® buoy data; and (4) use of a Kalman filtering formulation to

allow efficient updating of the estimated parameters as new data arrives.

Finally we note that an important topic that we have left unaddressed is

the connection between the estimated source parameters α and coastal inundation forecasts that use these parameters as inputs. The quality of these

forecasts depends on how well the source parameters can be used to successfully predict the tsunami event at locations far removed from the DART®

stations whose data formed the basis for the parameter estimates. For ex-
Extracting tsunami source parameters via inversion of DART® buoy data 21

ample, in the context of the Kuril Islands event focused on here, a further

test would be to see how well we can use the estimated α’s to predict what

was actually observed at far-field DART® stations along the West Coast

of the United States. An additional test could be based on tide gage data

at Crescent City and Hilo using nearshore standby inundation model (SIM)

forecasts (see Gica et al., 2009 for details). These additional tests need to

be done in order to fully validate the inversion procedure. In addition, work

is needed to determine how uncertainty in the estimated source parameters

translates into uncertainty in the SIM forecasts themselves.
22 Percival et al.

9. Acknowledgments

This work is funded by the Joint Institute for the Study of the Atmosphere

and Ocean (JISAO) under NOAA Cooperative Agreement No. NA17RJ1232

and is JISAO Contribution No. 1627. This work is also Contribution No. 3134

from NOAA/Pacific Marine Environmental Laboratory. The authors would

also like to thank R.L. Whitney for his help in putting this document together.

10. References

Draper, N.R., and H. Smith (1998): Applied Regression Analysis. 3rd edition, John

Wiley and Sons, New York, 736 pp.

Fletcher, R. (1987): Practical Methods of Optimization. 2nd edition, John Wiley

and Sons, Chichester and New York, 450 pp.

Gica, E., M. Spillane, D. Arcas, D.B. Percival, and V.V. Titov (2009): Development of NOAA’s Short-term Inundation Forecast for Tsunamis (SIFT) and a

case study of the 15 November 2006 Kuril tsunami. NOAA Tech. Memo. OAR

PMEL, NOAA/Pacific Marine Environmental Laboratory, Seattle, WA, to be

submitted.

Gica, E., M. Spillane, V.V. Titov, C.D. 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,

NOAA/Pacific Marine Environmental Laboratory, Seattle, WA, 89 pp.

Horrillo, J., W. Knight, and Z. Kowalik (2008): Kuril Islands tsunami of November

2006: 2. Impact at Crescent City by local enhancement. J. Geophys. Res., 113,

c01021, doi: 10.1029/2007JC004404.

Kowalik, Z., J. Horrillo, W. Knight, and T. Logan (2008): Kuril Islands tsunami of

November 2006: 1. Impact at Crescent City by distant scattering. J. Geophys.

Res., 113, c01020, doi: 10.1029/2007JC004402.