A wavefront orientation method for precise numerical determination of tsunami travel time

We present a highly accurate and computationally efficient method (herein, the “wavefront orientation method”) for determining the travel time of oceanic tsunamis. Based on Huygens’ Principle, the method uses an eight-point grid-point pattern and the most recent information on the orientation of the advancing wavefront to determine the time for a tsunami to travel to a specific oceanic location. The method is shown to provide improved accuracy and reduced anisotropy compared with the conventional multiple grid-point method presently in widespread use.


Introduction
Determining tsunami travel time is one of the fundamental roles of tsunami warning centres and other agencies tasked with estimating the possible impact of approaching tsunami waves on different coastal regions.In contrast to most other natural hazards, tsunami travel-time calculations can be performed well in advance.Because the travel time is the same when the source and the target are inverted, the inverse travel time for populated coastal sites have been computed and the corresponding databases of inverse travel-time maps created.When there is a tsunami event, the travel time to a specific coastal site is available immediately once information has been received about the initial parameters of the underwater earthquake and its location.Any initial travel-time estimation will be preliminary (it is not based on the actual extent of the tsunami source) and needs to be refined as new information regarding the tsunami source becomes available.
In addition to the need for accurate real-time travel-time estimates for tsunami alerts and warnings, highly reliable estimates of tsunami travel time (TTT) are required for optimizing the location of offshore tsunami warning stations used for refining information on the parameters of approaching tsunami waves (Poplavski et al., 1988).These stations require accurate advanced knowledge of the difference in wave travel time between specific coastal locations and the warning station (Titov et al., 2005).The greater the distance between the tsunami station and the coast, the greater the advanced forecast time for distant remotely generated tsunamis but the shorter the forecast time for locally generated tsunamis.The optimal warning system design should, therefore, be based on extensive determination of wave propagation times, which, in turn, requires an efficient and accurate method for calculating the tsunami travel time.Accurate estimation of observed wave arrival times is also needed for research delineating the tsunami source region.In particular, numerical modelling of the travel times for near-source locations can provide important information concerning the size and shape of the source region (e.g.Abe, 1973;Fine et al., 2005;Hayashi et al., 2012).Travel-time discrepancies of only a few minutes can lead to source-positioning errors of several tens of kilometres.Lastly, any discrepancy between the computed and observed tsunami travel times may indicate that the geometry of the source region is incorrect or that there may be physical effects, such as wave dispersion, nonlinearity, and coupling with elastic earthquake modes, which should be taken into account.These affects often require better wave travel-time accuracy than for tsunami warning purposes.
There are presently two primary approaches for calculating tsunami travel time from gridded bathymetry: (1) by kinematic wavefront propagation calculations based on Huygens' Principle (Shokin at al., 1987); and (2) by solving the dynamical equations of motion, typically using the finite-difference method (cf.Kowalik et al., 2005).The first method is the one most commonly used in modern tsunami travel-time calculations and is favoured by commercial software such as GEOWAVE (GEOWARE, 2013).Following Huygens' Principle, each of the points along a wavefront is a source for the tsunami waves.These points serve as start locations for travel-time computations to the next N neighbouring points.This set of points is referred to as "the pattern of neighbouring points".In GEOWAVE, N can have the values 8, 16, 32, 48 or 64, and the TTT between points is computed using the shallow water propagating speed formula, c = √ gh, where g is the gravitational acceleration and h is the water depth.Small meridional changes in gravity are incorporated in the methodology.Water depths are obtained from the gridded bathymetry provided by ETOPO (Amante and Eakins, 2009).If the arrival time to the nearest point is not determined, or if the previously computed arrival time is greater than the currently computed arrival time, the arrival time at the latter point is replaced with the newer value.At the end of the computational step, a new source point is specified as that frontal point that had the minimum travel (and arrival) times.This frontal point is converted to a permanent point to be used as source point for the next time step A critical step in the calculation is defining the spatial pattern of the N neighbouring points.The larger the value of N, the greater the isotropy of the wave field, and the better the accuracy of the timing calculation for open-ocean waves.On the other hand, a larger pattern also creates a broader region of uncertainty along the frontal zone, which, in addition to the computational overhead, creates problems in zones where the water depth is changing quickly.Moreover, a larger grid pattern can cause the tsunami arrival estimations along certain angles to "jump" over shallow water areas or over islands, leading to the need to integrate the travel time over the line connecting the two points, which significantly increases the computational complexity.In the case of a variable-depth ocean, this integration introduces an error arising from the fact that the actual path is no longer a straight line.
The wavefront orientation method proposed in this study uses a less extensive spatial grid pattern by taking advantage of the most recently gained information on the tsunami propagation in the neighbouring grid points.In particular, the method uses knowledge about the local direction of wavefront propagation obtained from previous calculations and, as a result, yields a more accurate travel-time calculation without the need to apply a large pattern of grid points.

The conventional method
The conventional method calculates tsunami travel time for each time step by dividing the spatial grid into three categories: (A) grid points for which the travel time has already been computed; (B) intermediate, or frontal, grid points where the travel time is to be evaluated; and (C) grid points that have yet to be reached by the advancing waves.At each time step, GEOWAVE starts the calculations from specific grid points in category A and then calculates arrival times for N neighbouring points in categories B and C (Fig. 1).New values replace the previously computed values if the new one is smaller.Figure 1 shows an example for a single step for the N = 32.In this example, the program computes the travel time for points 1 to 3 and points 20 to 32.Because points 4 to 19 belong to category A, they are not used in the current step calculation.
For any given grid pattern, there are errors in the conventional method due to the fact that not all propagation directions from the different source regions to the new grid point are covered by the pattern.Moreover, because of the structure of the rectangular pattern, the directional coverage is not uniformly distributed.As indicated by Fig. 1, the biggest gaps are in the directions along the grid axes, corresponding to those directions in line with the grid points (i.e. between the direction from centre O to point B1 and to point B2, or between directions from point O to point B24 and to point B25, and so on); consequently, the main arrival errors are in these directional sectors.Note that we are restricting our analysis to the case for uniformly propagating wave velocity, where the relative errors in travel time are equivalent to the relative path length errors.Directional errors for the conventional method can be estimated analytically.Each pattern has a fixed number (or list) of directions of propagation, α i , (i = 1, 2, . ..), corresponding to the grid pattern number N. If the directional angle β from the source to a given point is not in the list, the fastest way to reach a given point in the model is along one of the listed directions, α k , and then in the adjacent direction α k+1 such that α k < β < α k+1 . (1) The resulting accumulated distance for the two-direction pass will be greater than the distance between the source O and the given point A (Fig. 2).The relative increase in path length, δ (corresponding to an increase in the tsunami travel time), is The relative change in path length, δ, is a maximum when the direction of propagation, β, is equal to the mean of the bounding list directions, α k and α k+1 ; specifically, whereby The value of δ max is highest when either α k or α k+1 coincides with one of the coordinate axes.
The angle β max and statistics for the relative error err = δ− 1 are listed in Table 1.Results for the conventional method are listed as P8, P16, P32, P48 and P64 according to the number of points N used in the travel-time computations.Errors for the conventional method are always positive or zero, and only zero in those directions lying along the direct interconnection between spatial grid points.Based on Table 1, the directional errors are typically small and decrease with an increase in the number of points used in the computational grid pattern.However, because the relative error remains uniform with distance, the absolute error can become significant over long distances, such as in the case of trans-Pacific tsunami propagation.

The wavefront orientation method
We propose a new methodology for calculating tsunami travel time which, in addition to having a small spatial grid pattern, makes use of the most recently derived information regarding the tsunami travel time at the point closest to the current source position.Let O be the current source point and let A5 be the nearest point for which the travel time has been determined (Fig. 3).Because the time it has taken the tsunami to reach the current source position is, by definition, a maximum, the difference in travel (or arrival) time Table 2. Relative errors for the proposed wavefront orientation method and a plane square grid.Distance from the source location at point O (0, 0) is measured in terms of the specific point number on the grid.The larger the grid number, the greater the distance from the source region.RMS is the root mean square of the error, and STD is the standard deviation of the errors.
If we assume that the tsunami wavefront is a straight line in the vicinity of the source point, and that the tsunami arrival time, t, in the triangular region O-A5-B7 in Fig. 3 is a linear function of the local spatial coordinates x and y, then where c −1 is the wave slowness (the inverse of the wave speed, c) and the orientation angle, φ, of the tsunami wavefront can be found using known values of the time difference, t OA5 , corresponding to the arrival time between points O and A5.Specifically, where x is the incremental grid step in the x direction.The travel time to point B7 from point O will be where y is the grid step along the y direction.Equation ( 7) allows for more accurate computation of the travel time to point B7 than the conventional method.However, because it is based on previously computed travel times at two points, our method requires some initialization, which can be provided by the conventional method using the simple 8-point algorithm.Unlike the conventional method, this new approach has no fixed directional errors.As a consequence, relative errors decrease with distance from the start location (origin) of the tsunami event.This provides an improvement over the conventional method, whose level of error depends on the initialization pattern.
Table 3. Errors in tsunami travel time for a constant (4000 m)-depth ocean with bathymetry gridded at 1-arcminute steps for both the conventional and proposed methods.Values are compared to exact values obtained from an analytical solution.Results for the conventional method are denoted with a "P"; those for the waveform method with an "F".RMS is the root mean square of the error, and STD is the standard deviation of the errors.3 Method comparison

Constant-depth ocean on a spherical geographical grid
One of the primary concerns in tsunami research is determining the travel times for trans-oceanic (mostly trans-Pacific) tsunamis.It is, therefore, of interest to determine how our newly proposed method performs relative the conventional method for a spherical, constant-arc step grid for the Pacific Ocean.We first examine the method's performance for an ocean of uniform depth and then compare the results to those for an exact analytical solution.Table 3 shows the absolute errors in tsunami travel time obtained for an ocean of 4000 m depth for both methods.
According to Table 3, the maximum timing error decreases with grid pattern number but remains measurable even for the conventional P64 algorithm.The errors are significant for the P32 and, especially, the P16 algorithms.In addition, the conventional method has a positive bias, which needs to be corrected.The GEOWARE software decreases this bias by introducing a "correction coefficient".In contrast, the proposed wavefront orientation method is almost free from errors; the maximum differences between the numerical results and the theoretical (analytical solution) values do not exceed 12 s over the entire propagation distance of 15 000 km.
Figure 4 shows the distribution of errors for the different methods.As for the plane wave case discussed in Sect.2, the accumulated timing errors for the conventional method generally increase with distance, although the distribution of errors differs from that for the plane wave case.In particular, the errors concentrate near both sides of the latitude meridians passing through the source region but are almost absent near the parallels of longitude passing through the source region.This is because, for high latitudes, the wave propagation beams are not straight lines in the spherical coordinate system (except in the meridional direction).Specifically, the beams cannot be along parallels but must be along great circles.In the near-meridional direction, the relative error is similar to that for the plane wave case -i.e. about 0.5 %, 1.3 % and 2.7 % for the conventional grid cases P64, P32 and P16.For the new wavefront model, the error distribution has little directional bias mainly because it is an almost error-free calculation to begin with.

Pacific Ocean with realistic seafloor bathymetry
For estimations of tsunami travel time under realistic ocean conditions, accurate measurements of water depth are of major importance.In addition to its direct effect on the actual travel time, changes in water depth can affect the direction error as the wave beams undergo convergence and divergence.As shown by Satake (1988), these effects completely change the beam pattern of tsunami energy flux in the Pacific Ocean compared to that for an ocean of uniform depth.The difficulty with estimating the reliability of the tsunami traveltime calculation over a variable-depth ocean is that it cannot be compared to an analytical solution.Moreover, because of measurement uncertainties, the observed travel times to tsunami recording sites such as tide gauges and bottom pressure recorders are not sufficiently accurate to judge the numerically derived results.We are, therefore, limited to comparing differences between the different numerical methods, taking into account the hierarchy of the conventional method, denoted as models P8 to P64, for which the accuracy increases with increasing spatial pattern number.
It is important to understand how the conventional and wavefront methods compute the travel distance between two points over a variable-depth ocean.As noted earlier, GEOWARE computes the travel times over straight lines connecting pairs of points.For each line, the travel time is computed as the sum of the travel times for each individual grid cell that a particular line crosses.Inside each grid cell, the algorithm interpolates the inverse celerity so that the computed time between neighbouring points corresponds to the travel time based on the inverse celerity for individual grid cell regions of uniform depth.Thus, the uniform depth value used is always closest to that for the shallower point.
The wavefront method uses the same type of interpolation.However, because of the short spatial template (distance between neighbouring points), each inner-cell computation is limited to just one cell.We have run computations for the conventional and wavefront methods for a realistic ocean using the identical source position as for the uniform depth ocean (i.e.175 • W, 50 • N) and covering the same spatial domain (lower left corner: 95 • E, −72 • S; upper right corner: 80 • W, 65 • N).We have also used the same ETOPO-1 dataset for both the conventional and wavefront method.To avoid complications arising from different interpolation schemes in coastal areas, we have restricted our comparisons to openocean points only, for which water depths are greater than 3000 m. Results of these comparisons are presented in Table 4 and Fig. 5.
Because the conventional method always overestimates the tsunami travel time regardless of the number of grid points used, the travel-time errors using conventional model P16 versus using model P64 (written P16-P64) or model P32 versus P64 (P32-P64) are on the low side relative to the actual analytically derived timing errors.For this reason, the errors for conventional model P16 versus our wavefront Table 4. Statistical parameters for differences in computed trans-Pacific tsunami travel times using the ETOPO-1 dataset for the conventional and the proposed wavefront orientation methods.F8 denotes runs using the new waveform orientation method based on the 8-point GEOWAVE algorithm.RMS is the root mean square of the error, and STD is the standard deviation of the errors.3.

Methods Minimum difference
As shown by the bottom 3 rows of Table 4, the difference in travel time P16-F8 is close to the difference P16-P64, while the differences for P32 versus F8 are close to the P32 versus P64 differences.The last row in Table 4 shows differences between the best conventional method results and the wavefront orientation method results.These values are clearly quite low, with maximum and minimum differences only 1.6 min and −1.7 min, respectively.On the other hand, the differences of P16-F8, P32-F8 and P64-F8 are consistent with the first 3 rows of Table 3 (which compares P16, P32 and P64 with the exact solutions for a uniform depth ocean).This indicates that F8 (which combines the P8 algorithm with the wavefront orientation algorithm) outperforms all of the conventional runs, including P64, even for the variable-depth case.
Figure 5 shows the difference distributions presented in the same order as in Table 4. Generally, the distributions for the cases P16-F8, P32-F8 and P64-F8 (panels a1, b1 and c1) are similar to the error distributions shown in Fig. 4 (panels a, b, and c, which represent the differences between the conventional method and the exact solution for the uniform depth ocean case).This again confirms that the new method provides more accurate time-of-arrival values than the conventional method.

Discussion and conclusion
The previous results demonstrate that our proposed "wavefront orientation method" provides a more effective and accurate methodology for calculating the travel times for transoceanic tsunamis compared to the conventional method.Moreover, our method works especially well for large bathymetric grid regions; the larger the gridded array, the more advantageous the method.
It is important to note that the accuracy of our newly proposed method should not be confused with the general accuracy of tsunami travel-time calculations for specific events.Actual tsunami travel times depend on many factors that are independent of the particular algorithm being used (Wessel, 2009).With respect to the two methods discussed in this study, both assume small depth changes inside each grid cell.In practical terms, the actual bathymetry often varies more rapidly than is assumed by the numerical methods.This raises the question of whether any method is accurately representative of real situations.If the water depth (or, more to the point, the wave celerity) varies between neighbouring points, the travel-time error can be of the order of the computed values themselves (which in turn strongly depends on the type of interpolation used between points).In the deep ocean, the travel time between points can be small, and, accordingly, the error will be small.Unfortunately, the biggest change in water depth usually occurs in shallow areas, where the travel time between points is long and, accordingly, the errors are especially large.The only way to avoid such errors is to increase the grid resolution Actual travel-time accuracy depends on the accuracy of the bathymetric dataset.In the open ocean, present gridded global datasets are sufficient for most practical applications.However, this is not the case for coastal areas.Moreover, local, high-resolution bathymetric data are not fully incorporated in the ETOPO datasets provided with the GEOWARE software.This can lead to significant errors in the traveltime estimates.In addition, global datasets do not resolve narrow passes in coastal areas, necessitating the use of local nested bathymetric data.The limitations of the bathymetric data are compounded by the fact that the tsunami generation region is typically poorly defined.In a scientific twist, observed tsunami travel times are often needed to better define the source region that generated the waves in the first place.
There are other physical factors which can cause the actual speeds of tsunami waves to depart from wave speed estimates obtained from simple linear shallow water theory.One of the most important factors is frequency dispersion related to nonhydrostatic effects.Non-hydrostatic effects are especially important during far-field wave propagation than during 1 Figure 5. Maps of differences in tsunami travel time (minutes) using numerical solutions for 2 the ETOPO-1 grid for the Pacific Ocean.Left-hand panels: inner-method differences in travel 3 time for the conventional solutions only; (a1) between P16 and P64; and (b1) between P32 4 and P64.Right-hand panels: between model differences in travel time derived using the 5 conventional and waveform orientation methods; (a2) between conventional model P16 and 6 the waveform method; (b2) between conventional model P32 and the waveform method; and 7 (c2) between conventional model P64 and the waveform method.8 Fig. 5. Maps of differences in tsunami travel time (minutes) using numerical solutions for the ETOPO-1 grid for the Pacific Ocean.Lefthand panels: inner-method differences in travel time for the conventional solutions only; (a1) between P16 and P64; and (b1) between P32 and P64.Right-hand panels: between model differences in travel time derived using the conventional and waveform orientation methods; (a2) between conventional model P16 and the waveform method; (b2) between conventional model P32 and the waveform method; and (c2) between conventional model P64 and the waveform method.
near-field tsunami wave propagation, or when the tsunami source area is relatively small (shorter waves) and/or located in relatively deep water (González and Kulikov, 1993).Although non-hydrostatic effects modify tsunami wave propagation in a depth-variable ocean, the travel times for nonhydrostatic waves can be computed in the same way as for shallow water theory by examining each wave frequency separately.This is because, for a given frequency, the arrival time is determined by the group velocity of the dispersive waves, which depends on the depth only.
In shallow water regions, frequency dispersion becomes less important.However, non-linear effects arising from wave interaction with the tsunami wave field, and from interaction between the tsunami waves and tidal currents, can modify the travel times.In estuarine areas, tidal motions and river flow may also alter the travel times.
As shown by Watada (2013), ocean stratification, seawater compressibility and elasticity of the solid earth also affect the gravity waves.All these factors can generate noticeable changes to tsunami travel times.
Regardless of the factors affecting tsunami wave propagation, it is clear that accurate estimation of tsunami travel time is of considerable importance.Because Huygens' Principle does not depend on wave speed formulations, the above factors (excluding perhaps non-linear effects) can be included in our travel-time algorithm.Thus, in addition to providing improved computational efficiency and better accuracy compared to the conventional methodology, our proposed method can be used to enhance general tsunami research.

Fig. 1 .
Fig.1.Schematic showing one time step in the calculation of tsunami travel time for an N = 32 spatial pattern of neighbouring grid points for the conventional method.Grid points are divided into three categories: (A) points for which the travel time has been calculated (dark grey area); (B) points in the frontal zone (light grey area) for which the time is being calculated; and (C) points where no calculation has yet been made (unfilled area).

Figure 2 .Fig. 2 .
Figure 2. Schematic showing the source of the directional error arising from a calculation based on wave propagation along a path with two separate line segments (solid lines).The dashed line denotes the path for which there is no directional error but for which there is no intermediate point, B.

Figure 3 .Fig. 3 .
Figure 3. Schematic showing the N = 8 pattern of neighbouring spatial points 4 proposed wavefront orientation method.The shading is the same as for the c 5 method presented in Fig. 1.The group of three small arrows show the direction of 6 propagation.7 8

Figure 4 .Fig. 4 .
Figure 4. Spatial distributions of tsunami travel time errors (minutes) for a uniform depth 5 (4000 m) ocean using: (a, b, c) the conventional method for N = 16, 32, and 64 neighbouring 6 point patterns, respectively; and (d) the wavefront orientation method based on an N = 8 7 neighbouring point algorithm.8 9

Table 1 .
Directional error parameters for the conventional method (denoted by the letter "P") applied to a plane square grid.RMS is the root mean square of the error, and STD is the standard deviation of the errors.Column 2 gives the angle yielding the maximum error.