To avoid system errors, if Chrome is your preferred browser, please update to the latest version of Chrome (81 or higher) or use an alternative browser.
Click here to login if you're an NAE Member
Recover Your Account Information
Author: Philip L.-F. Liu
Computational models of tsunamis can be used to design and operate early warning systems. The Earthquake At 00:58:53UTC on December 26, 2004, a strong earthquake struck southwest of Banda Aceh in northern Sumatra (3.30?N, 95.78?E). The earthquake, which occurred near the junction of the Indian, Australian, and Burma plates, was caused by the sudden release of tectonic strain built up by the subduction of the Indian plate under the Burma microplate (at roughly 60 millimeters [mm]/year) along a fault line stretching from Indonesia in the south to the Andaman Islands in the north. The original earthquake and aftershocks indicate that approximately 1,200 kilometers (km) of this fault ruptured along the northern Sunda Trench. The initial estimate of the magnitude of the earthquake was MW = 9.0 based on the Harvard University centroid moment tensor (CMT) solution. The seismic moment released on the fault plane (strike = 320? and dip =11?) was estimated at 3.57 5 1029 dyne-cm (Ji, 2005). The duration of the associated total rupture was 200 seconds (sec), and the peak slip was about 20 meters (m). The rupture propagated northwestward for nearly 400 km at a speed of 2.0 km/sec. In the seismic body-wave analysis, only the first 220 sec of seismic data were used to constrain the slip. Hence, the body-wave analysis cannot resolve the later, smaller slip if it occurred further north. Analyzing the normal modes excited by this earthquake, Stein and Okal (2005) estimate that the seismic moment was as large as 1.0 5 1030 dyne-cm (MW = 9.3). They also suggest that the larger moment can be accounted for by 11 m of slip on a fault 1,200 km long and 200 km wide (down-slip dimension). This would indicate that the slow slip probably occurred over the northern part of the 1,200-km long rupture zone indicated by the aftershocks. It has also been suggested that the rupture propagated to the north relatively slowly, producing a long source duration of approximately 30 minutes (Dr. Kenji Satake, Active Fault Research Center, National Institute of Advanced Industrial Science and Technology, Japan, personal communication). However, the entire rupture process, including the cause of the suggested slow rupture propagation, remains unresolved. The Tsunami The earthquake triggered giant tsunami waves that propagated throughout the Indian Ocean, causing extreme inundation and destruction along the northern and western coasts of Sumatra. Within hours, the tsunami devastated the distant shores of Thailand to the east and Sri Lanka, India, and Maldives to the west. The tsunami also caused destruction in Somalia and other nations of east Africa. The highest recorded tsunami run-up occurred on the west coast of Sumatra, where the height of the water reached almost 35 m. In Indonesia alone, more than 230,000 are dead or missing. The death toll in Sri Lanka and India is more than 50,000. More than 1,500,000 people are believed to be homeless as a result of the tsunami. Post-Tsunami Survey In response to the Indian Ocean tsunami, the international community quickly initiated and organized post-tsunami survey plans. Scientists and engineers from Japan, Korea, New Zealand, and the United States participated in this effort. Within three days of the earthquake, three U.S. survey teams were organized and dispatched to Indonesia, India, and Sri Lanka and Maldives. Later, a small team was sent to east Africa. Most of the team members from the United States were supported by the Earthquake Engineering Research Institute (EERI), National Science Foundation (NSF), and U.S. Geological Survey (USGS). The U.S. teams included researchers from Cornell University, Georgia Institute of Technology, Oregon State University, Portland State University, Texas A&M University, University of Arizona, University of Southern California, and University of Washington; they were joined by researchers from USGS and local scientists from India, Indonesia, and Sri Lanka. The primary objective of the post-tsunami surveys was to collect physical evidence and eyewitness accounts on maximum tsunami heights, maximum run-up heights, inundation distances, areas of inundation, and tsunami deposits. These data are essential for validating numerical models for tsunami-wave propagation over the ocean and run-up/inundation in coastal areas and for improving our understanding of the fault-plane mechanism. Maximum tsunami height is defined relative to sea level and is measured by watermarks on buildings, scars on trees, and rafted debris. Inundation distance is the distance from the shoreline to the inland limit of tsunami flooding; maximum run-up height is the elevation at the inundation distance. Every mark used for these measurements was photographed and its location identified by GPS. Wave heights were plotted for Sri Lanka, the east coast of India, and the north coast of Aceh, Indonesia. Soil samples from tsunami deposits were also collected. The tsunami in India and Sri Lanka deposited sand from the beach and ocean floor in buildings, on top of boulders, and on the ground. Tsunami sand deposits were found at every visited site. The width of the tsunami deposit varied, depending on the characteristics of the deposit and the wave height. The thickness and grain size distribution of tsunami deposits will be correlated to the wave height and flow velocity associated with the tsunami to provide a scientific basis for analyzing paleo-tsunami deposits. These data will also be added to the database for estimating recurrences of tsunamis at given sites. Damage to houses and buildings was also recorded and documented. In general, most wood-framed struc-tures and straw houses were totally destroyed. Some masonry buildings survived, as did almost all reinforced concrete houses (Figure 5). One of the prominent effects of the tsunami was extensive scouring around the foundations of buildings (Figure 6). The scouring appeared to be the result of both incoming and returning flow and was significant enough to undermine foundations. More detailed reports on field surveys can be found in EERI Special Earthquake Report (March 2005) and Liu et al. (2005). Numerical Simulations of Tsunami Propagation Numerical simulations of tsunami propagation have been greatly improved in the last 30 years. Several computational models are being used in the National Tsunami Hazard Mitigation Program, sponsored by the National Oceanic and Atmospheric Administration (NOAA), to produce tsunami inundation and evacuation maps for Alaska, California, Hawaii, Oregon, and Washington. These models will also be used to design and operate early warning systems. The computational models include MOST (method of splitting tsunami), developed originally by researchers at the University of Southern California (Titov and Synolakis, 1998); COMCOT (Cornell multi-grid coupled tsunami model), developed at Cornell University (Liu et al., 1994); and TSUNAMI2, developed at Tohoko University in Japan (Imamura, 1996). All of these models require further validation by large-scale laboratory experiments and field data. A series of numerical simulations of tsunami propagation were performed for the December 26, 2004, earthquake using different fault-plane mechanisms to test their validity. Numerical results were compared with data from field surveys and available tidal-gauge records. In this paper, I describe numerical results based on two fault mechanisms, one provided by USGS (Ji, 2005) and one provided by NOAA (Dr. Vasily Titov, University of Washington, Joint Institute for Study of the Atmosphere and Ocean, personal communication). The USGS fault-plane parameters (Table 1) were calculated using the Harvard University CMT solution. The estimated length of the fault plane was 500 km; the estimated width was 150 km. Because the strike angle was 320?, which denotes the angle from the north to the fault line in a clockwise direction, the fault line pointed northwest. For the seismic moment 3.57 5 1029 dyne-cm, the rigidity of Earth’s crust can be estimated as ? = 4.3 5 1010 N/m2. The NOAA fault-plane parameters (Table 2) are based on suggestions by Stein and Okal (2005) that fault slips occurred over the northern part of the rupture zone. Three subfault planes are oriented along the northern Sunda Trench. The total length of the fault zone is 1,170 km; the width of the fault plane is kept as a constant at 150 km. Using ? = 4.3 5 1010 N/m2 for rigidity of Earth’s crust, the total seismic moment is estimated as 1.13 5 1030 dyne-cm, which is significantly bigger than the estimate with the USGS fault-plane mechanism. The seafloor displacements for both fault-plane mechanisms were calculated using the linear elastic-dislocation theory with a rectangular uniform fault plane (e.g., Mansinha and Smylie, 1971). Because the water is incompressible and the rupture is assumed to occur instantaneously, the water-surface deformation mimics the seafloor deformation. Figures 7 and 8, generated by the USGS and NOAA fault-plane mechanisms, respectively, show the contours of the initial free-surface elevations. The free-surface profiles exhibit N-wave characteristics, with an elevated waveform on the western side and a depression waveform on the eastern side of the fault plane. The maximum positive amplitude of initial free surface based on the NOAA fault-plane mechanism is about 7 m; the maximum based on the USGS fault-plane mechanism is about 4 m. The numerical model COMCOT was used to simulate tsunami propagation. COMCOT adopts a modified leapfrog finite-difference scheme to solve (both linear and nonlinear) shallow-water equations (Liu et al., 1994). (Only the linear shallow-water equations in spherical coordinates are shown here). In the model, the coastline is placed at a water depth equal to or less than 5 m; the model did not take into account the nearshore bathymetry or inland topography, which might change the direction of wave propagation and the overland flow pattern. The simulated domain covers almost the entire Indian Ocean, ranging from 30?E to 110?E longitude and –25?S to 23?N latitude, with a grid size of 2 minutes (ETOPO2). The dimension of grids is 2,401 by 1,441. The simulations were done on a desktop computer (1.5GB RAM and Athlon XP 2600+ CPU). Numerical Results According to the initial free-surface profiles calculated, the NOAA fault-plane mechanism is more energetic and would generate a much larger tsunami. Furthermore, because the NOAA fault plane stretched farther north, the generated tsunami would have much more severe impacts on the coast of Sri Lanka and the east coast of India. The calculated maximum tsunami wave amplitudes along shorelines of Sri Lanka, India, and Aceh, Indonesia, based on the NOAA fault-plane mechanism, are shown in Figures 2, 3, and 4, respectively (see PDF version). The overall agreement between the observation data and the numerical results is reasonable, considering that the present model does not resolve the nearshore bathymetry and topography. The numerical results based on the USGS fault-plane mechanism (not presented here) significantly underestimated wave heights everywhere. The numerical simulations showed that a large portion of wave energy was reflected from Maldives and propagated back towards the south and west coasts of Sri Lanka. This was confirmed by eyewitness accounts near Colombo, Sri Lanka. The University of Hawaii Sea Level Center (supported by NOAA) maintains a network of sea-level stations in the Indian Ocean that report real-time sea-level data on Global Sea Level Observing System (GLOSS) sites. The numerical results of the NOAA model were compared with the available tide-gauge data (with tidal fluctuations removed) at Colombo, Sri Lanka; Gan, Maldives; and Male, Maldives. Figure 9 shows typical results from these comparisons. The numerical results show the correct phase of leading waves, but the amplitudes are lower than the measured data. Moreover, the arrival time of leading waves is about 10 to 20 minutes earlier in the model. However, considering the low grid resolution (2 minutes) and uncertainty in the estimation of fault-plane parameters, the numerical results match the tide-gauge records fairly well. Conclusion Although a very coarse grid size was adopted for the numerical simulations, numerical results, based on the NOAA fault-plane mechanism, match the field survey data fairly well along the coastlines of Sri Lanka, India, and Indonesia. However, there are still discrepancies for both arrival time and wave amplitude between numerical results and tide-gauge data. The discrepancies could be caused by the large grid size used and uncertainties in the fault-plane mechanism. Keep in mind also that the effects of slow rupture propagation have not been considered in the present results. An examination of the data obtained by satellite radar altimeters, including the TOPEX/Poseidon satellite that passed across the Indian Ocean from southwest to northeast about two hours after the earthquake (TOPEX/Poseidon passed the equator at 03:01:57UTC on December 26, 2004) and the Jason-1 satellite that passed by 1:56 hours after the earthquake (Jason-1 passed the equator at 02:55:24UTC on December 26, 2004), can further our understanding of the characteristics of Indian Ocean tsunamis and the fault-plane mechanism. The accuracy of these sea surface elevation data is better than 5 cm (JPL, 2005). These data could be used to improve the parameterization of fault-plane mechanism. Finally, evaluating the impacts of tsunamis in affected coastal regions will require higher resolution bathymetry/topographic data. Nonlinear effects, including wave breaking, must also be clarified and included in the simulation model. Acknowledgment Work associated with this paper has been supported by grants from the National Science Foundation to Cornell University. Xiaoming Wang performed the numerical simulations used in this paper. References EERI (Earthquake Engineering Research nstitute). 2005. Distribution of the Tsunami Heights of the 2004 Sumatra Tsunami in Banda Aceh measured by the Tsunami Survey Team. Available online at: http://www.eri.u-tokyo.ac.jp/namegaya/sumatera/surveylog/ eindex.htm. EERI Special Earthquake Report. 2005. The Great Sumatra Earthquake and Indian Ocean Tsunami of December 26, 2004. Oakland, Calif.: Earthquake Engineering Research Institute. Imamura, F. 1996. Review of tsunami simulation with a finite difference method. Pp. 25–42 in Long Wave Runup Models, edited by H. Yeh, P. Liu, and C. Synolakis. Hackensack, N.J.: World Scientific Publishing Co. Ji, C. 2005. Preliminary Rupture Model. Available online at: http://neic.usgs.gov/neis/eq_depot/2004/eq_041226/neic_ slav_ff.html. JPL (Jet Propulsion Laboratory). 2005. Ocean Surface Topography from Space. Missions—Jason-1. Available online at: http://sealevel.jpl.nasa.gov/mission/jason-1.html. Liu, P.L.-F., Y.-S. Cho, S.B. Yoon, and S.N. Seo. 1994. Numerical simulations of the 1960 Chilean tsunami propagation and inundation at Hilo, Hawaii. Pp. 99–115 in Recent Developments in Tsunami Research, edited by M.I. El-Sabh. Dordrecht, Netherlands: Kluwer Academic Publishers. Liu, P.L.-F., Y-S. Cho, M.J. Briggs, U. Kanoglu, and C.E. Synolakis. 1995. Runup of solitary waves on a circular island. Journal of Fluid Mechanics 302: 259–285. Liu, P.L.-F., P. Lynett, J. Fernando, B.E. Jaffe, H. Fritz, B. Higman, R. Morton, J. Goff, and C. Synolakis. 2005. Observations by the international tsunami survey team in Sri Lanka. Science (in press). Mansinha, L., and D.E. Smylie. 1971. The displacement fields of inclined faults. Bulletin of the Seismological Society of America 61(1): 1433–1440. Stein, S., and E.A. Okal. 2005. Speed and size of the Sumatra earthquake. Nature 31(434): 581–582. Also available online at: http://www.earth.northwestern.edu/people/seth/research/ sumatra2.html. Titov, V.V., and C.E. Synolakis. 1998. Numerical modeling of tidal wave runup. Journal of Waterway, Port, Coastal and Ocean Engineering 124(4): 157–171. TABLES TABLE 1 Fault Parameters from USGS