A S S E S S M E N T OF T S U N A M I H A Z A R D S O N T H E BRITISH COLUMBIA COAST DUE TO A LOCAL MEGATHRUST SUBDUCTION E A R T H Q U A K E By Max Kin-Fat Ng B. Sc. (Physics and Mathematics) University of Victoria A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE m THE FACULTY OF GRADUATE STUDIES PHYSICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA January 1990 © Max Kin-Fat Ng, 1990 In presenting degree freely at this the available copying of department publication of in partial fulfilment University of British Columbia, for this or thesis reference thesis by for his this and scholarly or thesis for her of P H Y S I C S The University of British Columbia Vancouver, Canada Date DE-6 (2/88) A p r i l , 1990 I I further purposes gain the shall requirements agree that agree may representatives. financial permission. Department study. of be It not that the be Library an by understood allowed the advanced shall permission for granted is for make extensive head that without it of copying my my or written Abstract Strong evidence suggests that the Cascadia subduction zone, off the west coast of Canada and the United States, is strongly seismically-coupled and that a possible megathrust earthquake might occur in that area in the near future. A study of tsunami hazards along the Canadian west coast associated with such a hypothetical earthquake is presented in this report. Numerical simulations of tsunami generation and propagation have been carried out using three models based on shallow water wave theory. Three cases of ground motion representing the ruptures of different crustal segments in the area have been examined. Computed results provide information on tsunami arrival times and a general view of the wave height distribution. The outer coast of Vancouver Island was found to be the most strongly affected area. At the head of Alberni Inlet, wave amplitudes reached up to three times the source magnitude. Inside the Strait of Georgia, the wave heights are significant enough to receive closer attention, especially in low-lying areas. n Table of Contents Abstract ii List of Figures v Acknowledgement xiv 1 Introduction 1 2 Earthquake potential in the Pacific Northwest 4 3 Tsunami Generation and Propagation Models 12 3.1 3.2 The deep-ocean model 12 3.1.1 Tsunami generation and propagation 12 3.1.2 The geographical location of the deep-ocean model 15 3.1.3 The coarse grid 15 3.1.4 Numerical solution of the deep-ocean model 20 3.1.5 Boundary conditions 21 3.1.6 Model testing and reliability 22 Thefinegrid model 34 3.2.1 The governing equations 34 3.2.2 Numerical scheme 35 3.2.3 Boundary Conditions 36 3.2.4 Junctions between the coarse andfinegrids 36 in 3.3 4 5 Tsunami propagation in inlets 51 3.3.1 The inlet model 52 3.3.2 Numerical solutions 52 3.3.3 Boundary conditions 53 Results 57 4.1 Rupture of the Juan de Fuca, Winona and Explorer segments 58 4.2 Rupture of the Juan de Fuca segment alone 62 4.3 Rupture of the Winona and Explorer segments 63 4.4 Sources of different magnitudes: a test of model linearity 64 Conclusion 110 Bibliography 114 Appendices 116 A Approximation of a spherical polar grid by a rectangular grid 116 B 123 Constricted flow in Burrard Inlet iv List of Figures 2.1 The Cascadia subduction zone, showing the location of four different segments. The rupture or 'deformation' front is denoted by dashed lines. (After Rogers, 1988) 2.2 7 The three-stage deformation cycle for a subduction zone. The continental plate is shown on the right while the oceanic plate is on the left. (After Dragert and Rogers, 1988) 2.3 8 (a) Contours of sea-bottom displacement (in cm) with all three segments ruptured as a single unit: the Explorer segment is shown in dashed lines with the Winona segment in the north and the Juan de Fuca segment in the south. Ground motion includes a maximum uplift of 5m along the rupture front and a subsidence of about 0.5 m in coastal regions. . 9 2.3 (b) Two-dimensional cross-section view of the bottom displacement. . . 10 2.3 (c) Three-dimensional view of the bottom displacement 3.1 Location map of the study area for the west coast of British Columbia. 11 The 5-km coarse grid (A) is used for wave propagation modelling in the ocean-continental shelf region, while the 2-km fine grid (B) is for relatively shallower waters. (After Dunbar et al, 1989) v 18 3.2 The coarse grid used in the deep-ocean model. The dots in the grid layout represent vertices of grid elements. The central vertical line denotes the open boundary before extension. A flat bottom ocean of depth 2,700 m is assumed for the extended part. The shaded area represents a test tsunami source of upward displacement of 2 m. Marked stations are chosen in the order of increasing angles of wave incidence onto the open boundary: waves incident perpendicularly on the boundary for stations starting at grid point (1,75); for stations starting at (4,113) and continuing northwards, the incident angles become more oblique. Corresponding elevation time series are shown in Fig.3.4 and 3.7. ... 3.3 Arakawa-C-type grid. (After Henry, 1982) 3.4 A comparison of computed elevations between the original grid (solid 19 25 line) and the extended grid (dashed line) as a test of spurious reflections off the open boundary. Locations of stations are indicated in Fig.3.2. In (a), as the waves are incident normally, the leading waves almost coincide for the two cases, and sucessive waves also agree quite well. As the incident angles increase, as in (b), small differences are observed both for the leading and subsequent waves. In (c), as the angles deviate more from normal incidence, the entire wave train is modified by the spurious reflections. The vertical scale is in metres 26 3.4 Continued 27 3.4 Continued . 28 3.4 Continued 29 3.5 Finite-difference forms of the 2-d radiation boundary condition on the three seaward boundaries 30 vi 3.6 A comparison of computed elevations for three different boundary cases: original grid with 1-d radiation boundary conditions (solid line), original grid with 2-d conditions (dashed line), extended grid with 1-d conditions (dotted line). Clearly, the use of the 2-d boundary condition approximates the extended grid results better. The Coriolis terms are dropped in computations when applying the 2-d conditions. Stations are located in the northern half of the grid, as shown in Fig.3.2. The vertical scale is in metres 3.7 31 The percentage difference (vertical axis) between calculated elevations using linearized and non-linearized equations. Station locations are indicated in Fig.3.2 32 3.7 Continued 33 3.8 Geographical location of the fine grid model. (After Crean et al, 1988). 39 3.9 Computational array of the fine grid model. (After Crean et al, 1988). . 40 3.10 Numerical scheme used in the fine grid model. (After Crean et al, 1988). 41 3.11 Coarse grid layout at the Juan de Fuca Strait entrance, (a) Location of the junction between the coarse and fine grids, (b) Two configurations of open boundaries are tested at the Strait entrance of the coarse grid. 42 3.12 A comparison of sea levels calculated from implementing the open boundary through the u-v points (dashed line) and through only the u-points (dotted line). The difference is perceivable only for some brief moments at point (67,18). Results of the original solid boundary are shown in solid line. In panel (d), only the results through u-points are shown. Locations of stations are indicated in Fig.3.11a. The vertical scale is in metres 43 vn 3.12 Continued 44 3.13 The extended coarse grid layout. Solid lines denote the original boundaries(Fig.3.2), whereas the dashed lines indicate the fine grid boundaries. In order to have a sufficient number of depth points (and thus water surface elevation points) covering the junction and also to reduce possible spurious reflections from the open boundaries, the coarse grid was extended in two entrances: Juan de Fuca Strait and Johnstone Strait. Because the width of Johnstone Strait could not be resolved within the coarse grid dimensions, the extension made here is not exact. The approximation is however justifiable since the junction consists of one grid point only . Johnstone Strait is assumed to have a depth of 180m 45 3.14 The Juan de Fuca Strait junction between the coarse andfinegrid models. Juan de Fuca Strait is extended as far up the Strait head as possible in the coarse grid model. For 129 < I < 133, depth points are extrapolated from thefinegrid data as obtained from Crean et al (1988). For 134 < I < 143, all points are given a constant depth of 150 m. Such an approximation is acceptable as this part of the region will be repeated in thefinegrid calculation 46 3.15 A comparison of computed elevations (in m) in the Juan de Fuca Strait area. Solid and dashedfinesdenote respectively the results before and after the grid extension. Station locations are indicated in Fig.3.14. . . 47 3.15 Continued 48 3.16 Junction of thefinegrid and extended coarse grid at the Juan de Fuca Strait . 49 vm 3.17 Four-point interpolation scheme used in calculating thefinegrid open boundary values 50 3.18 The geographical location of Alberni Inlet. (After Marshal Macklin Monaghan, 1986) 55 3.19 Numerical scheme used in the inlet model 4.1 56 Tsunami travel times (in minutes) in the deep-ocean model domain, (a) All three segments rupture together, (b) Only the Juan de Fuca segment ruptures. Contour intervals are 10 minutes 66 4.1 Continued 67 4.2 The coarse grid layout showing the locations at which time-histories of water level are plotted. Station locations are mostly presented in geographical order, beginning at Port San Juan and going north 4.3 68 Computed elevations (in m) of one of the source points located above the deformation front for the case of all three segments ruptured together. The initial elevation is 5.0m. Station location is shown on Fig.4.2. . . . 69 4.4 Time-series of computed elevations in the deep-ocean model for the case of all three segments ruptured together. Station locations are indicated in Fig.4.2. The vertical scale is given in metres 70 4.4 Continued 71 4.4 Continued 72 4.4 Continued 73 ix 4.5 Spatial attenuation of wave energy. Leading wave amplitudes, given in centimetres at grid points located at the lower-left hand corners of the numerals, are seen to decrease with distance from the source. The deformation front is indicated with an uplift curve of 500 cm. Elevations at coastal points may be amplified and are not included here 4.6 74 The fine grid layout showing the locations at which time-histories of water level are plotted. Station locations are presented in geographical order, beginning at Juan de Fuca Strait and going into the Strait of Georgia-Puget Sound system 4.7 75 Time-series of computed elevations in thefinegrid model for the case of all three segments ruptured together. Station locations are indicated in Fig.4.6. The vertical scale is given in metres 76 4.7 Continued 77 4.7 Continued 78 4.7 Continued 79 4.8 A comparison of thefinegrid model results with and without the nonlinear advective terms, denoted respectively by solid and dashed lines. Station locations are indicated in Fig.4.6. The vertical scale is given in metres 80 4.8 Continued 81 4.9 Computed elevations at Alberni Inlet for the case of all three segments ruptured together. The model is driven at Barkley Sound. The vertical 4.9 scale is given in metres, a) Barkley Sound sea level displacements. . . . 82 b) Inlet head sea level displacements 83 x 4.10 Contours of sea-bottom displacement for the case when the Juan de Fuca segment is ruptured alone. Contour intervals are given in centimetres 84 4.11 Time-series of computed elevations in the deep-ocean model for the case of the rupture of the Juan de Fuca segment alone. Station locations are indicated in Fig.4.2. The vertical scale is given in metres 85 4.11 Continued 86 4.11 Continued 87 4.11 Continued 88 4.12 A comparison of the junction points at the Juan de Fuca Strait for the first two cases of source displacement: the rupture of all three segments (solid line) and of only the Juan de Fuca segment (dashed line). Station locations are indicated in Fig.4.2 89 4.12 Continued 90 4.13 Time-series of computed elevations in the deep-ocean model for the case of the Winona and Explorer segments ruptured together. Station locations are indicated in Fig.4.2. The vertical scale is given in metres. 91 4.13 Continued 92 4.13 Continued 93 4.13 Continued 94 4.14 Time-series of computed elevations in the fine grid model for the case of the Winona and Explorer segments ruptured together. Station locations are indicated in Fig.4.6. The vertical scale is given in metres 95 4.14 Continued. 96 4.14 Continued 97 xi 4.14 Continued 98 4.15 Computed elevations at Alberni Inlet for the case when the Winona and Explorer segments are ruptured together. The model is driven at Barkley Sound. The vertical scale is given in metres, a) Barkley Sound sea level displacements 99 4.15 b) Inlet head sea level displacements 100 4.16 Three-dimensional views of water surface elevations at various time intervals. The shoreline is represented as a 20-cm step. Maximum wave heights are indicated in metres at each plot 101 4.16 Continued 102 4.16 Continued 103 4.16 Continued 104 4.16 Continued 105 4.17 The near linearity of the wave equations is revealed by applying two sources with one's magnitude being double of the other's. Solid lines denote the maximum elevation of 5 m whereas dashed lines correspond to the 2.5 m source. All three segments are ruptured in both cases. Station locations are indicated in Fig.4.2 4.17 Continued 106 107 4.18 The near linearity of the wave equations is shown by comparing results of two sources: first for the case when all three segments are ruptured as a whole (solid lines), and in the other case by adding the calculations for each segment rupturing separately 108 4.18 Continued 109 A.l A parametric surface S 117 xii A.2 Earth represented as a spherical surface A.3 A spherical segment S of sides pAc/> and p sin </>A# 118 118 A. 4 A spherical surface approximated by a rectangular grid. The numerals indicate the latitude-longitude range of the coarse grid used in the deep- B. l ocean model 119 Burrard Inlet showing the three basins. (After Baines (1958)) 124 xm Acknowledgement I wish to thank Professor Paul LeBlond for his supervision and encouragement within these two years, and for his helpful suggestion and patience in reviewing various drafts of the manuscript. My thanks also go to Dr. Tad Murty and Dr. Faulkner Henry of the Institute of Ocean Sciences for their helpful suggestion and assistance, to Dr. Don Dunbar of Seaconsult Marine Research Ltd. for providing technical information on the shelf and inlet models, and to Dokyo Lee and Wing Quan of IOS for assistance at the early stage of applications of the shelf and GF7 models. Finally, I express my special thanks to Dr. Toshiyuki Hibiya for his stimulating discussion and useful suggestion. xiv Chapter 1 Introduction The atmosphere, the earth and the ocean constitute the essential elements which man depends upon to survive, yet they also give rise to natural disasters that claim thousands of lives. Among these hazards are earthquakes and tsunamis. The effects of earthquakes are often catastrophic. Not only do seismic waves directly cause devastating damage but associated effects such as tsunamis can sometimes be triggered and pose additional threats to life and property. These seismic sea waves, which are often incorrectly termed tidal waves, can propagate long distances from their sources and may be amplified in bays and harbours, causing severe flooding. Tsunamis of seismic origin with heights of up to 30 metres have been reported (Abe, 1979). From 1947 to 1981 alone, the global loss of life due to tsunamis has reached 8,568 persons, averaging 252 deaths per year (Thompson, 1982). Tsunamis occur frequently in the Pacific Ocean because of the large number of earthquake zones on its shores. Coastal regions on both sides of the ocean are susceptible to inundations by these sea waves originating in the Pacific basin. The west coast of Canada is no exception. Historical earthquakes such as the 1960 Chilean and the 1964 Alaskan earthquakes generated tsunamis that reached Canada's west coast. Studies of these tsunamis have been carried out by several authors: Loucks (1962), Murty and Boilard (1970) and more recently Dunbar et al (1989). In particular, the tsunami wave response along the British Columbia coast due to distant 1 Chapter 1. Introduction 2 sources was studied in detail in the latter paper. An assessment of tsunami hazards in this area associated with a local event has, however, not been treated. As wave characteristics change appreciably through propagation, wave signatures on the British Columbia coast might be quite different for the cases of distant and near-field sources. This paper presents an estimate of the behaviour of a tsunami generated locally by a hypothetical megathrust earthquake. A series of great earthquakes, of magnitudes of about 8, or a massive earthquake of magnitude over 9, has been predicted to take place in the Cascadia subduction zone of the Pacific Northwest (Heaton and Kanamori, 1984; Heaton and Hartzell, 1987; Rogers, 1988; Dragert and Rogers, 1988; Dragert and Horner, 1989). At the edge of the continental margin, sections of the ocean floor are slowly converging onto the North American Plate and descending beneath it (that is, subducting). Analysis of sea floor magnetic anomalies has suggested that this subduction started 40 million years ago (Rogers, 1988) and is continuing presently at a rate of 3 to 4 cm per year. One of the Cascadia subduction zone's most distinguishing features is the youthfulness of the oceanic lithosphere being subducted, varying from 6 to 10 million years. There are other areas in the world with tectonic features similar to the Cascadia subduction zone: the subduction zones in Southern Chile, southwestern Japan and Colombia (Heaton and Hartzell, 1987; Rogers, 1988). All of these zones have experienced large thrust earthquakes in historic times, of magnitudes measuring eight or above on the Richter scale. Large local tsunamis were generated in all these events. The maximum run-up heights of the tsunamis locally generated by the 1944 Tonanki (Japan), 1946 Nankaido (Japan), and 1960 Chile earthquakes were 7.5, 5.3 and 22.6 Chapter 1. Introduction 3 metres respectively (Abe, 1979). If a very large subduction earthquake were to occur in the Cascadia zone, there is almost no doubt that a large tsunami would be triggered. To assess the tsunami hazards along the British Columbia coastline, computations of tsunami generation and propagation were carried out. The object of this investigation is to present a numerical study of the tsunami generation based on the two-dimensional, partly linearized shallow water equations. Complexities of the ocean topography and bathymetry are included by specifying water depths at intervals of 2 and 5 km respectively in two grids. The tsunami propagation into Port Alberni Inlet, where a severe flooding resulted from the 1964 Alaskan earthquake, is also treated by using an one-dimensional model extension. In the following chapters, the tectonic setting of the Pacific Northwest is first discussed, followed by a description of the models being used and their numerical formulation. Results are then presented with a discussion of the tsunami response at various locations along the British Columbia coast. Chapter 2 Earthquake potential in the Pacific Northwest The Cascadia subduction zone off the west coast of Canada and the United States is a 1400-km-long region composed of four segments: the Winona, Explorer, Juan de Fuca and Gorda South plates (Fig.2.1). Subduction takes place about 100 km seaward of Vancouver Island, where the oceanic and continental plates meet. The nature of seismic activities and the physical characteristics across the plate boundary shed light on the likelihood of a large thrust event occuring here. Not all subduction zones experience great earthquakes. Variations in earthquake magnitudes depend on the degree of seismic coupling. When seismic coupling is strong, slip occurs suddenly during earthquakes; when it is weak, the plates glide past each other in a slow aseismic creep. Strongly seismically-coupled subduction zones are those capable of generating great earthquakes. There are two striking features of the Cascadia subduction zone, as pointed out by numerous authors (Heaton and Kanamori, 1984; Heaton and Hartzell, 1987; Rogers, 1988; Dragert and Rogers, 1988; Dragert and Horner, 1989). One is the absence of any historic shallow-thrust earthquakes across the plate interface since European settlements were established there about 200 years ago. Such a long-lasting seismic quiescence is shared by other subduction zones such as southern Chile, southwestern Japan and Colombia, where long periods of time elapse between violent events. 4 Chapter 2. Earthquake potential in the Pacific Northwest 5 Another is that the subducting segments are relatively young, ranging from 6 to 10 million years of age. Young lithosphere is considered to favor strong coupling because of its high buoyancy, permitting a large normal force across the plate interface, and also because of its flexibility, allowing a larger area of contact. Other physical features suggesting that the Cascadia subduction zone may be strongly coupled include the absence of active back-arc basins in the northwestern United States, the gentle dipping of the subducted slab (10° to 15° ), the shallow depth of oceanic trenches, and the smooth topography of the Juan de Fuca plate. More details can be found in Heaton and Kanamori (1984). The possibility of strong coupling between the plates implies an accumulation of elastic strain. The resulting cycle of ground motion pattern may be deduced by referring to the three-stage deformation cycle in Fig.2.2. At the pre-earthquake stage, the oceanic plate is locked to the continental plate while being pushed under, thus building up strain. Submergence is seen to occur seaward of the deformation front. Upon reaching the elastic limit, the strain is released and an earthquake occurs. A sudden uplift occurs along the rupture front due to the rebound of the strained layer. The general displacement pattern at the post-earthquake stage includes both uplift and subsidence of the earth surface as shown in Fig.2.2. This pattern, when transposed to the British Columbia coast, would result in a sea-bottom uplift along the deformation front off the west coast of the Vancouver Island, accompanied by coastal subsidence in the northeastern and southern parts of the island (Fig.2.3). Chapter 2. Earthquake potential in the Pacific Northwest 6 The discussion so far has been limited to a comparison of the Cascadia subduction zone with other subduction zones; it would be far more convincing, in terms of earthquake probability, if direct evidence for the strain accumulation and for prehistoric large earthquakes were available. Some of this information has been provided by recent geodetic surveys carried out in the region. It has been reported (Dragert and Rogers, 1988) that a vertical surface uplift of 2 mm per year is taking place from Campbell River to Neah Bay. Horizontally, contraction occurs at a rate of 2 to 3 mm per year along an axis directed parallel to the plate motion vector. Both of these findings suggest that strain is currently accumulating in a pattern consistent with that discussed earlier. Evidence for prehistoric large earthquakes has also been given by Heaton and Hartzell (1987) and Dragert and Rogers (1988). The deformation pattern given in Fig.2.3 summarizes the results of various authors' work mentioned in this chapter. The four segments which are of slightly different ages might rupture independently, resulting in a sequence of several great earthquakes with magnitudes of about 8. If all segments ruptured as a whole, a giant earthquake of magnitude over 9 would result. In all cases, Dragert and Horner (1989) estimated a maximum uplift of 5 m at the rupture front and a maximum subsidence of 0.5 m in coastal regions. The form of the bottom displacement acting as the source of a potential tsunami is shown in plane view in Fig.2.3a, in a vertical section in Fig.2.3b and in perspective in Fig.2.3c. Chapter 2. Earthquake potential in the Pacific Northwest Figure 2.1: The Cascadia subduction zone, showing the location of four different segments. The rupture or 'deformation' front is denoted by dashed lines. (After Rogers, 1988) Chapter 2. Earthquake potential in the Pacific Northwest 8 DEFORMATION C Y C L E FOR A SUBDUCTION T H R U S T Z O N E 1. STRAIN ACCUMULATION iri . .'-v -* > V. .. .:. / •> 2. THRUST RUPTURE 3. POST-SEISMIC READJUSTMENT £2& Figure 2.2: The three-stage deformation cycle for a subduction zone. The continental plate is shown on the right while the oceanic plate is on the left. (After Dragert and Rogers, 1988) Chapter 2. Earthquake potential in the Pacific Northwest 9 Figure 2.3: (a) Contours of sea-bottom displacement (in cm) with all three segments ruptured as a single unit: the Explorer segment is shown in dashed lines with the Winona segment in the north and the Juan de Fuca segment in the south. Ground motion includes a maximum uplift of 5m along the rupture front and a subsidence of about 0.5 m in coastal regions. Chapter 2. Earthquake potential in the Pacific Northwest Figure 2.3: (b) Two-dimensional cross-section view of the bottom displacement. 10 Figure 2.3: (c) Three-dimensional view of the bottom displacement Chapter 3 Tsunami Generation and Propagation Models To calculate the tsunami response on the British Columbia coast, three models covering different areas were developed: a two-dimensional, 5-km mesh deep-ocean model for tsunami generation and propagation across the continental shelf region; a two-dimensional fine grid model of 2-km grid spacing for wave propagation into Juan de Fuca Strait and the Strait of Georgia; and a one-dimensional inlet model, also of 2-km grid intervals, for tsunami propagation into Alberni Inlet. Numerical formulations of all three models, together with their couplings, are discussed in detail in the following sections. 3.1 The deep-ocean model 3.1.1 Tsunami generation and propagation 3.1.1.1 Tsunami generation The forcing term in the numerical modelling of tsunami comes in the form of water surface elevation caused by the underlying sea bottom vertical displacement. As discussed previously in the tectonic section, the net vertical motion for the hypothetical megathrust earthquake would include an uplift of about 5 m along the rupture front and a subsidence of 0.5 m in the coastal regions of Vancouver Island. In the simulated source region adopted here, the vertical sea bottom displacement 12 Chapter 3. Tsunami Generation and Propagation Models 13 is assumed to occur instantaneously and simultaneously at every depth point across the displacement zone. Strictly speaking, seismic waves which travel at speeds of 3 to 4 km/s take 1 to 2 minutes to propagate from the epicentre across the displacement zone, and vertical displacements which occur once the seismic waves arrive may take 10 seconds or so (Hwang & Divoky, 1970). If such a detailed time history of crustal motion were to be included, the continuity equation for the ocean layer could be written as (v-()t = -(du) -(dv) x y (3.1) where TJ is the elevation of water surface above mean level, d the mean water depth, and £ the bottom displacement. The £ term may however be dropped, as is in our case, without affecting the final tsunami wave form because the tsunami travel time over the bottom displacement area (at a speed of y/g~h) is of the order of 15 minutes, which is far longer than the crustal deformation time. On the time scale of water waves, the seismic disturbance thus appears instantaneous. Once the vertical sea bottom deformation occurs, the overlying water column gets displaced as a direct consequence of mass conservation. This results in a pressure gradient due to the sea-surface slopes, and thus a wave is generated and propagates away from the source. 3.1.1.2 Tsunami propagation A tsunami, with an amplitude small compared with both wavelength and ocean depth, is well described over relatively short propagation distances by the shallow Chapter 3. Tsunami Generation and Propagation Models 14 water equations: rj = -{du) t (3.2) - (dv). x ku\U\jd (3.3) v + uv + Wy = — gn — fu - kv\U\/d (3.4) u + uu + vuy — -gn t x t x x + fv - y where U=(u,v), with u, v the vertically averaged horizontal velocities in the x and y directions respectively, d the mean water depth, f the Coriolis parameter, and k the dimensionless quadratic friction coefficient. Because of the long wave nature of tsunamis, the gradients of the velocity vector in the horizontal direction, and hence the nonlinear acceleration terms in Eqs.(3.3) and (3.4), are anticipated not to be locally important. As to be verified later, the results of linear and nonlinear calculations indeed look almost identical. Because of the computational simplicity of the linear case, the nonlinear terms will be dropped. In summary, the governing equations for the local tsunami response are based on the linearized non dispersive shallow water equations: y (3.5) ku\U\/d (3.6) 7] = ~(du) - (dv) t u t'= x —dVx + fv - (3.7) Chapter 3. Tsunami Generation and Propagation Models 3.1.2 15 The geographical location of the deep-ocean model The configuration of ocean topography dictates the fate of wave propagation; for example, wave reflection, refraction and energy convergence result from changes in water depth. Ocean bathymetry is therefore explicitly included in numerical modelling of wave propagation. Our area of study is basically contained in two grids (Fig.3.1): a coarse grid encompassing the north coast of British Columbia's mainland, the Queen Charlotte Islands and the west coast of Vancouver Island, and a fine grid encompassing Juan de Fuca Strait, the Strait of Georgia and several narrow channels and fjord systems in the mainland. The fine grid model gives a better spatial resolution which is necessary to account for the complexity of coastal structures. Before the essential features of the coarse grid are to be outlined, some justification is offered for choosing a rectangular Cartesian grid,over a spherical polar grid, which better describes the earth's surface. Since only the local effect of the tsunami is going to be investigated, the area of interest is relatively small, ranging from 45.3° N to 56.9° N in latitude. Over this latitude range, there is some measurable curvature to the earth's surface. The approximation of a curvilinear surface by a Cartesian grid thus causes some error in the calculations of surface area and subsequently in wave amplitudes. The error on computing the latter, as is shown in Appendix A, is only about 6% within the specified latitude range. The use of a more complicated spherical polar grid therefore seems to be unnecessary in our tsunami near-field calculation. 3.1.3 The coarse grid The two dimensional coarse grid of resolution 5x5 km is centred at 129.5° W,51.5° 2 N and is oriented at an angle of 30° with respect to the East (Fig.3.1). The grid, Chapter 3. Tsunami Generation and Propagation Models 16 which has 70 columns and 220 rows, was obtained from Seaconsult Marine Research Ltd., which in turn extracted the water depth information from the U.S. National Geophysical Data Centre world bathymetry and topography database. For the present tsunami study, two changes were made to that grid: i/ The tsunami source region which is located approximately around the 2,000 m bathymetric line is quite close to the seaward boundary of the Seaconsult grid: about 20 grid spacings. Spurious reflected waves, as to be discussed later, are likely to arise unless the boundary condition is carefully treated. To delay the arrivals of possible unwanted reflections, the grid was extended to almost twice its original size by attaching an extra section of a flat-bottom ocean of 2,700 m in depth (Fig.3.2). The choice of such a water depth provides smooth continuity with the depth points in the original boundary and also introduces a time delay of about 61 minutes in the arrival of possible reflected waves from the western computational boundary. Such a time delay is found to be sufficient to recognize reflected waves, which would arrive much later than the leading waves of the tsunami. As we are not concerned with a precise calculation of the tsunami in the deep ocean, we need not worry about the artificiality of aflat-bottomedwestward grid extension. ii/ The second change made on the coarse grid is at the Vancouver Island's Juan de Fuca Strait entrance and at the northeast entrance to the Strait of Georgia. Originally two solid boundaries were imposed at these places. This choice is, of course, not applicable when the study area extends beyond the model boundary into the Strait of Georgia. The extension was done by connecting the deep-ocean model with afinegrid model and by driving the latter with the coarse grid results. The choice of boundary Chapter 3. Tsunami Generation and Propagation Models 17 conditions most suitable for matching the two grids at the mouth of the Strait of Juan de Fuca will be discussed in section 3.2.4 below. Chapter 3. Tsunami Generation and Propagation Models 18 Figure 3.1: Location map of the study area for the west coast of British Columbia. The 5-km coarse grid (A) is used for wave propagation modelling in the ocean-continental shelf region, while the 2-km fine grid (B) is for relatively shallower waters. (After Dunbar et al, 1989). Figure 3.2: The coarse grid used in the deep-ocean model. The dots in the grid layout represent vertices of grid elements. The central vertical line denotes the open boundary before extension. A flat bottom ocean of depth 2,700 m is assumed for the extended part. The shaded area represents a test tsunami source of upward displacement of 2 m. Marked stations are chosen in the order of increasing angles of wave incidence onto the open boundary: waves incident perpendicularly on the boundary for stations starting at grid point (1,75); for stations starting at (4,113) and continuing northwards, the incident angles become more oblique. Corresponding elevation time series are shown in Fig.3.4 and 3.7. Chapter 3. Tsunami Generation and Propagation Models 3.1.4 20 Numerical solution of the deep-ocean model The numerical scheme employed to solve Eqs.(3.5) to (3.7) is the explicit onesided difference scheme in space and time using an Arakawa-C-type grid (Fig.3.3). The explicit nature of the scheme requires the Courant-Friedrich-Lewy condition to be satisfied: < y - [gd ((AxY + (A»)')]i ( ) AxA A t 3 8 max A fixed time step of 10 sec is chosen to ensure stability. The finite difference forms of Eqs.(3.5) to (3.7) are summarized as follows: Vij ~ VH At Vij+ijdjj v-i+ij {dij + d ) - u^di-ij 2Ax i+lij + djj ) +1 - 2Ay v (di j_ ij i 1 -f- d^) + d^) (3.9) Primed and unprimed variables indicate those evaluated at the current and previous time steps respectively. Special attention is drawn to the Coriolis and friction terms in the momentum equations. If u, v are calculated in the order shown, the Coriolis term in Eq.(3.10) uses the unprimed v and the corresponding term in Eq.(3.11) Chapter 3. Tsunami Generation and Propagation Models 21 uses the primed u. This scheme is reported by Henry (1982) to be necessary for stability. Similarly, the quadratic friction term in Eq.(3.10) uses the unprimed values to calculate the modulus of U but primed value for the u-component, i.e. (u'\U\)/d, and correspondingly (v'\U\)/d in Eq.(3.11). To reduce possible bias, the model computes variables in the order v,u,v on odd time steps and T/,V,U on even steps. More de- tails on this scheme can be found in Henry (1982). Because bottom friction is not as important in the deep-ocean and shelf model as in the shallow near-shore areas, a constant friction coefficient of 0.0025 is assumed at all calculation points. 3.1.5 Boundary conditions At closed boundaries, a zero velocity is assumed, and hence no volume transport is allowed. This boundary condition is equivalent to complete reflection of the wave at the coastline. At the seaward open boundaries, the problem is not obvious because the boundaries do not exist in nature. Much work has been done, for instance, by Engquist (1977) and Camerlengo (1980) to ensure that the solution in the model interior is affected as little as possible by spurious reflection from the artificial open boundary. The approach taken by Henry (1982) is to let waves radiate out of the model domain from the interior: the standard radiation condition. He made the assumption that the velocity vector could be treated one-dimensionally on the boundary, and hence the radiation boundary condition reads u(or v) = depending on the direction of propagation. ±(^)^, (3.12) Chapter 3. Tsunami Generation and Propagation Models 3.1.6 Model testing and reliability 22 * The deep-ocean model applied on the coarse grid, as just mentioned, is adopted directly from Henry (1982), who applied it successfully to the Beaufort Sea, and hence sensitivity testing of the numerical scheme was not carried out extensively, except for the open boundary conditions to be discussed below. The model has also been used in previous tsunami applications, for example, by Murty and El-Sabh (1985) in the St. „ Lawrence estuary, and by Crean et al (1988) in modelling the 1946 earthquake near Comox, Vancouver Island. As a test source, an area of 30x30 km was taken to be 2 9 grid spacings away from the original seaward boundary (Fig.3.2); an instantaneous uniform uplift of 2m was assumed to occur over that whole area. Mass conservation was checked, and it was found that transports across the model boundaries satisfy continuity. The square area source was also used to test boundary conditions. 3.1.6.1 Testing of open boundary conditions As discussed earlier, the open boundary condition, as given in Eq.(3.12), should let waves pass out of the model domain. Complete transmission is, however, achieved only for waves incident perpendicularly on the boundary; unwanted reflection may arise for other angles of incidence. To see the significance of any spurious reflections under the 1-d radiation condition (Eq.3.12), the open boundary was shifted away from the test source by extending the grid westward (Fig.3.2), and results are compared before and after the extension (Fig.3.4). As indicated in Fig.3.2, a number of stations are chosen to test reflection at almost normal or oblique angles. Significant reflections are observed in the northern part of the grid (Fig.3.4c) due to the dominance of the velocity component parallel to the boundary, which is neglected in the 1-d condition (Eq.3.12). To suppress such unwanted reflections, Eq.3.12 is modified to include both C h a p t e r 3. T s u n a m i G e n e r a t i o n and P r o p a g a t i o n M o d e l s 23 velocity components. Thus at the western boundary, U = ±TJ^T) 2 — v (3.13) - v (3.14) 2 and in the northern and southern boundaries, v =± ^ T / 2 2 The ± signs of Eqs.(3.13) &; (3.14), together with the equations' finite difference forms, are described in detail in Fig.3.5. Attempts have been made to implement the extended radiation condition (3.13) and (3.14) on the appropriate boundaries. Instability was however experienced on the southern boundary when the square root in Eq.c (Fig.3.5) became imaginary. Why this instability occurs may be explained as follows. In Eq.c (Fig.3.5), the boundary variable t;^ is calculated from values in the interior, e.g. from Vi , 2 un, Ui+ij. These interior variables are in turn evaluated in Eqs.(3.10) & (3.11). Recalling the earlier discussion, the Coriolis and friction terms in Eqs.(3.10) & (3.11) are expressed in terms of both primed and unprimed values of u and v. The mixing of the new and old values may cause small phase shifts. This small phase difference may become important when dealing with quantities such as those inside the square root of Eq.c (Fig.3.5) where rjn, un and u i+1<1 have to be at the same time level. Results of elevation time histories when applying the 2-d boundary conditions (Eqs.3.13 & 3.14) are given in Fig.3.6 for the first 90 minutes (instability occurs around 100 min). The discrepancies observed previously at the grid's northern region (Fig.3.4c) are obviously removed to a large degree under the present 2-d radiation Chapter 3. Tsunami Generation and Propagation Models 24 scheme. The reflection off the western boundary seems to be greatly suppressed. Because of the numerical instability arising from the 2-d radiation condition, Eq.(3.12) is used instead for later model runs. Spurious reflections are delayed by extending the coarse grid westward as discussed earlier. Stable 2-d radiation conditions need to be developed in future studies in order to carry out computations more efficiently. 3.1.6.2 Comparison of linear and nonlinear terms It was anticipated earlier that the nonlinear terms in Eqs.(3.2) to (3.4) were not locally important. By using the test source in Fig.3.2, elevations were computed using linearized and non-linearized equations. The percentage difference of the two cases is shown in Fig.3.7 for some selected stations. This minimal difference in computed elevations is observed at most places, and hence the use of the linearized equations (Eqs.(3.5) to (3.7)) is justified. Chapter 3. Tsunami Generation and Propagation Models AX- Uj-|,j*| + Ui,j | + 0 + + Vi-i.j+i Vij+i X + 4- o UJ4l,j+l 3 I+,J+ ' + Vi+| j+| X f + o + + O + Vij X - Vi-1,1 x —— f 0 X o Ui-| j-| 77JJ+J O V M j H y.v, j i X,U,I o + x 17, elevation point u, velocity point v, velocity point Figure 3.3: Arakawa-C-type grid. (After Henry, 1982). Chapter 3. Tsunami Generation and Propagation Models 26 Figure 3.4: A comparison of computed elevations between the original grid (solid line) and the extended grid (dashed line) as a test of spurious reflections off the open boundary. Locations of stations are indicated in Fig.3.2. In (a), as the waves are incident normally, the leading waves almost coincide for the two cases, and sucessive waves also agree quite well. As the incident angles increase, as in (b), small differences are observed both for the leading and subsequent waves. In (c), as the angles deviate more from normal incidence, the entire wave train is modified by the spurious reflections. The vertical scale is in metres. Chapter 3. Tsunami Generation and Propagation Models Figure 3.4: Continued. 27 28 Chapter 3. Tsunami Generation and Propagation Models o" (20,139) ID (17,135) O (15,131) in (12,127) in o. i 0.0 30.0 60.0 90.0 T(MIN);LI NEAR: (b) Figure 3.4: Continued. 120.0 150.0 180.0 Chapter 3. Tsunami Generation and Propagation Models Figure 3.4: Continued. 30 Chapter 3. Tsunami Generation and Propagation Models l>i,221 *>i,221 = ±2y^/di,22cV?, ,220 ? ( .\220 + «i+l,22o)'/4 - _ u square root sign + if i?i,22o > 0 - uij = V 20 ii2 (Eqb) if Vi,220 < 0 ±2^ld )m\ x d - {v 1J+l Ut+l,220 Ui,220 Ui,220 + vuY/4 - u 2J (Eq.a) square root sign + if vij < 0 - if ViJ > 0 = ±2^/<*i.^, .i square root sign + if - if ? Ki+ ) / ~ »i.a 2 4 ( q- ) £ c < 0 >0 Figure 3.5: Finite-difference forms of the 2-d radiation boundary condition on the three seaward boundaries. 31 Chapter 3. Tsunami Generation and Propagation Models in o (21,169) - S ^ i v in o (20,169) in (19,169) CE > UJ ^^s^^^v (18,169) A^w^ in o" (17,169) ^Vw^/V^ in o I 0.0 30.0 I 60.0 I 90.0 T(MIN) ;LINEAR: 1 120.0 COMPARE 1 150.0 1 180.0 Figure 3.6: A comparison of computed elevations for three different boundary cases: original grid with 1-d radiation boundary conditions (solid line), original grid with 2-d conditions (dashed line), extended grid with 1-d conditions (dotted line). Clearly, the use of the 2-d boundary condition approximates the extended grid results better. The Coriolis terms are dropped in computations when applying the 2-d conditions. Stations are located in the northern half of the grid, as shown in Fig.3.2. The vertical scale is in metres. Chapter 3. Tsunami Generation and Propagation Models T 1 M I N I ; L INEflR L 32 NONLINEAR Figure 3.7: The percentage difference (vertical axis) between calculated elevations using linearized and non-linearized equations. Station locations are indicated in Fig.3.2 Chapter 3. Tsunami Generation and Propagation Models Figure 3.7: Continued. 33 Chapter 3. Tsunami Generation and Propagation Models 3.2 34 The fine grid model After the tsunami propagates into Juan de Fuca Strait, it goes through a complex series of narrow islands and channel divisions at the southern and northern entrances of the Strait of Georgia, where damping and nonlinear interactions might result. A relativelyfinergrid model than the deep-ocean model is obviously necessary to account for the irregularities of the topography and coastline configuration. A 2-km mesh vertically-integrated model adapted directly from the GF7 model of Crean et al (1988) serves this purpose. Thefinegrid model layout is shown in Fig.3.8. In practice, however, a much more compact array (Fig.3.9) is used, where empty spaces are closely filled with segments of the actual array. Junction points are denoted by letters as shown. 3.2.1 The governing equations The fine grid model is basically similar to the deep-ocean model, but some differences do exist. The vertically-integrated equations of continuity and momentum are (3.15) dqi 9q 2 0 9 q\ 8 qq x 2 d qq x q\ 2 fq2+9(d + v)^ L + fq (d 1+g + + r,)^- + kq\\q (d + v) 2 kq \q\ 2 (d + v) 2 0 (3.16) 0 (3.17) Chapter 3. Tsunami Generation and Propagation Models 35 where 171 = u(d + rj), tj2 = v(d + 77), with u,v the vertically averaged velocities as defined in the deep-ocean model. The previous assumption of (d + rj) — d in the deepocean model does not apply here because the magnitudes of 77 and d may become comparable in shallow waters. This difference results in additional nonlinear terms, such as, in the momentum equations: gwT) x and gnvy. The friction terms also change to (kq \q\)/(d + rj) and (kq \q\)/(d + r}) . The advective accelerations, though included 2 2 1 2 here, are later seen to have very little effect on the solutions as is the case in the deep-ocean model. Another major change from the deep-ocean model is the variability of bottom friction coefficients. Since friction effects and their local variation are more important in shallow depths, a regional distribution of dissipation was allowed by varying coefficients of bottom friction. As discussed in Crean et al (1988), byfittingthe computed amplitudes of the M tidal constituent to the measured values, the overall friction 2 coefficient was set to 0.003, with the following exceptions: the Puget Sound system and Rosario Strait (0.006); Haro Strait(0.03); Active, Porlier, Gabriola Passes, and Dodd Narrows in the southern Gulf Islands (0.03); Seymour Narrows (0.006), Okisollo Channel (0.03), and Yuculta Rapids (0.03) of the northern passages. Locations of the channels and passes are shown in Fig.3.8. 3.2.2 Numerical scheme The numerical scheme employed here is again similar to that in the deep-ocean model: explicit one-sided differencing in space and time. However, slight differences are also present. Thefinite-differenceforms of Eqs.(3.15) to (3.17), based on the computational grid in Fig.3.10, can be found in Crean et al (1988). We point out that Chapter 3. Tsunami Generation and Propagation Models 36 both the friction and Coriolis terms here use the 171 and qi computed at the previous time step, whereas in the deep-ocean model, these terms, as emphasized previously, involve quantities in mixed time steps (Eqs.(3.10) and (3.11)). To satisfy the CFL criterion (Eq.(3.8)) for stable solutions, a fixed time step was chosen to be 20 sec. 3.2.3 B o u n d a r y Conditions Coastal boundaries are characterized, as usual, by zero normal components of velocity vectors. Open boundaries are taken along the junctions between the coarse and fine grids where prescribed elevations are obtained from the deep-ocean model. 3.2.4 Junctions between the coarse and fine grids Junctions of the coarse andfinegrids are located at Juan de Fuca Strait in the south and Johnstone Strait in the north (Fig.3.1). The fine grid model is driven by the coarse grid junction values, and thus the model accuracy within the Georgia StraitPuget Sound system is affected by the nature of boundary conditions chosen for the coarse grid. The influence of the shape of the boundaries was first investigated. By using the test source of uniform uplift of 2m in Fig.3.2, two Juan de Fuca Strait open boundaries have been tested on the coarse grid: one running through the u-component of the velocity vector and one running alternately through the u, v points as shown in Fig.3.11. Time histories of elevation at the neighborhood of the junction are plotted for the two cases (Fig.3.12). Results of the original solid boundary are also shown for comparison. It was found that the two open boundaries show very little difference. The boundary through the u-points was chosen for later model runs because of simplicity Chapter 3. Tsunami Generation and Propagation Models 37 in coding. Secondly, to reduce possible spurious reflections from the open boundaries and to have sufficient elevation points covering the junctions, the u-point boundaries were extended eastward to include Johnstone Strait in the north and a major part of Juan de Fuca Strait in the south. The resulting grid is shown in Fig.3.13 and Fig.3.14. Comparisons of modelled water levels for the extended and original coarse grids are given in Fig.3.15. It is observed that the results agree fairly well for the two cases even for locations which are only 3 to 4 grid spacings away from the original boundary. The present extension scheme, though only resolving the Juan de Fuca Strait dynamics partially, should provide enough accuracy for the junction points, which are more than 14 grid intervals from the extended boundary, in driving the fine grid model. It is assumed that the deep-ocean and fine grid models are coupled independently; that is, the deep-ocean model is operated as if the waves reflected from the fine grid domain did not enter the coarse grid region. This assumption is adequate if the tsunami is significantly absorbed, rather than reflected, within the Strait of GeorgiaPuget Sound system. Prescribed boundary values on thefinegrid model are obtained directly from the deep-ocean model and are interpolated spatially and temporally before application. The spatial arrangement of the fine and extended coarse grids at the Juan de Fuca Strait junction is shown in Fig.3.16. If vector quantities (tjj or q velocities) were 2 chosen to be the boundary values, interpolations would become complicated because of the angle between the grids. Instead elevations are chosen to be prescribed on the open boundaries. A four-point linear interpolation scheme is used to estimate the fine Chapter 3. Tsunami Generation and Propagation Models 38 grid boundary elevation points, as is illustrated in Fig.3.17. At the northern passage, the junction consists of only one grid point and is set to the nearest elevation point available in the coarse grid. No spatial averaging is necessary. The stored coarse grid data which are given at 2-min intervals are linearly interpolated so as to be compatible with the fine grid calculations of every 20 sec. Chapter 3. I 3? Tsunami Generation and Propagation Models 5= 3? 39 I Figure 3.8: Geographical location of the fine grid model. (After Crean et al, 1988). Chapter 3. Tsunami Generation and Propagation Models 40 Figure 3.9: Computational array of the fine grid model. (After Crean et al, 1988). Chapter 3. Tsunami Generation and Propagation Models • £i-n-l - V i-n-l U i-n-l 41 f £i-n U i-n U i-n*l i-n*l i - n - V; i-l J+r» U J- V y.v. [ i+n-r Un Cl+n+l l+n+l U i+n+l V x,u. Figure 3.10: Numerical scheme used in thefinegrid model. (After Crean et al, 1988). Chapter 3. Tsunami Generation and Propagation Models 42 Figure 3.11: Coarse grid layout at the Juan de Fuca Strait entrance, (a) Location of the junction between the coarse and fine grids, (b) Two configurations of open boundaries are tested at the Strait entrance of the coarse grid. Chapter 3. Tsunami Generation and Propagation Models 43 Figure 3.12: A comparison of sea levels calculated from implementing the open boundary through the u-v points (dashed line) and through only the u-points (dotted line). The difference is perceivable only for some brief moments at point (67,18). Results of the original solid boundary are shown in solid line. In panel (d), only the results through u-points are shown. Locations of stations are indicated in Fig.3.11a. The vertical scale is in metres. Chapter 3. Tsunami Generation and Propagation Models 44 •o. CD ' CD (C) fM O . 0.0 20.0 -I 40.0 1 60.0 TIME in -o—o-o-o- 1 80.0 ~ i 100.0 (MIN) s o l i d boundary open boundary on u p o i n t s 120.0 only CD° CD in CM (d) in 0.0 20.0 I 40.0 TIME 1 60.0 (MIN) Figure 3.12: Continued. 1 80.0 l 100.0 1 120.0 Chapter 3. Tsunami Generation and Propagation Models 45 Figure 3.13: The extended coarse grid layout. Solid lines denote the original boundaries(Fig.3.2), whereas the dashed lines indicate the line grid boundaries. In order to have a sufficient number of depth points (and thus water surface elevation points) covering the junction and also to reduce possible spurious reflections from the open boundaries, the coarse grid was extended in two entrances: Juan de Fuca Strait and Johnstone Strait. Because the width of Johnstone Strait could not be resolved within the coarse grid dimensions, the extension made here is not exact. The approximation is however justifiable since the junction consists of one grid point only. Johnstone Strait is assumed to have a depth of 180m. Chapter 3. Tsunami Generation and Propagation Models 46 > I Figure 3.14: T h e J u a n de Fuca Strait junction between the coarse and fine grid models. J u a n de Fuca Strait is extended as far up the Strait head as possible i n the coarse grid model. For 129 < I < 133, depth points are extrapolated from the fine grid data as obtained from Crean et al (1988). For 134 < I < 143, all points are given a constant depth of 150 m . Such an approximation is acceptable as this part of the region will be repeated in the fine grid calculation. Chapter 3. Tsunami Generation and Propagation Models 47 (M O Figure 3.15: A comparison of computed elevations (in m) in the Juan de Fuca Strait area. Solid and dashed lines denote respectively the results before and after the grid extension. Station locations are indicated in Fig.3.14. Chapter 3. Tsunami Generation and Propagation Models Figure 3.15: Continued. 48 Chapter 3. Tsunami Generation and Propagation Models 49 Figure 3.16: Junction of the fine grid and extended coarse grid at the Juan de Fuca Strait. Chapter 3. Tsunami Generation and Propagation Models 50 Ax Ax Vf = <f>2+n((p -(j) )l/\ l where <pi = rj-i 4>2 — 2 x +m(v —T]i)/Ax 2 Vi+nT-iVs-V^/Ax Figure 3.17: Four-point interpolation scheme used in calculating the fine grid open boundary values. Chapter 3. Tsunami Generation and Propagation Models 3.3 51 Tsunami propagation in inlets Wave magnification in coastal waters arises from a number of factors, such as re- flection, refraction, funnelling and resonance. Relatively large wave amplitudes may result if some of these factors come into effect together. An example in the west coast of Canada is given by Port Alberni on Vancouver Island. The City of Port Alberni which is located at the head of Alberni Inlet has experienced a number of tsunamis in historic times, with accounts dating back to 1896 (Marshal Macklin Monaghan, 1986). Of particular interest is a 4-m tsunami triggered by a distant source during the 1964 Alaskan earthquake. The damage resulting from this event totalled $5 million, representing the most destructive tsunami recorded on the west coast of Canada and U.S.A.. The large-amplitude waves observed at Port Alberni are due to the combined effects of the topographic amplification in Barkley Sound and the resonance magnification in the inlet (Marshal Macklin Monaghan, 1986; Murty and Boilard, 1970). Of the numerous urban settlements located at the inlet heads along the British Columbia coast, Port Alberni poses the highest risk to tsunami threats. For this reason, Alberni Inlet is chosen for some detailed study. Fig.3.18 shows the location of the City of Port Alberni. For simplicity, the inlet is represented one-dimensionally by a straight line. This representation is justified on the grounds that the width of the inlet is much smaller than its length. Branchings of channels are neglected in the study. It is further assumed that the inlet model is operated independently from the deep-ocean model. Chapter 3. Tsunami Generation and Propagation Models 3.3.1 52 The inlet model The inlet model is based on the one-dimensional, depth averaged equations of continuity and momentum (3.18) Wr,t + QX = 0 (3.19) where Q is the transport (m /s) given by Q=uWd, with u the vertically averaged s velocity as defined in the deep-ocean model, W the surface width of the inlet, d the water depth, A the cross sectional area of the inlet, the frictional shear stress, and p the density of water. These equations have been applied by several authors in previous studies of tsunami response in inlets, as referenced in Dunbar et al (1989). The geometry of the inlet is described in details by specifying widths, depths, surface and cross-sectional areas at model calculation points. Because of the importance of friction dissipation in shallow waters, friction coefficients are allowed to vary at each grid point. 3.3.2 Numerical solutions Eqs.(3.18) and (3.19) are solved numerically in space and time using an explicit, one-sided differencing scheme, as in the previous two models. A time step of 20 sec was chosen to ensure solution stability. Based on the staggered grid in Fig.3.19, the numerical forms of Eqs.(3.18) and (3.19) are given as follows: 53 Chapter 3. Tsunami Generation and Propagation Models l + (3.20) ( At\Q?\)/(A?C?d?) 9 A*(Q? ~ Qi-i) +1 (3.21) Si-l where Q is the transport (m /s), A? the cross-sectional area (m. ) at a Q point [A™ = 3 2 A® + Wi(rj™ +r)i')/2] with A° the cross-sectional area at low water and W{ the surface +1 width at a Q point; d? is the depth (m) at a Q point [d? = d° + (rjf + V?)/^} ^h w +1 d° the mean depth; Si is the inlet surface area (m ) between adjacent Q points; Ci is 2 the Chezy friction coefficient at a Q point, Ci = yu\u\pg/rh, with the frictional stress Th given in Eq.(3.19). The bottom friction coefficients are adjusted by matching the computed tidal response of the M , S and K\ constituents with the measured values 2 2 (Dunbar et al, 1989). 3.3.3 Boundary conditions 11 Closed boundaries As in the previous cases, a condition of zero normal velocity is imposed at the closed boundary. This complete-reflection assumption is, however, true only for steep coasts. On gently sloping beaches, as is our case, the wave runs inland and dissipates much of its energy through inundation of the surrounding land. As pointed out repeatedly by Dunbar et al (1989), the zero-velocity boundary condition would over-estimate the computed water-levels at the inlet head by neglecting the wave dissipation from flooding. The accuracy of modelling wave response at the inlet head could be improved by taking into account the run-up process, including possible bore formation and Chapter 3. Tsunami Generation and Propagation Models 54 rushup, through use offineresolution models with more realistic boundary conditions, for instance, the moving water line condition. ii/ Open boundaries Prescribed open boundary values at the mouth of the inlet are provided by the deep-ocean model at Barkley Sound. The stored coarse grid data are linearly interpolated tofixedtime steps of 20 sec. It is noted that because of the narrow width of the inlet compared to the wavelength, not all wave energy at Barkley Sound would propagate up the inlet. As a result, modelled wave amplitudes along the inlet could be improved by separating the wave at Barkley Sound into incoming and reflected outgoing components. These improvements are left to further studies. Figure 3.18: The geographical location of Alberni Inlet. (After Marshal, 1986). Chapter 3. Tsunami Generation and Propagation Models INLET HEAD n.t INLET MOUTH -V n+l t=(n+l)At V\ 56 Qr 1 n+ l Vi Q? +1 i,x v?£ Vl9 n+l Ql9 — i — At _n Qi Vi Q? t=nAt Vl9 —e- Ax I Prescribed elevations Figure 3.19: Numerical scheme used in the inlet model. 19 Chapter 4 Results Three cases of sea-bottom displacement in the Cascadia subduction zone were studied for the tsunami generation and propagation problem: the rupture of, first, the three upper segments all at once; second, of only the largest segment, the Juan de Fuca segment; and, third, of both the Winona and Explorer segments but not of the Juan de Fuca segment. A rupture of all three segments would represent the worst scenario studied here, which is equivalent to an event of seismic moment magnitude (Mw) of over 9.2 (Rogers, 1988). This was, however, reported to be unlikely as different subducting segments are of slightly different ages and may rupture independently. If the rupture is confined to the Juan de Fuca segment, it would cause a magnitude 9.1 earthquake. If the Winona and Explorer segments rupture together, a magnitude 8.7 event would result. All three cases include regions of uplift and subsidence, the first case of which is shown in Fig.2.3, with a maximum elevation of 5 m along the deformation front and a maximum 0.5 iri submergence along the northeast and southwest coast of Vancouver Island. The length of the deformation front is about 500 km for the first case, 290 km for the second and 210 km for the third. Numerical simulations using the three models discussed in the previous chapter were performed for the associated tsunami problem, and results are presented in the 57 Chapter 4. Results 58 forms of (a) a first-wave travel time chart to give an overall view of the propagation of the tsunami wavefronts; (b) a spatial map of the leading wave amplitudes to show the spatial decay of wave energy when there are no coastal interactions; and (c) plots of wave amplitude time-series at various coastal locations to illustrate the local tsunami response. The selected regions include some major tide-gauge locations and populated areas along the western British Columbia coast. Finally, three-dimensional snapshots of sea surface displacement in the deep-ocean model region are given at selected time intervals to provide another view of the propagating wavefronts. The tsunami leading-wave travel time charts are shown in Fig.4.1 for the first two cases of bottom displacement. Two things are evident from the figures: (i) the waves are refracted as they approach the shores; this wave refraction is caused by the decreasing water depth, as predicted by Snell's Law. (ii) Waves travel more slowly in shallow water, as is observed both on the western coast of Queen Charlotte Island where circular deep-ocean wavefronts are clearly distorted and in Hecate Strait area where contours are more closely packed. 4.1 Rupture of the Juan de Fuca, Winona and Explorer segments Time-series plots of simulated tsunami waves for the first case of bottom disturbance, that is, the rupture of all three segments as a single unit, are given in Figs.4.3, 4.4, 4.7 & 4.9. Several stations, with their locations shown on Fig.4.2, are chosen here for discussion. Fig.4.3 shows the sea-level variation at a point located above the rupture front. The elevation of water level is seen to drop rapidly from the initial 5m uplift as the Chapter 4. Results 59 wave energy propagates away, and the sea surface appears to be relatively calm after twelve hours, confirming the lack of reflection from the seaward edge of the numerical model. The rupture front is located around the 2,000-m bathymetric line, and the wave response discussed here displays some common features usually observed at other deep-ocean points in the model: a broad spectrum, without resonances and with small amplitudes. Fig.4.4J shows the time-series of sea level displacements at Cape St. James, whose high frequency appearance is very similar to Fig.4.3. Cape St. James, situated at the southern tip of the Queen Charlotte Island, is exposed directly to the open sea and has a narrow continental shelf that keeps wave shoaling to a minimum. Amplification of the wave through resonance does not take place. The response at Cape St.James is therefore very similar to that of the open ocean. In general, wave amplitudes, without coastal interactions, decrease with distance from the source due to the spatial dispersion of wave energy, as is shown in Fig.4.5. But wave amplitudes are seen to increase again, as in Rennel Sound and Queen Charlotte City (Fig.4.4H &; I), both located north of Cape St. James, when the waves reach the shallow shorelines. This amplitude build-up is commonly observed in coastal regions where the continental shelf is wide enough to cause wave shoaling. Shoaling results from a decrease of water depth and a slow down in wave velocity which, when conserving energy flux, would imply an increase of wave amplitudes. On the other hand, waves may also be amplified by refraction. As stated earlier, refraction bends the wave when approaching the shore. The wave number vector continuously changes direction so as to approach the shore normally and so does the energyflowvector. Chapter 4. Results 60 Energy can therefore be focussed. These two effects may explain responses of some other locations along the western Vancouver Island, for instance, Bamfield and Torino (Fig.4.4B and D). In addition to shoaling and refraction, wave amplitudes at Barkley Sound and Quatsino Sound (Fig.4.4C and E) appear to be further enhanced by resonant oscillations. By comparing Fig.4.4El and E2, it is evident that certain frequencies are preferentially excited after the tsunami enters Quatsino Sound. An eigenmode of about 60 minutes seems to be present in both Barkley Sound and Quatsino Sound. Port San Juan which is located at the northwest entrance of the Strait of Juan de Fuca shows well tuned oscillations of relatively high frequency (Fig.4.4A). Seiche excitation in Port San Juan has been studied by Lemon et al (1979), who reported a fundamental resonant period of 15 minutes. This is the period which is prominent in Fig.4.4A. After the tsunami propagates into the Strait of Georgia, the wave amplitude is greatly attenuated partly because the wave energy continues to be dispersed over a greater area along the path of propagation and, additionally, because energy is lost to friction which becomes important in shallow waters, notably in the complex Gulf Islands system. Such a decrease in amplitude is seen in the sequence of time-series from Victoria (Fig.4.7.4) to Vancouver (Fig.4.7.13) where the wave is attenuated by a factor of two. The tsunami action near Vancouver, though small, is still noticeable, with an amplitude of about 60cm. As the oscillations continue for many hours, a tsunami would be superimposed upon a high-tide and might pose some danger to low-lying areas of the Eraser River delta. Chapter 4. Results 61 In general, elevations in the Strait of Georgia are below a metre. An apparent resonance in Burrard Inlet deserves further consideration. Fig.4.7.16 shows a growing, well tuned oscillation in the inner basin of Burrard Inlet, suggesting the presence of resonance at a period of approximately 3 hours. A simple but adequate model of tides in Burrard Inlet has been presented by Baines (1958) and may serve as a guide to understand the resonance phenomenon. It is seen in Fig.Bl of Appendix B that Burrard Inlet is composed of three basins inter-connected by two narrow channels. Because of their restricted entrances, the two inner basins behave like Helmholtz resonators. The oscillations are controlled by friction through the two narrows and the geometry of both the basins and channels. Following the work of Baines (1958), a resonant period of approximately 4.1 hours is obtained using the constants suggested in the paper, with details of calculations outlined in Appendix B. Considering the dependence of the resonant period on the exact geometry and frictional values used, this resonant period is in fairly good agreement with the one (3 hours) computed in the GF7 fine grid model. It is noted that the nonlinear advective terms, though included in the fine grid GF7 model (Eq.3.16 and 3.17), play only a small role in computing the water level elevations, as is illustrated in Fig.4.8, where calculations with and without those terms are compared. The maximum positive and negative sea surface displacements in the deep-ocean model are 8.5 and -9.6 m respectively, occuring near Quatsino Sound (Fig.4.4Ml and M2). The corresponding values for thefinegrid model are 6.8 and -6.4 m, close to the northern junction (Johnstone Strait). Chapter 4. Results 62 The propagation of the tsunami into Alberni Inlet is determined by applying the inlet model, with results shown in Fig.4.9. The wave height at the head of the inlet is seen to be amplified about three times above that at the mouth in Barkley Sound. By comparing the maximum wave heights at the head of Alberni Inlet and Tofi.no (Fig.4.4D), the ratio is approximately 4 to 1. The same ratio for the 1960 Chilean earthquake was reported by Marshall Macklin Monaghan (1986) to be 3.2. The resonance characteristics of the inlet were also studied by Murty and Boilard (1970) with and without the inclusion of some narrow channels along the inlet. By treating the inlet as a system by itself, the fundamental mode of the free oscillations was found to be 85 minutes. In Fig.4.9b, a resonance at about 75 min appears to be present. The exact period calculated for the resonance depends rather critically on the exact geometry and frictional values chosen. It is noted that the wave height at the inlet head might be over- estimated for the following reasons: i/ the total-reflection boundary condition at the coast does not allow the dissipation of wave energy through running up the beach; ii/ the wave energy at Barkley Sound, where it was taken to drive the inlet model, was assumed to propagate into the inlet without reflection. Reflection back to sea would remove some of the energy. 4.2 Rupture of the Juan de Fuca segment alone In the second numerical experiment, only the Juan de Fuca segment is ruptured (Fig.4.10). The forcing terms at the southern half of the Vancouver Island are therefore identical to the first case, resulting in wave reponses (Fig.4.11A to D) that are Chapter 4. Results 63 very similar to those seen i n Fig.4.4A to D . T h e stations i n the north, on the other hand, experience a further drop i n amplitudes due to the larger spatial attenuation ( F i g . 4 . 1 1 E l to J and F i g . 4 . 4 E l to J ) . Particularly in Rennel Sound and Queen Charlotte City (Fig.4.11H and I), the waves are seen to be almost twice as small as earlier (Fig.4.4H and I). The largest waves are found near the Tofi.no area, with maximum elevations of 6.7m and depressions of sea level of 9.0m. Time series at the locus of the m a x i m u m positive wave is shown i n Fig.4.11M3, with the station location indicated in Fig.4.2. Because of the similar wave signatures i n the south compared to the first case of the source displacement, calculations i n the fine grid G F 7 and i n the inlet models are not repeated i n this run. Fig.4.12 shows a comparison of the junction points at the J u a n de Fuca Strait entrance for the two cases, confirming the similarity of tsunami effects. 4.3 Rupture of the Winona and Explorer segments In the third calculation, only the W i n o n a and Explorer segments are assumed to rupture together. The forcing terms are then located in the Vancouver Island's northern half, and hence the wave characteristics there (Fig.4.13El to J ) look similar to F i g . 4 . 4 E l to J . In the southern half, the waves arrive at a later time and, as expected, are smaller i n amplitude (Fig.4.13A to D and Fig.4.4A to D ) . T h e Strait of Georgia system and the Alberni Inlet follow the same trend: a later arrival of waves with smaller amplitudes; the patterns of the waves are, however, Chapter 4. Results 64 preserved in general (Fig.4.14 and 4.15). The largest waves in this case occur at the same locations as in the source case (i), with magnitudes 8.5 and -11.3 m in the deep-ocean model area(Fig.4.13Ml and M2) and 6.5 and -6.0 m in the fine grid model area. In Fig.4.16, three-dimensional snapshots of sea-surface displacements in the deepocean model domain are shown at selected time intervals for this particular source rupture. These pictures complement the travel time curves of Fig.4.1 and give another view of the propagating wavefronts. The large amplitudes are seen to concentrate along the coast of the Vancouver Island, as described in the local time series (Fig.4.13); the seaward propagation is now evident. 4.4 Sources of different magnitudes: a test of model linearity Since the governing equations for the tsunami generation and propagation are linear, doubling of the initial wave height would result in final amplitudes being twice as large. This expectation was confirmed by repeating the first case of numerical calculations, the rupture of all three segments together, with now a source magnitude one-half the earlier one. Results (Fig.4.17) show that the amplitudes are generally decreased by a factor of one-half with a few exceptions, notably at Quatsino Sound where a small phase shift is observed, but the general wave form is preserved. This small discrepancy may arise from the small nonlinear bottom drag. Furthermore, the linearity is revealed by adding results from the second and third cases of the ground motion and comparing with the equivalent counterpart, namely, the first case when Chapter 4. Results 65 all three segments rupture as a single unit. Fig.4.18 summarizes the equivalence of the two situations. Chapter 4. Results 66 (a) Figure 4.1: Tsunami travel times (in minutes) in the deep-ocean model domain, (a) All three segments rupture together, (b) Only the Juan de Fuca segment ruptures. Contour intervals are 10 minutes. Figure 4.2: The coarse grid layout showing the locations at which time-histories of water level are plotted. Station locations are mostly presented in geographical order, beginning at Port San Juan and going north. Chapter 4. Results 69 Figure 4.3: Computed elevations (in m) of one of the source points located above the deformation front for the case of all three segments ruptured together. The initial elevation is 5.0m. Station location is shown on Fig.4.2. Chapter 4. Results 70 Figure 4.4: Time-series of computed elevations in the deep-ocean model for the case of all three segments ruptured together. Station locations are indicated in Fig.4.2. The vertical scale is given in metres. Chapter 4. Results 71 o in ~i j j • n 0.0 F.Port Hardy l 2.0 l 4.0 l 6.0 l 8.0 T(HR):HLL SEGMENTS RUPTRE Figure 4.4: Continued. l 10.0 I 12.0 J «Cape St. James I.Queen Charlotte City H.Rennel Sound 1 0.0 1 2.0 1 1 1 4.0 6.0 8.0 T(HR):RLL SEGMENTS RUPTRE Figure 4.4: Continued. 1 1 10.0 12 Chapter 4. Results Chapter 4. Results 74 Figure 4.5: Spatial attenuation of wave energy. Leading wave amplitudes, given in centimetres at grid points located at the lower-left hand corners of the numerals, are seen to decrease with distance from the source. The deformation front is indicated with an uplift curve of 500 cm. Elevations at coastal points may be amplified and are not included here. r Chapter 4. Results 75 Figure 4.6: The fine grid layout showing the locations at which time-histories of water level are plotted. Station locations are presented in geographical order, beginning at Juan de Fuca Strait and going into the Strait of Georgia-Puget Sound system. Chapter 4. Results 76 Figure 4.7: Time-series of computed elevations in the fine grid model for the case of all three segments ruptured together. Station locations are indicated in Fig.4.6. The vertical scale is given in metres. Chapter 4. Results 77 l6.Burrard I n l e t 14.Point Atkinson • H 0.0 I 2.0 1 4.0 I 6.0 I 8.0 T(HR):GF: RLL SEGMENTS Figure 4.7: Continued. I 10.0 I 12.0 Chapter 4. Results Figure 4.7: Continued. Chapter 4. Results Figure 4.7: Continued. Chapter 4. Results 80 Figure 4.8: A comparison of the fine grid model results with and without the nonlinear advective terms, denoted respectively by solid and dashed lines. Station locations are indicated in Fig.4.6. The vertical scale is given in metres. Chapter 4. Results 81 7. Sidney o 8.Swartz Bay Figure 4.8: Continued. Chapter 4. Results 82 Figure 4.9: Computed elevations at Alberni Inlet for the case of all three segments ruptured together. The model is driven at Barkley Sound. The vertical scale is given in metres, a) Barkley Sound sea. level displacements. 7-| nn 1 1 2 0 4.0 1—6.0 J~ 8.0 T(HR): ALL SEGMENTS RUPTURE Figure 4.9: b) Inlet head sea level displacements. Chapter 4. Results 84 Figure 4.10: Contours of sea-bottom displacement for the case when the Juan de Fuca segment is ruptured alone. Contour intervals are given in centimetres. Chapter 4. Results 85 Figure 4.11: Time-series of computed elevations in the deep-ocean model for the case of the rupture of the Juan de Fuca segment alone. Station locations are indicated in Fig.4.2. The vertical scale is given in metres. Chapter 4. Results 86 o in F . P o r t Hardy o o o in' E 2 . Q u a t s i n o Sound o o o in' El o in" D.Tofino in. i 0.0 -1 2.0 1 4.0 1 6.0 1 8.0 T(HR);JUflN DE FUCR RUPTRE Figure 4.11: Continued. 10.0 12.0 Chapter 4. Results 87 Figure 4.11: Continued. Chapter 4. Results 88 Figure 4.11: Continued. 89 Chapter 4. Results 0.0 1.0 T i l 2.0 3.0 4.0 T(HR):SEGMENTS COMPARISON r 5.0 6.0 Figure 4.12: A comparison of the junction points at the Juan de Fuca Strait for the first two cases of source displacement: the rupture of all three segments (solid lines) and of only the Juan de Fuca segment (dashed lines). Station locations are indicated in Fig.4.2. Chapter 4. Results 0.0 1.0 90 2.0 i 3.0 r 4.0 T(HR){SEGMENTS COMPARISON Figure 4.12: Continued. 5.0 6.0 Chapter 4. Results 91 Figure 4.13: Time-series of computed elevations in the deep-ocean model for the case of the Winona and Explorer segments ruptured together. Station locations are indicated in Fig.4.2. The vertical scale is given in metres. Chapter 4. Results 92 F.Port Hardy o o o E2.Quatsino Sound o TI o o o D.Tofino o o" _ I o.o 1 2.0 1 1 4.0 T(HR); 6.0 WINQNH I 1 8.0 EXPLORER Figure 4.13: Continued. 10.0 12.0 Chapter 4. Results 93 o (NJ J-Cape St. James O O o I.Queen Charlotte C i t y CM o o CM H.Rennel Sound o o" G.Kelsey Bay o CM . o.o -1 2.0 1 1 1 4.0 6.0 8.0 T(HR): WINQNfi I EXPLORER Figure 4.13: Continued. 10.0 12.0 Chapter 4. Results 94 o Ml fM ' ID O rsi. M2 1 0.0 0.5 1.0 1 1 1 T 1.5 2.0 2.5 3.0 T(HR):MAX D1SP.:WIN0NR I EXPLR Figure 4.13: Continued. 3.5 Chapter 4. Results 95 8. Swart z Bay LV": 7.Sidney 4.Victoria 1. J.De Fuca junction o.o -I 2.0 1 4.0 T(HR):GF; 1 6.0 WINDNfl I 1 8.0 10.0 12.0 EXPLOR Figure 4.14: Time-series of computed elevations in the fine grid model for the case of the W i n o n a and Explorer segments ruptured together. Station locations are indicated in Fig.4.6. T h e vertical scale is given in metres. Chapter 4. Results 16.Burrard Inlet 14.Point Atkinson 13. 1 2 . P o i n t Grey 11 .Richmond 10.Delta 9.Nanaimo o.o I 2.0 1 1 1 4.0 6.0 8.0 T(HR);GF: WINONA I EXPLOR Figure 4.14: Continued. 10.0 12. T(HR);GF; WINONH I EXPLOR Figure 4.14: Continued. Chapter 4. Results 98 26.Comox 7-1 0 0 1 1 I I I 2.0 4.0 6.0 8.0 10.0 TIHR);GF: WINONA I EXPLOR Figure 4.14: Continued. l 12.0 Chapter 4. Results 99 Figure 4.15: Computed elevations at Alberni Inlet for the case when the W i n o n a and Explorer segments are ruptured together. The model is driven at Barkley Sound. T h e vertical scale is given i n metres, a) Barkley Sound sea level displacements. Chapter 4. Results o cr> — Figure 4.15: b) Inlet head sea level displacements. Chapter 4. Results 101 Figure 4.16: Three-dimensional views of water surface elevations at various time intervals. The shoreline is represented as a 20-cm step. Maximum wave heights are indicated in metres at each plot. Chapter 4. Results Figure 4.16: Continued. Figure 4.16: Continued. Figure 4.16: Continued. Figure 4.16: Continued. Chapter 4. Results 106 o B.Bamfield • i 0.0 1 30.0 1 i 1 60.0 90.0 T(MIN);LINEAR: 1 120.0 150.0 COMPARE 1 i 180.0 210.0 (a) Figure 4.17: The near linearity of the wave equations is revealed by applying two sources with one's magnitude being double of the other's. Solid lines denote the maximum elevation of 5 m whereas dashed lines correspond to the 2.5 m source. All three segments are ruptured in both cases. Station locations are indicated in Fig.4.2. Chapter 4. Results 107 Chapter 4. Results 108 o Figure 4.18: The near linearity of the wave equations is shown by comparing results of two sources: first for the case when all three segments are ruptured as a whole (solid lines), and in the other case by adding the calculations for each segment rupturing separatelj'. Chapter 4. Results 109 7 2.0 0.0 i 2.0 T 4.0 T(HR); 1 6.0 LINEARITY r 4.0 6.0 T(HR); LINEARITY Figure 4.18: Continued. 12.0 12.0 Chapter 5 Conclusion T h e similarity of the Cascadia subduction zone's tectonic features to those of other subduction areas i n the world, together with results from recent geodetic surveys of the ground motion as well as records of prehistoric large earthquakes i n the area, suggest that a megathrust earthquake may occur in western British Columbia in the near future. T h e compelling evidence for a large seismic event off British Columbia has prompted the study of the local tsunami response. Three numerical hydrodynamical models covering different areas of the studied region were used to calculate tsunami generation and propagation. Computed results include tsunami travel times and sea surface displacements and may provide useful information for the emergency planning response to tsunami hazards. Three cases of source motion were studied. Simulated wave amplitudes have been plotted for major tide gauge locations and some densely populated areas along the British Columbia coast. T h e wave pattern observed at each location depends strongly on local topography, such as the width of the continental shelf, and the shoreline configuration, both of which affect wave refraction and reflection, the amount of shoaling and the excitation of resonance. 110 Chapter 5. Conclusion 111 Results have shown that the most affected area was along the outer coast of Vancouver Island because of its proximity to the source. The maximum positive and negative sea surface displacements in that area are respectively 8.5 and -9.6 m when the zone's three upper segments rupture together in a 5m maximum upthrust; 6.7 and -9.0 m for the breaking of the Juan de Fuca segment alone; and 8.5 and -11.3 m for the rupture of the Winona and Explorer segments together. The travel times of the first positive wave to Victoria and Vancouver, the two main business and commercial centres in British Columbia, are 94 and 187 minutes respectively if all three segments break as a whole and are correspondingly 127 and 226 minutes if only the Winona and Explorer segments rupture together. With the exception of those locations that are very near the source, the largest waves usually come after the leading waves, as is seen in the Strait of Georgia system. The wave amplitudes observed there, although the area is sheltered by Vancouver Island and the wave is attenuated by friction, are still significant. Around the Vancouver area, the maximum wave amplitude is approximately 1 m. Low-lying areas such as Richmond are therefore subject to threats of flooding when the tsunami arrives at high tide. The maximum positive and negative displacements in the Strait of Georgia fine grid model are respectively about 6.8 and -6.4 m and occur near the north-end junction (Johnstone Strait) in all three cases of source motion. At the head of Alberni Inlet, the waves are seen to be drastically magnified through resonance, with magnitudes about three times those in Barkley Sound. The major limitation on the tsunami amplitude computations comes from the uncertainty in specifying the source. A more confident description of source parameters Chapter 5. Conclusion 112 such as the dimension of the displacement zone would enhance the accuracy of the computed results. The wave models, on the other hand, could also be improved in a number of ways: 1 / The area of computation should be extended southwards to include the entire Cascadia zone. The Gorda plate and the lower portion of the Juan de Fuca plate were not considered in this study, which might result in different wave forms if any of these two segments were to rupture. 2/ The present model has omitted a detailed study of most of the inlet systems along the British Columbia coast. Only Alberni Inlet was treated in detail in the inlet computation. The results from the deep-ocean model therefore represent only the wave response adjacent to the mouths of the inlets. Finer resolution models are required to improve the results' accuracy in bays and inlets. 3/ The treatment of Alberni Inlet under the present scheme is, however, not yet complete. The results only provide the wave amplitudes at the shoreline; it does not give the run-up heights. The latter information is considered to be more important becausefloodingby tsunami waves causes most of the damage. The modelling of the run-up process, which may include bore formation, breaking and rushing-up, requires modification of the present assumption of zero mass transport at the coastline and development of high-resolution models to account for complex structures such as small islands and narrow branches at the head of the inlet. A better scheme of coupling the shelf and inlet models may also be required. Chapter 5. Conclusion 113 4/ Open boundary conditions in the deep-ocean model could be improved by developing more stable two-dimensional radiation conditions. Although spurious reflections do not seem to appear significantly in most stations of interest, the computation could be carried out more efficiently and economically without attaching the extra section of the extended grid which is almost the size of the original grid. 5/ Sea level changes and currents due to astronomical tides might have to be taken into consideration. In inlets and low-lying areas such as Richmond, the inclusion of tidal effects could provide a better estimate of the maximum wave heights. Including tidal currents with those of the tsunami wave could also give a better estimate of energy dissipation in near-shore areas. Bibliography Abe, K., 1979, Size of great earthquakes of 1837-1974 inferred from tsunami data, Journal of Geophysical Research, Vol 84, no. B4, pl561-1568. Baines, W.D., 1958, Tidal currents in constricted inlets, Proc. 1958 Coastal Engineering Conference, Miami, p545-561. Camerlengo, A.L., and O'Brien, J.J., 1980, Open boundary conditions in rotating fluids, Journal of Computational Physics, pl2-35. Crean, P.B., Murty, T., and Stonach, J., 1988, Mathematical modelling of Tides and Estuarine Circulation, Lectures notes on coastal and Estuarine studies #30, Springer-Verlag, 471pp. Dragert, H., and Horner, R., 1989, private communication. Dragert, H., and Rogers, G.C., 1988, Could a megathrust earthquake strike southern British Columbia, GEOS, Vol 17, no.3, p5-8. Dunbar, D., LeBlond, P., and Murty, T., 1989, Maximum tsunami amplitudes and associated currents on the coast of British Columbia, Science of Tsunami Harzards, vol 7, no.l, p3-44. Engquist, B., and Majda, A., 1977, Absorbing boundary conditions for the numerical simulation of waves, Mathematics of Computation, vol 31, no.139, pl-24. Heaton, T., and Kanamori, H., 1984, Seismic potential associated with subduction in the northwestern United States, Bulletin of the Seismological Society of America, vol 74, no.3, p933-941. Heaton, T., and Hartzell, S.H., 1987, Earthquake hazards on the Cascadia subduction zone, Science, Vol 236, pl62-168. Henry, R.F., 1982, Automated programming of explicit shallow-water models, Part 1, Canadian Technical Report of Hydrography and Ocean Sciences, No 3. Hwang, L.S., and Divoky, D., 1970, Tsunami generation, Journal of Geophysical Research, vol 75, no 33, p6802-6817. 114 Bibliography 115 [13] Lemon, D.D., LeBlond, P.H., and Osborn, T.R., 1979, Seiche excitation in Port San Juan, British Columbia, Journal of the Fisheries Research Board of Canada, Vol 36, no. 10, 1223-1227. P [14] Loucks, R.H., 1962, Investigation of the 1960 Chilean tsunami on the Pacific coast of Canada, M.Sc. Thesis, University of British Columbia, Vancouver, B.C., 21p. [15] Marshall Macklin Monaghan, 1986, Development management in tsunami hazard areas of Port Alberni, Consultants' report prepared for the City of Port Alberni. [16] Murty, T.S., 1977, Seismic sea waves - Tsunami, Bull. Fish. Research Board of Canada, 198, 337pp. [17] Murty, T.S., and Boilard, L., 1970, The tsunami in Alberni Inlet caused by the Alaska earthquake of March, 1964, pl65-187, in W.M. Adams (ed.) Tsunamis in the Pacific Ocean, East West Centre Press. Honolulu, Hawaii. [18] Murty, T.S., and El-Sabh, M.L, 1985, Numerical simulation of the tsunami due to a predicted large earthquake in the St. Lawrence estuary, Proceedings International Tsunami Symposium 1985, p75-81. [19] Rogers, G.C., 1988, An assessment of the megathrust earthquake potential of the Cascadia subduction zone, Canadian Journal of Earth Sciences, vol 25, no.6, 844-852. P [20] Thompson, S.A., 1982, Trends and developments in global natural disasters, 1947 to 1981, Natural Hazards Research Working paper #45, Boulder: University of Colorado, Institute of Behavioral Science. Appendix A Approximation of a spherical polar grid by a rectangular grid Given a parametric surface S (Fig.A.l) with parameters u and v, the position vector of a point lying on the surface S is denoted as (u,v) = (x(u,v), y(u,v), z(u,v)). The surface area of S can be shown to be (A.l) where R is the domain of integration. For a spherical surface with parameters 6 and 4> (Fig.A.2), the position vector of a point on the sphere is given by: r(c/>,6) = (psinc/>cos#, psin</>sin#, pcosc/>). The surface area of a spherical segment S, bounded by fa < (f> < fa + Ac^> and 6 < 6 < 6 + A#, 1 1 can then be determined from Eq(A.l) as follows: (A.2) It is easy to show that Eq(A.2) reduces to (A.3) where R is bounded by fa < <p < fa + A<£ and 0 < 6 < 6-^ + A0. X 116 Appendix A. Approximation of a spherical polar grid by a rectangular grid 117 Figure A.l: A parametric surface S The integral in Eq(A.3) could be evaluated by using the Intermediate Value Theorem, and it is found to be area(S) = p sin d> A(pA6 2 (A.4) A for some d> in (<j>i,<f>i + Ad>). Alternatively, Eq(A.4) could be derived directly from simple geometry as is shown in Fig. A.3. For small Ad> and A9, it is obvious from Fig. A.3 that the-surface area of S could be approximated by: area(S) = p sin(f>A<f)A6 2 (A.5) which is equivalent to Eq(A.4). In terms of the latitude £, where £ = 90° — d>, both Eqs(A.4) and (A.5) could be Appendix A. Approximation of a spherical polar grid by a rectangular grid Figure A.3: A spherical segment S of sides p&<f> and psin^Af?. 118 Figure A.4: A spherical surface approximated by a rectangular grid. The numerals indicate the latitude-longitude range of the coarse grid used in the deep-ocean model. rewritten as area(S) = p cos(A<f>A6 2 (A.6) Consider now two small spherical cells Si and S2 (Fig.A.4), both of sides A<f> and A9, located respectively at the centre and the northeastern corner of the spherical polar grid. From Eq(A.6), the areas of SI and S2 are area(Sl) . = p cos 51.5° A<j>A6 (A.7) area(S2) = p cos 56.9°<A<j>A6 (A.8) 2 2 Because the rectangular grid is tangential to Earth at point C, the small curvilinear cell Si centred at C could be closely approximated by its projected image onto the Appendix A. Approximation of a spherical polar grid by a rectangular grid 120 rectangular grid. In other words, as A<p and A # are small enough, Si could be treated as if it were a rectangular element. The surface area given in Eq(A.7) is then viewed as the surface area of a rectangular grid cell. The distortion in surface area of representing a curved surface, such as the one in the northern corner, by a rectangular grid element is thus given by the ratio of Eqs(A.7) and (A.8): area(S2) _ cos56.9° area(Sl) ~ co 51.5° ~ ° ( 9 ) S The distortion is about 12%. Along the southern boundary, the corresponding ratio is: cos46.4° ,SSTi= = 1 1 1 <'> A 10 with the distortion about 11%. Hence, the error of using a rectangular grid in computing surface areas over the specified latitude range is about 12%. The surface area distortion would subsequently cause error in computing the water surface elevation 77, which is determined next. We have energy flux = \j gh.area.E = constant (A.ll) where E is the energy density. Or Igh.area.n as rj is proportional to energy. 2 2 = constant (A. 12) Appendix A. Approximation of a spherical polar grid by a rectangular grid 121 Since the energy density is not the same for those grid elements located at the centre and along the boundaries (for instance, SI and S2), 77 has to be adjusted so as to conserve energy within each grid element. Hence in the northern boundary, the 77 calculated in the rectangular grid would be smaller than that in the spherical grid. On the other hand, 77 would be larger in the southern boundary. To calculate the percentage of difference, we have, from Eq(A.12): with k the proportional constant. By including errors in A and 77, where eA77 and e A j 4 indicate the perturbations in surface elevation and surface area respectively. Simplifying, i f (1 + 2e-^ + 0(e )) 2 = j(l - e—) (A.17) To the first order of e, AA \2r,A \ = \k— I v 1^1 - ^ 1 (A.18) ( - » Appendix A. Approximation of a spherical polar grid by a rectangular grid Since rj - A- ' 122 2 , AT/, .AA. A 2—- = — H v ( - ) A 20 thus, . An, ' n 1 1. AA 2 A (A.21) Hence in the northern region, the computed sea surface elevation in a rectangular grid could be under-estimated by at most 6%. Correspondingly, the amplitude is over-estimated by 6% in the south. Appendix B Constricted flow in Burrard Inlet Burrard Inlet which is situated on the southern mainland coast of British Columbia has three deep basins ( F i g . B . l ) . W i t h the exception of English Bay, the outer basin which is open to the Strait of Georgia without an entrance channel, the two inner basins are connected by narrow, shallow channels. Because of their constricted entrances, the inner basins behave like Helmholtz resonators. A tidal model i n Burrard Inlet presented by Baines (1958) estimated the flow i n the area under the interaction of the two inner basins. T h e model solutions consist of the flow velocities in the two narrows: Ci(C T/ a ~ (C - 1)(C 3 7 7 - 1) cosa - 1)(1 - tan a tan a ) - C,C 3 * cos a 5 B 2 X 3 . 2jd ^ T l " 2 (B.l) V b = -TF> run ivi V 1 \ 7TFT ( 6 — 1 ) ( C — 1)(1 — t a n a t a n a ) — C C 3 7 3 5 4 5 x c o s a 2 x sin(— - a ) 1 2 (B.2) where the subscripts a and b refer to the First and Second Narrows respectively, and all of the constants C and a are defined as follows: 'a 123 Appendix B. Constrictedflowin Burrard Inlet 124 BURRARD INLET WITH APPROXIMATE DEPTHS Figure B.l: Burrard Inlet showing the three basins. (After Baines (1958)) 1V 1 fgT C = 167ri2 2 (BA) a Q 9T a ±* L A 2 (B.5) a 2 a 2 gT a 47T L A 2 C7 = 4 (B.6) b 2 0 2 gT a 4ir*L A 2 (B.7) a b C = 2 ] 14 ] fT b 6 4ir LbA 2 2 A 3 J (B.8) (B.9) Appendix B. Constricted How in Burrard Inlet 125 C = /PV (B.10) 8 where T is the period of tide, A and A3 are surface areas of the second and third 2 basins respectively, L is the channel length, a is the channel cross sectional area, R is the channel hydraulic radius, with the subscripts a and b referring to the First and Second Narrows as stated earlier. The as are, in turns, defined as follows: C4C5 tan c t = (tana + tana )/[l - tana tana3 - — -— —] (C — 1)(G7 — 1J 2 5 3 5 (B.ll) 3 tana = ^ ^ - (B.12) 3 tan a = 4 c fi (B.13) C« — 1 tana = ^ - y (B.14) 6 More details can be found in Baines (1958). For both basins, resonance occurs when the denominators of Eqs(B.l) and (B.2) go to zero, that is, ( C - 1)(C - 1)(1 -tan a tana ) - C C = 0 s 7 3 5 4 5 (B.15) Using the Inlet dimension values suggested by Baines (1958), together with the flow velocities V = V =3.1m/s, the period T could be easily solved from Eq(B.15) a b by direct substitution into Eqs(B.4) to (B.14). T was found to be about 4.1 hours.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Assessment of tsunami hazards on the British Columbia...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Assessment of tsunami hazards on the British Columbia coast due to a local megathrust subduction earthquake Ng, Max Kin-Fat 1990
pdf
Page Metadata
Item Metadata
Title | Assessment of tsunami hazards on the British Columbia coast due to a local megathrust subduction earthquake |
Creator |
Ng, Max Kin-Fat |
Publisher | University of British Columbia |
Date Issued | 1990 |
Description | Strong evidence suggests that the Cascadia subduction zone, off the west coast of Canada and the United States, is strongly seismically-coupled and that a possible megathrust earthquake might occur in that area in the near future. A study of tsunami hazards along the Canadian west coast associated with such a hypothetical earthquake is presented in this report. Numerical simulations of tsunami generation and propagation have been carried out using three models based on shallow water wave theory. Three cases of ground motion representing the ruptures of different crustal segments in the area have been examined. Computed results provide information on tsunami arrival times and a general view of the wave height distribution. The outer coast of Vancouver Island was found to be the most strongly affected area. At the head of Alberni Inlet, wave amplitudes reached up to three times the source magnitude. Inside the Strait of Georgia, the wave heights are significant enough to receive closer attention, especially in low-lying areas. |
Subject |
Tsunamis -- British Columbia -- Mathematical models Earthquakes -- British Columbia -- Mathematical models Plate tectonics -- British Columbia -- Mathematical models |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-10-28 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085017 |
URI | http://hdl.handle.net/2429/29633 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1990_A6_7 N42.pdf [ 5.09MB ]
- Metadata
- JSON: 831-1.0085017.json
- JSON-LD: 831-1.0085017-ld.json
- RDF/XML (Pretty): 831-1.0085017-rdf.xml
- RDF/JSON: 831-1.0085017-rdf.json
- Turtle: 831-1.0085017-turtle.txt
- N-Triples: 831-1.0085017-rdf-ntriples.txt
- Original Record: 831-1.0085017-source.json
- Full Text
- 831-1.0085017-fulltext.txt
- Citation
- 831-1.0085017.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085017/manifest