Tsunami Methodology

Methodology

    In order to model tsunami propagation in the Pacific Ocean at both the scale of the ocean basin and at the local scale of coastal inundation, I needed to use a series of nested grids. Kochimura et. al. note that the appropriate temporal resolution of a spatial model is dependent on the grid spacing, through the relationship delta(x)² = 4h + (gh)delta(t)², where delta x is the grid spacing and delta t is the appropriate time step. In a sensitivity analysis of tsunami models at three different spatial scales, Myers and Baptista (2001) found that even a grid with 10 meter spacing can significantly alter the accuracy of bathymetry data and the location of features such as river deltas. Such horizontal and vertical uncertainty can lead to discrepancies as large as ± 6 meters in runup estimates.
    I structured my project after Myers and Baptista (2001) using three different grid resolutions for my model. However, the resolution of my data sets was no where near 10 meters: the finest resolution was 0.0007 degrees latitude, which is ~77 meters. The Smith and Sandwell data provides global coverage between 72 degrees north and south: I used a window from 95 East to 295 East longitude and from 72 North latitude to 72 South latitude. The data is stored as a Mercator projection: I used E.R. Mapper to convert it to geodetic coordinates, based on a spherical datum. To increase the resolution of the Smith & Sandwell data to the resolution of the BC DEM (1000 meter (30 arc second) resolution), I converted the raster data to vector points, treating it as point data for interpolation. Unfortunately, IDRISI was unable to handle this interpolation, so instead I mapped the data as raster points onto the blank (below sea surface) portion of the BC DEM, and then used a 7x7 minimum filter to interpolate the data, with fairly accurate results. I took a random sample of ten overlapping points and ran a regression analysis, which resulted in the relationship y = 1.006x, r²= 0.9961.
    IDRISI was likewise unable to handle interpolating the maximum waveheights, so I instead interpolated contours from the points in ArcView, then imported the contours to IDRISI, where I created a TIN and a raster surface from the contours (using points alone resulted in very rigid lines).
    To convert the hillshade JPEGs, I first separated them into colour bands 1, 2 and 3 using the separate module, then subtracted band 1 from band 3, and divided by band 2. This resulted in a relatively scaled DEM, with values near ± 1. I created a window using the southern Vancouver Island DEM, and tested only the above sea level portion to see if it could be linearly scaled. I wasn't able to run a complete spatial regression in IDRISI because the data were at different resolutions, but using the display scales (over exactly the same spatial extent), I arrived at the relationship y = 581.22x, r² = 1.
    To model tsunami waves, I used a simple linear periodic function:

Wave elevation (n) = A*sin(kx-ot), where A is the amplitude, k is the wavenumber (k = 2*π/T*C,where T is the wave period and C is the celerity, or wave speed), o (denoting sigma) equals 2*π/T, x is the displacement, and t is time (Rahman, 1994).

    In linear shallow water wave theory, waves are assumed to be non-dispersive, meaning that their phase velocity equals their group velocity, and their wavelengths are constant over time. This holds well for tsunamis in the deep ocean, as their wavelengths are much longer than the depth of the ocean, causing them to act like shallow water waves. Ideal shallow water waves move at a speed that is independent of wavelength, depending only on the depth: C = (h*g)^0.5, where h is the water depth, g is the acceleration of gravity (9.8065 ms^-2). However, as the waves move in to shore, and start to interact with irregular bathymetry, linear theory starts to break down, and the Boussinesq, or cnoidal wavetheory becomes more appropriate (Rahman, 1994). I did not, however, succeed in my attempts at applying non-linear wave equations.
    To compute distances around objects such as islands, which the spherical distance function won't do, I calculated both a spherical distance and a cost grow distance, using a land mask (ocean = 1; land = 0) as an absolute barrier. I then ran a regression module to determine how to scale the cost distance to real units. Depending on the distance involved and the amount of contortioning around islands, the regression usually fell between r²= 0.95-0.99.
    To model reflection and refraction, I modified the amplitude and wavelength following Rahman (1994), where the ratio of new speed or wavelength to the original speed or wavelength = 2*π*(h/g*T²)^1/2, and the amplitude ratio is (((1-sin²a_o*4*π²*(h/g*T²))/cos²a_o)^-1/4*(16*π²*h/(g*T²))^-1/2, where a_o is the angle of incidence between the wavefront and the bathymetric contours, which I approximated using the aspect, and force direction of the wave.
    For runup calculations, I used R = πA(2L/λ₀), where A is the amplitude of the last wave peak offshore, L is the distance from that point to the shore, and λ₀ is the initial wavelength (Choi et. al., 2002). To account for wave attenuation with distance, I used the formula H = H_o*(R_o/R)^1/2 from Choi et. al. (2002). Figure 1 gives an abbreviated conceptual model for wave design.
    For the earthquake parameters, there are several different methods for estimating initial wave heights from earthquake magnitude. Pelinovsky et. al. (2002) suggest log η_o = −4.31 − 4.36 log h_f+ 1.45M_o, where η_o is the initial water displacement, h_fis the depth of the earthquake focus, and M_ois the magnitude of the earthquake. They also give a relationshiop for the length of the earthquake fault: log l_o = 0.5M_o − 1.8, which leads to a hypothetical ellipsoidal source area with semi-major axis a = l_o + 2h_f, and semi-minor axis b= h_f, where lo is the length of the fault area, and h_f is the depth of the focus. Tsunami magnitude (m) can be related to earthquake magnitude through the equation m = (2.61 ± 0.22)*M - (18.44 ± 0.52), which is in turn related to maximum wave height m = log₂n_max. Tsunami energy is generally 1/10 the seismic energy (Murty, 1977).

To keep my comparison of the Alaska tsunami and the hypothetical Cascadia tsunami simple, I focused on the difference in angle of incidence and the resulting difference in the angle of reflection, using the resultant module in IDRISI to combine the incident wave speed and direction with the bathymetric aspect to determine the resulting direction of wave motion. I also calculated progressive wavefronts for both cases using the cost grow module using a friction surface based on the linear shallow water wave speed (C = (h*g)^0.5).