Natural Hazards and Earth System Sciences Modeling propagation and inundation of the 11 March 2011 Tohoku tsunami

On 11 March 2011 the Tohoku tsunami devastated the east coast of Japan, claiming thousands of casualties and destroying coastal settlements and infrastructure. In this paper tsunami generation, propagation, and inundation are modeled to hindcast the event. Earthquake source models with heterogeneous slips are developed in order to match tsunami observations, including a best fit initial sea surface elevation with water levels up to 8 m. Tsunami simulations were compared to buoys in the Pacific, showing good agreement. In the far field the frequency dispersion provided a significant reduction even for the leading wave. Furthermore, inundation simulations were performed for ten different study areas. The results compared well with run-up measurements available and trim lines derived from satellite images, but with some overestimation of the modeled surface elevation in the northern part of the Sanriku coast. For inundation modeling this work aimed at using freely available, medium-resolution data for topography, bottom friction, and bathymetry, which are easily accessible in the framework of a rapid assessment. Although these data come along with some inaccuracies, the results of the tsunami simulations suggest that their use is feasible for application in rapid tsunami hazard assessments. A heterogeneous source model is essential to simulate the observed distribution of the run-up correctly, though.


Introduction
On 11 March 2011 a M w = 9.0 earthquake occurred 130 km east of the Sendai coast, Japan (Lay et al., 2011).Its epicenter was located offshore Minamisanriku, with an extension of about 400 km along its strike direction and a maximum slip reported by some to exceed 60 m, being 15-20 m on average (Lay et al., 2011;Ozawa et al., 2011).The earthquake triggered a tsunami that reached the coastline 30-40 min later with run-up heights up to about 40 m, causing enormous destruction in the prefectures Iwate, Miyagi, and Fukushima (EERI, 2011;Mimura et al., 2011b;Mori et al., 2011).Around 20 000 people lost their lives, more than 90 % of them due to the tsunami, and 400 000 were reported to be homeless (EERI, 2011;Vervaeck and Daniell, 2011).76 000 houses were completely destroyed and damage costs add up to 300 billion US$ (Dunbar et al., 2011;EM-DAT, 2011;Mimura et al., 2011a;Mori et al., 2011).In addition, the tsunami caused extensive damage to the Fukushima Nuclear Power Plant followed by a radiation leak (Mimura et al., 2011a;Brumfiel, 2011).The earthquake and tsunami were also stronger than what some of the world's largest tsunami barriers were designed for (Cyranoski, 2011).
Through the recorded history, the Pacific coastline of northeastern Honshu has faced a number of destructive tsunamis (see e.g.NGDC, 2011).Two of the most damaging events in recent history were the 1896 and 1933 Sanriku events.The 1896 event was reported as a M w = 8.0 "tsunami earthquake" (Tanioka and Satake, 1996), causing more than 27 000 fatalities.The 1933 event (Kanamori, 1971) was a M w = 8.4 outer rise earthquake, but with considerably less casualties than the 1896 event.The wake-up call following the massive 2004 Indian Ocean disaster led to a revision also of past megathrust events, suggesting that megathrust events in excess of M w = 9 should not be ruled out along major subduction zones (Stein and Okal, 2007).On this basis, tsunamis due to megathrust earthquakes should not come as a surprise.However, there are reports suggesting that tsunamigenic earthquakes of magnitudes greatly exceeding 8 were not expected in the area (Geller, 2011;Monastersky, 2011).
A hindcast of the 2011 Tohoku tsunami is important not only to understand this particular tsunami, but also to prepare for possible future events.Rapid tsunami hazard assessments, although often based on coarse data, particularly valuable for risk assessment and mitigation, since they allow for a fast and cost-efficient hazard mapping.Therefore, the aim of this investigation is to simulate the generation, near and far field propagation, as well as the local and regional inundation of the 2011 Tohoku tsunami event to explore the possibility to model and map the tsunami event accurately based on medium-resolution and freely available data and within a short period of time.This study is not supposed to be a detailed local assessment.

Methodology
The earthquake rupture is instantaneous and purely dip-slip.It is noted that the finite duration of the earthquake rupture was not expected to largely influence the simulation results.The vertical seabed surface displacement was superimposed from a range of earthquake segments calculated by the formula of Okada (1985).Variable slip was enabled by allowing different slip values for each segment, varying the slip linearly along the strike direction.The vertical seabed surface displacement was copied to the sea surface and smoothed, applying a low-pass filter based on linear full potential theory (Pedersen, 2001), which is in close agreement with the well-known model of Kajiura (1963).
To compute the tsunami propagation, we applied the Boussinesq model labeled GloBouss (Pedersen and Løvholt, 2008;Løvholt et al., 2008), allowing for wavelength dependent wave speed (frequency dispersion).The bathymetric data were obtained from The General Bathymetry Chart of the Oceans and converted to a resolution of 4'.To simulate tsunami inundation the nonlinear shallow-water wave (NLSW) equations were used together with the Community Model Interface for Tsunami (ComMIT) (Titov andSynolakis, 1995, 1998;Synolakis et al., 2008;Titov et al., 2011).Inundation modeling in ComMIT was based on three, nested, rectangular computational grids (here called A-, B-, and C-grid) with bathymetry and topography as underlying data sets.By using a one way nesting procedure, ComMIT was coupled with GloBouss, interpolating the output from GloBouss over the A-grid boundary at each time step, using the global propagation simulation to drive the local inundation model (Løvholt et al., 2010).
A range of tsunami scenario simulations were conducted for different slip distributions, both in the dip and strike direction.By comparing simulation results with recorded surface elevations offshore and inundation data, we continuously adjusted the slip distribution to provide new scenarios in order to improve the agreement with observations.In the source refinement process, emphasis was given on comparing the simulated surface elevation with recorded surface elevations from Deep-ocean Assessment and Reporting of Tsunamis (DART) buoys, due to the limited quality of the topographic and bathymetric near shore data.A total of 17 scenarios simulations were conducted during a period of about one month.Out of these, four different labeled scenarios A, B, C, and D are presented below, with emphasis on the best fit scenario D.

Source description
Presently, there are a number of available seismic and geodetic inversions of the earthquake rupture (e.g.Ozawa et al., 2011;Lay et al., 2011;Ammon et al., 2011;Saito et al., 2011); see also GEO Geohazards Supersite (2011) for a more extensive compilation.Most of the inversions have in common that the maximum slip is located in the area east offshore from Minamisanriku in the South to Miyako in the North.As an example Ozawa et al. (2011) report a maximum slip of 27 m close to the epicenter just north of the 38 • latitude, whereas Lay et al. (2011) report a maximum slip of 63 m.From the GEO Geohazards Supersite (2011) compilation, we found that the maximum slip varied between 10-40 m, and the range of inverted slip distribution also appears differently.Hence, using such inversions in tsunami simulations would imply largely different water levels.Nevertheless, the earthquake source applied to the tsunami simulations is implemented to roughly comply with the main trends from the slip inversions, but even more so to make the modeled tsunami match the one observed from the offshore wave gauges and the inundation pattern.Scenario slip distributions and resulting initial surface elevations are displayed in Figs. 1 and 2. Whereas scenario A represents one of the simplest cases with a uniform slip, the other scenarios are all heterogeneous.Moment magnitudes range from 8.85 to 8.95 (see also Fig. 1), and maximum initial elevation range from 4 to 13 m.In all cases, the scenarios are surface rupturing, and a shear stiffness of 40 GPa is assumed.Dip angles are 25 • except for scenario A, which has a dip angle of 15 • .The simulation that overall matched the field observations best is described below.The best fit source (scenario D) comprises 6 × 6 segments along the strike and dip direction, respectively.Source D has a maximum slip of 15-20 m located between latitudes 38 • and 40 • and a total width of 150 km.The total length of the modeled source rupture is approximately 400 km.Compared to the co-seismic slip distribution of Ozawa et al. (2011), the slip distribution of scenario D is extending slightly further north.Some of the segments display a slight overlap due to the bending of the source, resulting in small areas of locally increased slip.However, the overlapping areas are mostly displaying short wavelengths and therefore vanish due to the smoothing.The maximum sea surface response is slightly above 8 m. 4 Results from offshore tsunami simulation Scenario simulations were conducted for the main part of the Pacific Ocean.In Fig. 3 the maximum surface elevation for scenario D is depicted.Outside the eastern coast of Japan, the maximum surface elevation is up to 10 m, while the directivity of the source is about WNW-ESE.Comparison between the surface elevations for the different scenario simulations A-D were evaluated at the DART buoys given in Fig. 3.The comparisons are given in Fig. 4 (dispersive solutions).As shown, scenario A provides too small amplitudes compared to the DART gauges.Scenarios B and C both provide relatively good agreement, but particularly scenario C tends to overestimate the surface elevation.In addition, scenario B gives too early of an arrival for a number of locations.Generally, scenario D gives the best agreement, both with respect to the elevation and arrival time, with errors seldom exceeding 10 % for the maximum leading wave.In summary, the numerical results for the best fit scenario D compare very closely with the recorded data, both in shape and height.Between the simulated data and the recorded data from the DART buoys, there is a small time shift for a few locations.The reason for this shift may be an incorrect location of the buoys in the numerical model, inaccuracies in the bathymetric data, or minor errors in location of the source.Hence, for Figs. 4 and 5 the DART data were given a small time shift for easier comparison between the leading waves of the recorded and modeled waveforms.In this manner the leading waves in the DART data and those from the simulations arrive at the gauges at the same time.The time shift (100 to 200 s) corresponds to an EW displacement of the location of the source of only 15-30 km.
For scenario D, the numerical model is applied in both linear non-dispersive (linear shallow water approximation -"LSW") and linear dispersive mode ("disp") (Fig. 5).By comparing the solutions for the different models, the effect of dispersion becomes evident.As earthquake tsunamis are generally long-crested, the effect must accumulate over long propagation distances to be noticeable.In the present case the wavelength is several hundred kilometers (at a depth of 4-6 km), so the dispersive effect is not found for the closest DART buoys (21 413 and 21 415).However, for the DART buoy close to Hawaii (51 407) and especially for the buoy offshore Mexico/Guatemala (43 413), the leading waves are markedly influenced by dispersion, resulting in a smoother and lower leading wave as evident from comparison with the shallow water simulations.Restating the simple rule of thumb proposed by Kajiura (1963) and later by Shuto (1991), we find that the tsunami may be influenced by dispersions after travelling a distance d = 6 h • (λ/ h) 3 • P −3 .
Here, λ is the tsunami wavelength, and h the still water depth.
In the following, we choose P = 4 and h = 5 km as examples.As the wavelength is not easily defined, we exemplify that the necessary travel distance, by assuming initial wavelengths of 100, 150, and 200 km, gives estimated travel distances of 0.6R, 2R, and 4.7R, respectively, R being the Earth's radius.The distance between the Japan Trench and DART 43413 is about 1.8R.Given that the source width is 150 km, it seems that the rule of thumb above expression gives reasonable estimates for the necessary travel time needed for significant effect of dispersion.
The convergence for the simulations of the tsunami propagation stage is assured by sensitivity tests by using grids with resolution of 4', 2', and 1'.For the closest DARTs the deviation of the leading peak, compared to the solution on the finest grid (1'), is 0.5-1.1 and 0.1-0.3% for the resolutions 4' and 2', respectively.

Inundation simulations
Inundation simulations were conducted for ten study areas along the eastern Japan, with all characterized by a different shape of the coastline and varying topography.The study areas are located in the floodplains of Minamisoma, Soma and Sendai, in Ishinomaki city, and along Minamisanriku, Kesennuma, Rikuzentakata, Otsuchi, Miyako and Kuji (Fig. 6).The size of the study areas varies from 298 km 2 to 1969 km 2 .For the inundation simulations the ComMIT model was used with varying resolution for the A-grid (700 m × 500 m), Bgrid (178 m × 140 m) and C grid (90 m × 90 m).GEBCO data with a resolution of 30"(≈910 m) were applied and 90 m resolution data from the Shuttle Radar Topography Mission, SRTM (Jarvis et al., 2008) were used to represent land elevation in the model.It is stressed that the SRTM data set is not corrected for the subsidence due to the earthquake, which is reported to be up to 1.2 m (Mimura et al., 2011a).This gives rise to a slight underestimation of the run-up expected to be of the same order of the subsidence, which is typically an order of magnitude less than the typical mean run-up.Several studies on tsunami inundation modeling have shown the need for high resolution topographic and bathymetric elevation data to simulate local inundation patterns sufficiently well (e.g.Taubenböck et al., 2009).However, high resolution data are not available in many parts of the world and often not suitable for rapid assessments due to the demands on time and effort to handle huge data sets and related computational requirements.For these reasons SRTM is an important data source for many tsunami risk studies.As other elevation models derived from space or airborne, the SRTM data describe the Earth's surface, including vegetation and settlement structures in the height description (Sun et al., 2003;Hofton et al., 2006).In some of the study areas, this required some manual adjustments of the SRTM data to remove significant offsets (i.e.coastal forests), which would influence the results of inundation modeling considerably and lead to false (underestimated) water levels.This was done in a Geographic Information System (GIS) by correcting offset pixels to the surrounding ground level pixels, based on satellite images derived from Google Earth.Since no detailed land use data were available, this correction was done very conservatively and only for obvious green belts and artifacts.General offsets due to housing in urbanized areas were not removed.Despite removing the most prominent forests, there are still high inaccuracies related to the horizontal and the vertical resolution of the elevation data, which have to be kept in mind when analyzing the results.
Different friction coefficients were tested in ComMIT to account for wave attenuation caused by land cover roughness such as forests or human-made structures, which are supposed to influence inundation characteristics (Gayer et al., 2010;Mimura et al., 2011a).From the Manning coefficients n 2 = 0.0009, n 2 = 0.0017, and n 2 = 0.0012, the latter provides the best results when validating the modeled surface elevation against run-up measurements from the field.The time step for modeling with ComMIT was set to 0.5 s in Otsuchi, Rikuzentakata and Kesennuma and to 1 s in all other study areas.Tides were not considered in the analyses.
For scenarios B, C, and D, source sensitivity was investigated for the locations of Sendai, Rikuzentakata, and Ishinomaki by comparing the observed trim lines derived from satellite data by the Tsunami Damage Mapping Team, Association of Japanese Geographers (2011) with simulated trim lines; results are depicted in Fig. 7. Given the quality of the local data, the scenarios compare fairly well with the observations.However, the simulated trim lines for the individual scenarios may not be as easily separated as the corresponding offshore elevations.For Sendai, the simulated trim lines are hard to distinguish and they compare equally well with the observations.For Rikuzentakata scenario B gives best agreement, whereas scenario D provides best agreement in Ishinomaki.Based on the favorable comparisons for scenario D with the DART observations, further analysis is devoted to this scenario.
The results of the ComMIT simulations for scenario D were compared to (a) 442 run-up measurements derived from The 2011 Tohoku Earthquake Tsunami Joint Survey Group ( 2011) and (b) trim lines derived from the Tsunami Damage Mapping Team, Association of Japanese Geographers (2011).The field survey, which is now published in Mori et al. (2011), comprises more than 5300 field measure-ments in its final version.However, during our investigation we could only use the first measurements of this study available online in the early phase after the tsunami.The final measurements differ in the order of a few decimeters from the first version points used here, which is however considered to be negligible due to the overall accuracy of the underlying elevation data in the model.For each study area (Fig. 6), a mean value was created from the measured points in the field and compared to a corresponding mean value of the modeled surface elevation in these points.In doing so, extreme water levels were smoothed for both the measurements and the modeled inundation.Thus, the water levels in the study areas are directly comparable.
Results from the ComMIT simulations show mean maximum water levels of 5.5 m in the Sendai area, whereas mean maximum water levels of approximately 20 m could be observed in the steep valleys of Rikuzentakata and Otsuchi (Fig. 6).The maximum modeled run-up value occurred in Otsuchi with 49.8 m.The maximum modeled inundation distance was approximately 6 km in Otsuchi, Rikuzentakata and Sendai, which agrees with observations provided by Mimura et al. (2011a).In the Kitakami River close to Minamisanriku the water reached 15 km inland.
In total in Minamisoma, Soma, Kesennuma, Miyako and Kuji, good results have been achieved with a deviation of the mean modeled values from the measurements of 5-13 %.In Sendai and Otsuchi, however, there is a significant overestimation (191 %, 157 %) in the modeled maximum water level, which likely results from inaccuracies in the topography/ bathymetry, the disregard of coastal structures, as well as a local source peak southeast off Otuschi (see Fig. 2).A slip consentration was causing a local source peak offshore Otsuchi.However, the overestimation of the maximum water level in Otsuchi suggests that this source peak is artificial.
A comparison with the trim lines (Fig. 8) shows good agreement in almost all study areas, even in those where the point measurements are dissent (e.g. in Rikuzentakata and Sendai).However, a good match with the trim lines in the steep valleys of Rikuzentakata, and Otsuschi was expected, since the topography changes significantly over a short distance and thus directs water flow narrowly through the valley.

Concluding remarks
A rapid assessment study comparing scenario simulations with observations was conducted and the process of source refinement was demonstrated by displaying the results for a limited selection of the scenarios.Results from the numerical simulations show a good match with the field data available, particularly compared to the offshore DART buoys (Figs. 4  and 5).The influence of dispersion is shown to be important for the transoceanic propagation, significantly reducing the leading wave compared to traditional shallow water models.This effect is known to be important for tsunamis of landslide origin (e.g.Lynett et al., 2003;Løvholt et al., 2008), but is often neglected for tsunamis of earthquake origin, despite simulations suggesting distinct dispersive far field effects for the 2004 Indian Ocean tsunami (Glimsdal et al., 2006).There are some uncertainties and sources of errors in the model and the data, which should be considered when evaluating the results, including the earthquake source model, the horizontal and vertical resolution of the topography, the limited consideration of bottom friction and flotsam, as well as the neglecting of coastal structures like walls or revetments, which might have attenuated the wave impact (EERI, 2011;Mori et al., 2011;Mimura et al., 2011a).Despite these uncertainties the study highlights the feasibility of a combined modeling of earthquake seabed displacement, tsunami generation, propagation and inundation in a rapid tsunami hazard assessment based on medium-resolution data.Some information on land cover needs to be available though, in order to obtain and interpret inaccuracies in digital elevation models and to account for the influence of bottom roughness caused by settlement and vegetation.The approach is considered to be useful for modeling tsunami scenarios and impacts in areas potentially exposed to tsunamis and thus support risk mitigation.

Fig. 2 .
Fig. 2. Initial surface elevations for the four different earthquake source scenarios A-D.

Fig. 3 .
Fig. 3. Maximum surface elevation for scenario D and DART locations for entire Pacific Ocean.Contour lines are drawn for every 20 cm.

Fig. 4 .
Fig. 4. Source sensitivity for a range of DART locations exemplified for scenarios A-D.

Fig. 5 .
Fig. 5. Comparisons of the time history of the surface elevation at eight DART buoys for scenario D only."LSW" and "disp" refer to the linear non-dispersive and linear dispersive solutions.

Fig. 6 .
Fig. 6.Comparison of the modeled maximum surface elevation with field measurements published by the 2011 Tohoku Earthquake Tsunami Joint Survey Group.Mean values summarizing point measurements in each study area and the corresponding mean values of the modeled surface elevation in these points are shown.

Fig. 7 .
Fig. 7. Trim lines for scenarios B, C, and D for Sendai, Rikuzentakata, and Ishinomaki compared to the trim lines derived from pre-/posttsunami satellite images by the Tsunami Damage Mapping Team, Association of Japanese Geographers (2011).

Fig. 8 .
Fig. 8.Comparison of the maximum modeled inundation extent with trim lines derived from pre-/post-tsunami satellite images by the Tsunami Damage Mapping Team, Association of Japanese Geographers (2011).Results are shown for the Kitakami River close to Minamisanriku, Soma, Kesennuma, Minamisoma, and Sendai.