Evaluation of Seismic Rupture Models for the 2011 Tohoku-Oki Earthquake Using Tsunami Simulation

Developing a realistic, three-dimensional rupture model of the large offshore earthquake is difficult to accomplish directly through band-limited ground-motion observations. A potential indirect method is using a tsunami simulation to verify the rupture model in reverse because the initial conditions of the associated tsunamis are caused by a coseismic seafloor displacement correlating to the rupture pattern along the main faulting. In this study, five well-developed rupture models for the 2011 Tohoku-Oki earthquake were adopted to evaluate differences in simulated tsunamis and various rupture asperities. The leading wave of the simulated tsunamis triggered by the seafloor displacement in Yamazaki et al. (2011) model resulted in the smallest root-mean-squared difference (~0.082 m on average) from the records of the eight DART (Deep-ocean Assessment and Reporting of Tsunamis) stations. This indicates that the main seismic rupture during the 2011 Tohoku earthquake should occur in a large shallow slip in a narrow range adjacent to the Japan trench. This study also quantified the influences of ocean stratification and tides which are normally overlooked in tsunami simulations. The discrepancy between the simulations with and without stratification was less than 5% of the first peak wave height at the eight DART stations. The simulations, run with and without the presence of tides, resulted in a ~1% discrepancy in the height of the leading wave. Because simulations accounting for tides and stratification are time-consuming and their influences are negligible, particularly in the first tsunami wave, the two factors can be ignored in a tsunami prediction for practical purposes.


InTROdUcTIOn
On 11 March 2011, the M w 9.0 Tohoku-Oki earthquake, with an epicenter close to the conjoint site of four tectonic plates (Fig. 1), caused considerable ground tremors and huge tsunami waves in and around northern Honshu, Japan.Based on over 100 years of seismic observations, the northern Honshu area is considered to be a relatively calm area compared to the active seismogenic zone of southern Honshu (e.g., Ando 1975).Hence, this catastrophic event, which caused extensive and severe structural damage in northeastern Japan, challenges the current understanding of seismotectonics for plate convergence around Japan.The findings from the Global Centroid Moment Tensor Solutions indicate that the mainshock of the 2011 Tohoku earthquakes was characterized as a reverse fault with a compression axis in the WNW-ESE direction (Nettles et al. 2011).Most of the aftershocks occurred in areas of persistent background seismicity, but were clearly subject to the Pacific Plate subducts under the Okhotsk Plate and were distributed in an approximately 200-km wide and 500-km long area (cf. the seismic catalogue of Japan Meteorological Agency -JMA website http://www.jma.go.jp).
To study the deformation processes of this earthquake, many rupture models of the mainshock have already been created using teleseismic and regional seismometer records, GPS observations, and tsunami wave height records (e.g., Ammon et al. 2011;Fujii et al. 2011;Ide et al. 2011;Lay et al. 2011a, b;Lee et al. 2011;Maeda et al. 2011;Ozawa et al. 2011;Simons et al. 2011;Tsushima et al. 2011;Yamazaki et al. 2011;Yue and Lay 2011;Zhao et al. 2011).All of these studies explain the faulting of the mainshock in different frequency bands extracted from the great complexity of the 2011 Tohoku earthquake.A realistic, three-dimensional coseismic rupture model, in fact, is difficult to create directly through band-limited ground-motion observations.A potential indirect method for improved evaluation is to conduct a tsunami simulation to verify the rupture model reversely.It is known that tsunamis are of relatively long wave length and insensitive to the rising time of an earthquake rupture.The initial sea surface deformation is caused by a coseismic seafloor deformation and is highly correlated to the rupture pattern of the earthquake.Because the initial condition dominates the tsunami waveforms and their propagation in time and space, an accurate simulation should presumably be able to examine the focal rupture pattern obtained from various seismic approaches.The primary objective of this study is aimed at differentiating the impacts from a range of mainshock rupture models in tsunami observations.
In addition to the initial conditions, there are other factors that may have different influences on the accuracy of simulated tsunamis.Sensitivity studies using numerical models help quantify these effects.For example, Løvholt et al. (2012) conducted a numerical study using the Boussinesq model (Løvholt et al. 2008;Pedersen and Løvholt 2008) and the Community Model Interface for tsunamis (e.g., Titov andSynolakis 1995, 1998) to examine the propagation and inundation of the 11 March 2011 Tohoku tsunami, and evaluated the sensitivity of the simulated tsunamis to define the horizontal grid resolution of the model.The authors found that the grid resolution of 4 arc-min and 2 arc-min had approximately 0.5 -1.1% and 0.1 -0.3% deviations, respectively, from the leading peak of the 2011 Tohoku tsunami simulated by the 1 arc-min resolution model (Løvholt et al. 2012).They also studied the dispersive effect of the tsunami wave propagation and found that it may not be included in the data recorded by DART buoys for the 2011 Tohoku tsunami event.
As a supplement to the study of Løvholt et al. (2012), the present work selects five distinct rupture models to discuss the rupture distributions in the 2011 Tohoku earthquake and to discriminate between critical factors.These models, denoted as Model I, II, III, IV, and V, are from Ammon et al. (2011), Fujii et al. (2011), Lee et al. (2011), Maeda et al. (2011), andYamazaki et al. (2011) respectively.The five seismic rupture models considered in this study are summarized in Table 1 and will be described in the following section.A three-dimensional numerical model revised from the Princeton Ocean Model (Blumberg and Mellor 1987) was adopted to simulate the tsunamis induced by this earthquake.The wave height and timing of the simulated tsunamis were compared with the observations at eight DART stations in the Western Pacific.
This study also quantified the influence of ocean stratification and tides on the tsunami propagation.These two oceanic conditions are normally disregarded in simulation.This practice may be acceptable because the influence of stratification may only become significant at abrupt topographies due to the energy transfer from surface waves to the generation of internal waves (e.g., Jan et al. 2008).Tsunamis should be less affected by ocean stratification in deep waters even in regions of abrupt topography.This is also because the time and length scales of tides are dramatically larger than that of tsunamis, and thus the interactions between tides and tsunamis can hardly affect the characteristics of tsunamis.Nevertheless, their influences should be quantified by, e.g., percentage of fitness, instead of a qualitative description.
The tsunamis triggered by the 2011 Tohoku earthquake were mostly observed through ocean bottom pressure recordings by DART stations located in the Pacific Ocean and around the Pacific Rim.The DART system is a component of the larger US National Tsunami Hazard Mitigation Program, which is designed to directly measure the tsunami height toward a far-field community (Milburn et al. 1996; see more information at DART website, http://www.ndbc.noaa.gov/dart.shtml).Generally, the standard sampling interval of a DART time series is 1-minute.Once a seismic event is detected, the sampling interval is switched to 15-second intervals so that the tsunami waveform can be recorded in a high temporal resolution.In this study, we collected the sea-surface height time series from eight DART stations in the Western Pacific (red squares in Fig. 1), which comprise near-field and far-field observations.After removing the tidal amplitudes by using the least-square fitting with well-known tidal periods (here we consider the diurnal, semi-diurnal tidal constituents and the seasonal tides), the sea-surface data from the eight DART stations clearly revealed the seismic Rayleigh waves as well as tsunami waveforms (right panel in Fig. 1).The speed of Rayleigh waves can be determined from the seismic-phase arrival time versus the spreading distance from the eight DART records, which is roughly 3.62 km s -1 .Tsunami waves propagate much slower than seismic waves and were simulated in a three-dimensional, quasi-realistic ocean model described in the next section.

SEISMIc RUpTURE MOdELS And InITIAL cOndITIOnS In TSUnAMI SIMULATIOn
As mentioned above, the focal mechanism and slip history of the 2011 Tohoku earthquake have been explained by numerous rupture models.The hypocenter of the mainshock was established at 38.1035°N, 142.861°E and 24 km for the focal depth (JMA catalogue).The major rupture of the mainshock is considered to have taken place at a fault plane approximately 200 ~ 300 km in length and 120 ~ 200 km wide (e.g., Ammon et al. 2011;Lay et al. 2011a;Simons et al. 2011;Yamazaki et al. 2011).Yamazaki et al. (2011) found that slight increases in hypocenter depth and average fault dip have a negligible effect on near-field tsunami propagation.The primary differences in simulations are typically caused by differences in the various features among the rupture models.Five rupture models were selected in this study to represent the typical determinations of focal character based on various combinations of seismological and geodetic observations.The five models are considered as Model I: the main rupture area is centered at the fault plane and elongating along the strike of the fault plane (Ammon et al. 2011), Model II: the main slip takes place around the mainshock's hypocenter, ranging from a very shallow depth down to 40 km along the fault dip (Lee et al. 2011), Model III: the largest slip is described at the shallow depth adjacent to the trench (Yamazaki et al. 2011), and Models IV and V are two simplified rupture models both showing a large slip at the trench area (Fujii et al. 2011;Maeda et al. 2011).The detailed focal mechanisms and fault plane parameters of Models I thru V are listed in Table 1.Among these models, Model I was derived from the teleseismic observations and a few high-rate GPS records.Model II was developed by a parallel inversion technique which combined the teleseismic, local strong motion and near-field coseismic geodetic data.Models IV and V used the data from Japanese sea-bottom pressure gauges (TM1 and TM2). Figure 2 shows the initial sea-surface conditions derived from the five rupture models for the tsunami   For the simulation in the stratified ocean, 25-yearmean temperature and salinity (T and S) profiles obtained from historical CTD data managed by the Ocean Data Bank of Taiwan's National Science Council were used as the initial fields.Figure 3 shows the mean T and S profiles over the water depth from 0 to 4000 m.For depths below 4000 m,  1).Tsunamis are generated by seafloor deformations caused by earthquakes.Here, the initial conditions of tsunami wave propagation were converted from the seismic rupture models by using the Okada's dislocation model (Okada 1985).By comparison, Models II, III, and V had larger maximal initial sea-surface height, which was up to 14, 16, and 14 m, respectively.The maximal initial sea-surface height in both Models I and IV was merely ~6 m (Fig. 2).The initial sea surface conditions contributed from the two aftershocks were relatively small, which may not generate significant tsunami waves in this event.
The far-field tsunami waveform may be modified by seafloor geometry, bottom friction, tides, and the baroclinic effect from the stratified ocean.The influence of the latter two factors on the tsunami simulation is discussed in this study.

MOdELIng OF TSUnAMI pROpAgATIOn
The three-dimensional, hydrostatic, primitive equation Princeton Ocean Model (POM) delineated in Blumberg and Mellor (1987) was revised and adopted to simulate tsunamis.For brevity, the governing equations are not demonstrated here.The model uses a split time step separating the two-dimensional external mode from the three-dimensional internal mode of motions to improve computational efficiency.The staggered Arakawa C grid was used for the variable arrangement.The vertical axis of POM was transformed to the σ-coordinate by where z is positive upward with the origin placed at the mean sea level, η is the sea level fluctuation, and H is the mean water depth.The horizontal viscosity and diffusivity were computed according to the formula described in Smagorinsky (1963).The vertical viscosity and diffusivity were determined by a level 2.5 turbulence-closure scheme (Mellor and Yamada 1982).
The computational domain was bounded between 115 -175.5°E and 10°S -23°N.The grid resolution in lateral space was approximately 2 arc-min (~3.7 km).The model topography was constructed using the revised ETOPO2 (US Department of Commerce et al. 2006)  both T and S values were set at 4000 m.The horizontal variation in the T and S fields was disregarded for excluding baroclinic flows.
Astronomical tides were incorporated into the model through prescribed tidal sea levels and tidal currents at the open ocean boundaries.The tidal sea levels and tidal currents at the open boundaries were composed by the OSU TOPEX/Poseidon Global Inverse Solution (Egbert and Erofeeva 2002).The constituents used to compose the tidal sea levels include eight primary constituents (M 2 , S 2 , N 2 , K 2 , K 1 , O 1 , P 1 , and Q 1 ), two long-period constituents (Mf and Mm), and three minor constituents (M 4 , MS 4 , and MN 4 ).The model run with tides was driven by tidal sea levels ten days before the 2011 Tohoku-Oki earthquake.
Ten experiments, Cases 1 through 10, using different model settings, were conducted to accomplish the objec- tive of this study.Table 3 lists the experiment settings for the simulation.The performance of the five seismic rupture models was tested in the first five experiments.The two significant aftershocks were included in Cases 1 to 6. Finally, the impact of stratification, tides, and the vertical resolution of the model itself were tested in Cases 7 to 10, respectively.

Tsunami Simulations based on the Varied Rupture Models
Figure 4 shows the maximal wave height of the simulated tsunamis from Models I to V. All five simulations reproduce a common steep peak near the epicentral region of the earthquake, and the maximal tsunami wave height forms  increases noticeably with the propagation distance (see Fig. 5), which may be caused by wave dispersion and diffusion effects (Burwell et al. 2007;Løvholt et al. 2010).The dispersive effect, which is nearly absent in the simulation of the 2011 Tohoku tsunami event by Løvholt et al. (2011), is similarly not included in the present model.The variation in tsunami waveform and associated wavelength from nearfield to far-field may be caused by wave diffusions (Burwell et al. 2007).
In comparison to the three models with well-determined rupture processes, the two simplified Models IV and V (Cases 4 and 5) show a limited scope for large tsunami heights (Fig. 4).Tsunami heights significantly decrease in approximately 1000 km.In fact, Models IV and V exhibit wider rupture surfaces than that of the other three models; the rupture amounts are consequently averaged out by a larger surface area and the lack of a precise spatial resolution.Among the five rupture models, the maximum slip in Model IV is similar to Model II.Notably, the tsunami propagation in Model IV is much more dispersive than in Model II.This is attributed to the rough spatial resolution in the rupture process.In comparing the tsunamis generated by the mainshock the tsunamis triggered by the two largest aftershocks, which are obtained by subtracting simulated tsunamis in Case 10 from that in Case 6, are relatively small (also see Fig. 4f).The maximal difference in simulated tsu- nami waves between Cases 1 to 5 is approximately 0.5 m.This is the value recorded in the vicinity of the mainshock's epicenter.
The bathymetric effect can be clearly observed from the simulations of the mainshock (Cases 1 to 5).As tsunamis propagate across the trench, the abrupt change in topography can cause a rapid change in phase speed and therefore cause a focusing/defocusing of the propagating waves.It was noted that along the Kuril, Japan, Izu Bonin, and Mariana Trenches, the tsunami height reaches approximately 0.6 m, which reflects the topographic effect in tsunami wave propagation.
Figure 5 illustrates the observed and simulated tsunami waveforms at the eight DART stations.Note that the sea-level variations contributed from tides are excluded from the observed and simulated data.It indicates that the simulated tsunamis triggered by the five rupture models in Cases 1 to Case 5 are mostly well-reproduced at seven of the eight DART stations, with the exception of Station 52406 (Fig. 1) where the arrival time of first tsunami peak is 7 min earlier than that of the first observed peak wave.The agreement in the observed and simulated tsunami phase and wave height at seven of the eight DART stations (500 to 3800 km from the epicentral region) confirms that the model effec-tively reproduces the wave propagation of the tsunamis.The remarkable model-data misfit in the arrival time at only the farthest Station, 52406, is presumably caused by the energy loss from tsunami to internal waves in the stratified ocean after such a long distance (~5300 km from the epicentral area) propagation, and by inaccurate tidal sea level conditions close to this station.Regardless of the overestimate in tsunami wave height at Stations 21416 and 21419, the tsunamis generated by Model III (Case 3) fit the observed tsunami waveforms better than that generated by Models I, II and IV, and V at the eight DART stations.The simulations in Cases 1, 2, 4, and 5 tend to underestimate the tsunami wave height at all eight DART stations.The initial pattern of sea-level displacement from Models I, II, IV, and V may be responsible for the higher level of model-data misfit than that from Model III.
Table 4 summarizes the first peak of the tsunamis from the observations and the simulations at the eight DART station for all cases in Table 3.The root-mean-square errors (RMSE) between the observations and simulations are calculated and averaged over the eight DART stations.According to the RMSE of Cases 1 to 5, Model III performs the best with a mean RMSE ~0.082 m in the tsunami simulation among the five selected rupture models.Cases 6 and 10 are designed to examine the contribution of the two aftershock-generated tsunamis to the observed tsunamis.Not surprisingly, the simulated tsunami wave height is small compared to that caused by the mainshock.The maximal tsunami wave height at the eight DART stations is smaller than 0.1 m, which accounts for only 5% of the observed maximal tsunami wave heights in the open ocean.Note that the RMSE is the same in Cases 6 and 10 because the timing of aftershocks lags behind the mainshock (i.e., aftershocks do not affect the first peak of the tsunami wave).The tsunamis contribute from the two aftershocks may be discerned only in the epicentral area, as shown in Fig. 4f.

Influences from Ocean Tides and Stratification with
Regard to the Simulation The long wavelength nature of tsunamis means that their propagation in the ocean mainly obeys the dynamics of shallow water waves.The phase speed of the waves, , where g is the gravitational acceleration and H is water depth, is approximately 170 m s -1 in water over 3000 m deep.This speed is much greater than the characteristic speed of many oceanic processes such as western boundary currents, tidal currents, etc.; consequently, tsunamis rarely interact with these processes, particularly in the vastness of the open ocean.However, in certain regions with abruptly changing topography, strong stratification, and shallow depths, for example, the nonlinear effects caused by the interactions between tsunami and oceanic processes may be important in modifying the tsunami waveforms (Kowalik et al. 2007;Santek and Winguth 2007;Kowalik and Proshutinsky 2010).
Based on Model III, we further examine the influences from the stratification (Case 7) and ocean tides (Case 8) to the simulated tsunamis.Case 9 is a model sensitivity run in which the vertical layers are set to 51 instead of 11 to check the influence of model's vertical resolution on the simulated tsunamis and also to verify the effects of the two oceanic process on tsunami propagation.The spatial variation in the simulated tsunami wave height, particularly in the near-field, in Cases 7, 8, or 9 is indistinguishable from the one shown in Fig. 4c and is thus not shown here.The simulated tsunamis are therefore insensitive to various vertical grid resolutions prescribed in the model for this study.The difference in the simulated wave height between Cases 7 and 8 and Case 3 occurs mainly at the far-field observations (Stations 52402, 52403, and 52405) where the simulated tsunami wave height is seemingly improved by < 0.01 m (see Table 4).The improvement, although minor, is in opposition to our general intuition as the oceanic parameters in Case 3 are more realistic than that in Cases 7 or 8, and thus the results of Case 3 are expected to be better than the latter two cases.It should be noted that all of the far-field stations localize to a complex bathymetry (Fig. 1), and the energy transfer from tsunamis to the generation of internal waves may decrease tsunami wave height after the abrupt change in topography.The wave height reduction due to tsunami-topography interactions under stratified ocean conditions may have led to an underestimation in far-field tsunami wave heights after a long-distance tsunami propagation in Case 3. The modeldata misfit in the tsunami simulation may have also resulted from the accuracy of the bathymetry and truncation errors in topography caused by the limited grid resolution.Summing up all of the results from Case 3 and Cases 5 to 9, we infer that the influence from varied factors of ocean conditions, which should include the bathymetry, can slightly affect the tsunami wave.The interactions compiled from all the factors may frequently be a result of setting off the amplification of tsunami waves caused by an individual influence.
The model-data misfit between Case 3 and Cases 7 and 8 is evaluated by a percentage difference.Case 3 was set as a standard run.The simulated first tsunami waveform in Cases 7 and 8 is compared to that in Case 3. The percent difference in simulated tsunami wave height is defined as: Comparing run Standard run Standard run -^h × 100%.Note that the difference between the simulated tsunami waveform with and without tides and stratifications is visibly indistinguishable, particularly for the first tsunami wave, at the eight DART stations.Figure 6 shows the distribution in the percent difference for Cases 7 and 8.Both results are in acceptable agreement for deep water regions.The largest differences in both experiments are found in shallow water regions, along the coasts, and in the nearly enclosed seas (Fig. 6).The influence of tides consists of tidal sea level variations and tidal currents, which may be more important than that of stratification.The maximal difference induced by the presence of tides reaches approximately 5% over shallow seas in the model domain and ~3% at DART Station 52406.In comparison, with the influence on the tsunami simulation by tides, the simulated tsunamis are generally less affected by stratification (< 1%) in the entire model domain (Fig. 6b).

SUMMARy
The five focal rupture models of the 2011 Tohoku earthquake determined from various seismic and geodetic observations were inversely examined using tsunami simulations.The influences of various oceanic factors were also evaluated in this study using a comparison between the simulated tsunamis and eight well-documented DART tsunami waveforms in the western Pacific.In the discussion of the initial conditions prior to the propagation of tsunami waves, those five distinctive seismic rupture models were adopted to generate the coseismic seafloor displacement.The consequent tsunami provocation was examined at the eight DART stations.In general, tsunami propagation is relatively insensitive to the temporal evolution of fault-plane rupture, but considerably influenced by the total amount and spatial pattern of slip in earthquakes, a carefully evaluated model can serve to efficiently differentiate the rupture model triggering the tsunamis.Consequently, Model III represented the tsunami waveforms better than the other four rupture models and implies that the main rupture should occur in a narrow range adjacent to the trench, rather than within the smooth asperities at relatively deep depths.
It was also observed that the far-field tsunami waveforms could be slightly disturbed by locally complex bathymetry (e.g., submarine ridges, continental shelf break, and nearly enclosed seas), stratified ocean and tides.In general, the latter two factors had rarely been considered in tsunami simulations.A series of numerical experiments quantified the model-data misfit caused by these two factors.The assumption of a homogeneous ocean in the simulation caused a < 1% discrepancy from the simulation under a stratified ocean.Alternately, simulations with the presence of ocean tides resulted in an approximately 3 -5 % difference in absolute wave height, compared to that without the influence of ocean tides.In comparison to the simulated tsunami waves converted from different rupture models, the influences of these two factors were small and can generally be disregarded.However, the combined influence of each factor in tsunami simulation is not linear and could have cas-cading effects.The combined effects in the tsunami simulation merit further study.

Fig. 1 .
Fig. 1.Plan view of the epicenters and DART stations.Red dots denote the numerous aftershocks in the 2011 Tohoku-Oki earthquake.The CMT focal mechanism solutions of the mainshock (05:45) and two large aftershocks (06:15 and 06:25) are presented in yellow.Four plates are joined around the epicentral region, which are the OK-Okhotsk Plate, the AM-Amur Plate, the PS-Philippine Sea Plate, and the PA-Pacific Plate.Eight DART (see red squares in the left panel for their locations) sea-surface height time series are arranged from top to bottom as their distance to the mainshock on the right panel.

Fig. 2 .
Fig. 2. Initial conditions of the tsunami sea-surface geopotential height.(a) to (e) are converted for Models I to V, and (f), the two aftershocks.The yellow star is the epicenter of the seismic source.The left and right color codes in meters are associated with the plots for Models I to V and for the aftershocks, respectively.
to propagate out of the computational domain.
Fig. 3. Temperature (blue) and salinity (green) profiles for the simulation with the ocean stratification.

Fig. 4 .
Fig. 4. Spatial distribution of the maximal tsunami wave height derived from the simulated tsunamis generated by the initial sea-level displacement for (a) Case 1, (b) Case 2, (c) Case 3, (d) Case 4, and (e) Case 5.The results from the tsunami simulation with the two largest aftershocks are shown in (f).The first peak level is indicated by the color code on the right of each panel.The topography is illustrated by gray shadings on each map for reference.The dashed lines represent the arrival time of the first tsunami peak.

Fig. 5 .
Fig. 5. Comparison of the observed and simulated tsunami time series at the eight DART stations.The black, yellow, red, blue, green, and purple lines denote the sea level variations caused by tsunamis from DART observations, and Cases 1 to 5, respectively.

Fig. 6 .
Fig. 6.Percentage difference for the simulated tsunami wave height between Case 3 and (a) Case 7 without the presence of stratification and (b) Case 9 without the presence of tides.The percentages are indicated by the color code on the right of each panel.The red squares indicate the locations of the eight DART stations.

Table 1 .
Parameters of the five rupture models selected for this study.

Table 2 .
Parameters used in the numerical model for the tsunami simulation.
simulation.Two large aftershocks, a M w 7.9 under-thrusting aftershock occurring at 06:15:41 and a M w 7.7 extensional faulting earthquake occurring at 06:25:51, were also considered in this study to assess the combined effects from multiple seismic sources in tsunami propagation (seeFigs.1,  2, and Table

Table 3 .
Settings of the numerical experiments.

Table 4 .
Characteristics of observed and simulated tsunamis.