Sensitivity of tsunami wave profiles and inundation simulations to earthquake slip and fault geometry for the 2011 Tohoku earthquake

In this study, we develop stochastic random-field slip models for the 2011 Tohoku earthquake and conduct a rigorous sensitivity analysis of tsunami hazards with respect to the uncertainty of earthquake slip and fault geometry. Synthetic earthquake slip distributions generated from the modified Mai-Beroza method captured key features of inversion-based source representations of the mega-thrust event, which were calibrated against rich geophysical observations of this event. Using original and synthesised earthquake source models (varied for strike, dip, and slip distributions), tsunami simulations were carried out and the resulting variability in tsunami hazard estimates was investigated. The results highlight significant sensitivity of the tsunami wave profiles and inundation heights to the coastal location and the slip characteristics, and indicate that earthquake slip characteristics are a major source of uncertainty in predicting tsunami risks due to future mega-thrust events.


Background
Massive tsunamis can destroy coastal cities and towns. Recent catastrophic events include the 2004 Indian Ocean tsunami and the tsunami generated by the 2011 off the Pacific coast of Tohoku earthquake (the 2011 Tohoku earthquake), which caused a large number of fatalities, damage to infrastructure, and huge economic loss (Murata et al. 2010;Fraser et al. 2013;Suppasri et al. 2013). To mitigate tsunami risk for coastal communities and built environments and to implement effective warning and evacuation protocols, accurate and reliable tsunami hazard assessments are essential. Techniques and methods for tsunami simulations have greatly advanced over recent years, and this has facilitated the development of operational tsunami early warning systems (Shuto 1991;Titov and Synolakis 1998;Synolakis et al. 2008;Tsushima et al. 2009). Tsunami inversion and simulation studies for the 2011 Tohoku event have demonstrated that observed tsunami wave and inundation characteristics can be reproduced by current tsunami simulation tools Lay et al. 2011;Yamazaki et al. 2011Yamazaki et al. , 2013Gusman et al. 2012;Mori et al. 2012;Grilli et al. 2013;Satake et al. 2013).
The major challenge in tsunami hazard and risk assessments is to reliably predict the occurrence and properties of future tsunamigenic earthquakes. In particular, earthquake slip has significant influence on the tsunami wave profile and the extent of inundation and run-up (Geist 2002;McCloskey et al. 2008;Løvholt et al. 2012;MacInnes et al. 2013;Satake et al. 2013). The uncertainty associated with earthquake source properties for future events is large, and thus, the assessment should include multiple scenarios that are likely to occur in light of available scientific evidence of past events. Such information on the possible range of predictions is valuable in developing robust design methods for coastal structures, tsunami-proof buildings, evacuation plans, and land use policies (FEMA 2008;Murata et al. 2010).
Earthquake slip is a complex phenomenon, governed by pre-rupture stress conditions, geometrical settings, and frictional properties of the fault that are largely unknown. Earthquake source inversions image the space-time evolution of the rupture by matching key features of synthetic data with observations. Typically, such slip distributions are obtained as a set of slip vectors over multiple subrupture sources, each of which is characterised by its fault geometry and temporal rupture process. The input data to source inversions vary, depending on specific studies, and are affected by available observational networks and their recordings (e.g. strong motion data, teleseismic data, geodetic measurements, and sea surface observations; Yokota et al. 2011). The quality and availability of data determine the practical limits of the developed slip model in terms of the usable frequency band and spatial resolution.
In the past, certain properties of slip models for various earthquakes, each derived differently using non-uniform data, alternative analysis methods, and variable parameterisations, have been investigated (Somerville et al. 1999;Mai and Beroza 2002;Lavallee et al. 2006). These studies show that a random-field model provides a powerful approach to parameterise the complexity of earthquake slip. In this case, the spatial slip distribution is conveniently characterised as a power spectral density in the wavenumber domain. For instance, Mai and Beroza (2002) investigated the validity of Gaussian, exponential, von Kármán, and fractal autocorrelation models by considering 44 slip distributions for worldwide crustal earthquakes. In their study, parameters of a random-field model (fractal dimension and correlation length) were related to macroscopic earthquake source parameters such as moment magnitude and fault length/width. The key idea for such analyses is to derive predictive models of complex earthquake slip distributions for future earthquakes that can be used for forward modelling of both tsunamis and strong motion.
The main objectives of the present study are (1) to develop stochastic earthquake slip models for the 2011 Tohoku mainshock based on the spectral analysis approach by Mai and Beroza (2002), and (2) to evaluate the impact of earthquake slip and fault geometry on tsunami simulation results in terms of near-shore sea surface profiles and inundation height. The investigation focuses on the 2011 Tohoku tsunami because rich datasets from global positioning system (GPS), ocean bottom pressure gauges, tidal measurements, and tsunami inundation/run-up survey results are available for validation of tsunami simulations Maeda et al. 2011;Mori et al. 2011Mori et al. , 2012Ozawa et al. 2011;Goto et al. 2012;Kawai et al. 2013). In addition, numerous source inversion studies have been conducted using tsunami, teleseismic, and geodetic data to constrain the rupture kinematics of the earthquake Fujii et al. 2011;Hayes 2011;Iinuma et al. 2011Iinuma et al. , 2012Lay et al. 2011;Shao et al. 2011;Yamazaki et al. 2011Yamazaki et al. , 2013Gusman et al. 2012;Sugino et al. 2012;Satake et al. 2013). This event therefore provides a unique opportunity to quantify epistemic uncertainty associated with predictions of the tsunami extent for different scenarios. We utilise slip models from different inversion studies and generate stochastic slip realisations for specific inverted models. We then compute initial wave profiles for the tsunami simulation (Okada 1985;Tanioka and Satake 1996) and finally obtain multiple realisations of maximum tsunami wave heights at various locations along the coast. Subsequently, tsunami hazard estimates such as near-shore wave profiles and inundation heights are assessed statistically, and their variability is evaluated for a given magnitude. The information on quantified variability is valuable for benchmarking the prediction uncertainty for future scenarios (for which no direct observations are available), thereby facilitating the rigorous validation of the tsunami simulation models and the quantification of uncertainty related to the predictions.
This study is organised as follows. The 'Methods' section provides a summary of the earthquake source models that were used in the tsunami simulation. The models include 11 published slip models for the 2011 Tohoku earthquake as well as realisations of stochastic slip models. The 'Results and discussion' section presents the tsunami simulation results using different earthquake slip models as well as slip distributions with different geometry parameters (strike and dip). The tsunami simulation was conducted by evaluating nonlinear shallow water equations considering run-up (Goto and Shuto 1983;Shuto 1991); the inundation calculation was performed by the moving boundary approach of Iwasaki and Mano (1979). The simulated results are compared with observed near-shore tsunami profiles at four GPS buoys and with the tsunami inundation/runup survey results along the Tohoku coast. Finally, the 'Conclusions' section provides the main conclusions of this study and describes future investigations that will be conducted.

Methods
This section presents an earthquake slip modelling procedure based on spectral analysis of an inversion-based earthquake slip distribution. The parameters derived from spectral analysis are then used in a spectral synthesis approach to generate random fields that represent earthquake slip with statistical properties equivalent to the inverted slip distribution. The analysis procedure is based on Mai and Beroza (2002). The 'Inversion-based source models' section provides a summary of the published inversion-based earthquake slip models for the 2011 Tohoku earthquake that we used for the subsequent spectral analysis and random-field generation. The 'Spectral analysis and synthesis of earthquake slip' section describes the details of the stochastic source modelling of earthquake slip distribution. Several modifications to the original method proposed by Mai and Beroza (2002) have been incorporated so that the method would be applicable to the Tohoku earthquake. These were necessary for three reasons: (i) the sub-fault size is larger (typically 10 to 50 km) than those for the 44 crustal earthquakes analysed by Mai and Beroza (2002); (ii) large slip occurs near the top edge of the fault plane along the Japan Trench; and (iii) we found that the marginal distributions of earthquake slip of inversion-based models for the Tohoku earthquake have a heavy right tail (i.e. relatively high probability at larger values, in comparison with the normal distribution; see Lavallee et al. 2006). The developed random-field generation method is designed to balance similarity in key features of the inversion-based models (e.g. overall slip distribution and its spectral characteristics) with dissimilarity of fine details (e.g. locations of large slip patches). Our approach is thus practical for studying the sensitivity of tsunami predictions with respect to earthquake source characteristics.

Inversion-based source models
Multiple earthquake slip models were collected from the literature, and these were applied to tsunami simulations to specify various initial tsunami profiles for the 2011 Tohoku event. A summary of the 11 models is given in Table 1; the corresponding slip distributions are shown in Figure 1. Some models (e.g. models 2, 6, and 11) have concentrated large slip patches (i.e. asperity) along the eastern edge of the fault plane, whereas others have large slip values near the epicentre. The maximum slip values for the 11 models vary significantly (from 35 to 75 m). The areas/shapes of the main slip asperity differ among the various models; for instance, in models 3, 5, and 7, slip concentrations extend primarily along strike, whereas models 1, 4, and 8 have circular slip concentrations. It is noteworthy that these slip models were developed using different types of data (teleseismic/tsunami/geodetic) to constrain the slip distribution over the fault plane. Therefore, details of the slip models and fault geometry may differ considerably (e.g. down-dip and southern extent of the fault plane).
These slip models reflect not only aspects of the complex rupture process of the 2011 Tohoku mainshock, but also the uncertainty in estimating the resulting slip distributions. For specific observations (e.g. tsunami waveform at a particular tidal gauge), some models may perform better than others. In predicting future tsunami scenarios, post-event model evaluations are not possible and a priori performance tests of individual models are not always feasible to perform. Therefore, it is reasonable to treat all models equally once they are calibrated against observed geophysical data. However, for some cases (e.g. where there is higher priority on tsunami inversion studies over seismic inversion studies), unequal weighting of the candidate models may be considered. The resulting variability in the tsunami predictions can be interpreted as epistemic uncertainty due to modelling errors when using different data and methods. Recognising that the 2011 Tohoku earthquake is extremely well recorded in comparison with other major earthquakes, the epistemic uncertainty for the 2011 Tohoku event may be considered relatively small given the current observation networks and inference techniques (i.e. empirical benchmark of epistemic uncertainty). Nevertheless, it is important to investigate the uncertainty of tsunami predictions using existing slip models. In addition, we generate stochastic variation of the inverted earthquake slip models to assess the sensitivity of tsunami simulations to the slip distribution (note that this is also applicable to strong motion simulations).

Spectral analysis and synthesis of earthquake slip
The Mai-Beroza method (Mai and Beroza 2002), which was originally developed for M w 6-8 crustal earthquakes, was modified to address issues specific to the M w 9 Tohoku event. The method consists of two parts: (1) spectral analysis of a slip model to characterise the spatial complexity of earthquake slip and (2) spectral synthesis to generate multiple realisations of the random-field-type slip distribution (Pardo-Iguzquiza and Chica-Olmo 1993). In the present spectral analysis, we used the anisotropic two-dimensional (2D) von Kármán auto-correlation function to capture the spatial characteristics of a given slip distribution; the power  Figure 1 Eleven inversion-based earthquake slip models with colour-coded final displacement (slip in metres) over the fault. The gridding displays the fault plane discretisation used in the slip inversion. Note the variability in the location and size of the rupture plane (the Pacific coastline of Japan is shown as a reference), and the differences in the location, size, and displacement of the large slip patches. spectrum P(k) of the von Kármán auto-correlation function is given by: where k is the wavenumber; k = (A z 2 k z 2 + A x 2 k x 2 ) 0.5 ; A z and A x are the correlation lengths for the down-dip and along-strike directions, respectively; and H is the Hurst number. This representation describes slip heterogeneity as a function of the wavenumber; A x and A z control the absolute level of the power spectrum in the low wavenumber range (i.e. when k ≪ 1), while H determines the slope of the power spectral decay in the high wavenumber range (Mai and Beroza 2002).
Our modifications to the Mai-Beroza approach include the use of a larger grid spacing (i.e. 10 to 25 km) in comparison with the original method (i.e. 1 km) and nonlinear transformation of the synthesised slip distributions to generate heavy right-tail characteristics found in the original slip distributions (i.e. source inversion models as published in the literature; Figure 1). Below, we provide further explanations for the modified spectral analysis and synthesis method. In particular, we have tested the effects of tapering (i.e. adding rows and columns of slip at the edges of the fault plane to achieve smooth transition at the fault boundary) and interpolation of the slip distribution. By investigating four auto-correlation models (exponential, Gaussian, von Kármán, and fractal), we have found that the anisotropic 2D von Kármán model provides the required flexibility to characterise the properties of the observed power spectra of the inverted slip distributions. A step-by-step procedure for spectral analysis and synthesis is described in the Appendix.
Once the parameters of the von Kármán model are obtained from the spectral analysis of a slip model, we deploy a stochastic random-field generator based on a Fourier integral method (Pardo-Iguzquiza and Chica-Olmo 1993). This approach generates a second-order stationary random field whose covariance structure is defined by the von Kármán model. The synthesised random field has quasi-normal marginal distribution with a desired spatial correlation structure, while the original slip distribution for the Tohoku earthquake has heavier right-tail characteristics than the normal distribution. To equalise the characteristics of the original and synthesised slip distributions, we implemented a scaling and modification procedure for the synthesised slip values by (i) using the Box-Cox power transformation and (ii) adjusting the mean and standard deviation of the transformed slip distribution. We thus ensured that the statistics of the synthesised slip model were similar to those of the original slip model. The Box-Cox analysis identifies the best power parameter such that the non-normal random variable becomes similar to the normal random variable after transformation. The Box-Cox transformation of a variable x is defined as: where λ is the power parameter (note that when λ = 0, y = log(x), i.e. the logarithmic transformation). An optimal value of λ can be estimated by calculating the linear correlation coefficient between the transformed variable y (for a given λ) and the standard normal variate. We applied nonlinear scaling to the synthesised slip distributions (generated in the wavenumber domain with the estimated von Kármán parameters and then transformed into the spatial domain via a 2D inverse fast Fourier transform (FFT)) using the Box-Cox transformation. As this transformation may lead to locally very large slip values, we constrained its upper bound by the maximum of the original model. The corresponding mean and standard deviation were then adjusted to ensure that the seismic moment for the simulated slip field was similar to that for the original slip model. Moreover, to ensure that the synthesised slip distributions resemble the original slip models, we added a constraint on the location and extent of the rectangular high-slip patch (asperity zone). Through a pattern-matching algorithm, we selected the slip realisations from a large suite of randomised slip models that satisfied the given criteria on the spatial slip concentration (see the Appendix for further details).
The spectral analysis and synthesis of a random field are illustrated for the slip distribution obtained by Satake et al. (i.e. model 2) in Figures 2, 3, 4, 5. Their slip model has a grid spacing of 25/50 and 50 km and fault dimensions of 200 and 550 km in the down-dip and alongstrike directions, respectively ( Figure 1). By adopting a spatial spacing of 25 km for interpolation and tapering, and then carrying out the 2D FFT, we obtained the normalised power spectrum (Figure 2a). Its usable wavenumber ranges extend from 0.01 to 0.126 km −1 and from 0.01 to 0.063 km −1 for the down-dip and alongstrike directions, respectively. Figure 2b displays the parameter estimation approach for the von Kármán model, returning an overall estimate of H = 0.82 based on the spectral decay fitting (leading to H = 0.61) and a 2D search of H and A for the circular average (leading to H = 1.1). Subsequently, we found values for A z and A x that were 56 and 107 km in the down-dip and alongstrike directions, respectively ( Figure 2c). Hence, the correlation length in the along-strike direction is greater than that in the down-dip direction (as was anticipated from visual inspections of the slip model). Generally, the misfit function between the observed and modelled power spectra shows a clear minimum (Figure 2c), which indicates that the estimated correlation lengths are stably estimated for a given Hurst number. On the other hand, the value of the Hurst number is less well constrained by the data. This is why we applied two estimation methods and then adopted the geometric mean of the two estimates of H as the representative value (see step 4 in the Appendix). The marginal distribution of the original slip model is shown in Figure 3a, highlighting the heavy right tail of this slip model. The Box-Cox analysis for the original model ( Figure 3b) yields an optimal power parameter of λ = 0.2. Figure 4 shows one random-field realisation that has the desirable characteristics of the slip distribution (after nonlinear scaling and adjustments of the mean and standard deviation). Notice that the final synthesised slip model has an asperity cell (star) within the specified asperity zone of the original slip model. Here, we define the asperity zone as sub-faults having slips greater than three times the average slip (black rectangles). More than 30% of the total slip is concentrated around the asperity that occupies only about 10% of the total rupture area, which is consistent with the results of Mai et al. (2005). Note that the dimensions of the asperity rectangle and the slip concentration are determined based on the slip distribution of the original model. Using the above-mentioned procedure, we generated 50 realisations of the target slip distribution; their average probability mass function (Figure 3a) indicates that the marginal distributions of earthquake slip for the original and synthesised models are similar. To demonstrate that the synthesised slip distributions exhibit similar spectral characteristics to the original model, Figure 5 compares normalised power spectra of the original, fitted, and simulated spectra. These results indicate that our modified spectral synthesis procedure achieves satisfactory performance.
Finally, we apply our modified spectral analysis and synthesis procedure to the 11 inversion-based slip models (Table 2). Table 2 also includes the Box-Cox parameter and other relevant statistics and information for spectral synthesis. The majority of the Hurst numbers are indicated as 0.99 in order to constrain H so that it falls within a physically meaningful range, although a better fit can be achieved by allowing H > 1.0. The Box-Cox parameter is relatively stable across different slip models; the logarithmic transformation case, having a heavier right tail than λ = 0.2 to 0.3, was obtained from model 11 (which has large slip values that exceed 70 m for some cells).

Tsunami simulation setting
We utilised a tsunami simulation tool (Goto and Shuto 1983;Shuto 1991) that is capable of generating off-shore tsunami propagation and inundation profiles by evaluating nonlinear shallow water equations considering runup using a leap-frog staggered finite difference scheme. For simulating the 2011 Tohoku tsunami, information on bathymetry, surface roughness, and coastal defence structures was obtained from the Cabinet Office of the Japanese Government. To compute initial water surface elevation for the tsunami simulation, analytical equations for elastic dislocation by Okada (1985) together with the equation by Tanioka and Satake (1996) were used. The latter takes into account the effects of horizontal movements of the steep seafloor on the vertical water dislocation. For all cases, fault rupture is assumed to occur instantaneously (note that information on the detailed space-time rupture evolution is available only for a subset of the inversion-based models listed in Table 1). We carried out tsunami simulations considering three nested domains (one 1,350-m domain and two 450-m domains) and wave-propagation duration of 2 h (integration time step Δt = 2 s). Note that we also examined tsunami simulations for several cases using finer grid sizes (i.e. 150 and 50 m); the results in terms of near-shore wave profiles at the GPS buoys and inundation heights along the coast were similar. Computational cells include those on land, and coastal defence structures are taken into account using an overflowing formula as a sub-grid model. Based on our simulations, we then assessed the sensitivity of wave profiles at four GPS buoys (Kawai et al. 2013) and inundation heights along the Tohoku coast (Mori et al. 2011  Colour-coded in blue and red are the spectra of the inverted model and the corresponding fitted spectra, respectively. Thin grey lines depict the spectra of 50 simulated slip realisations with the corresponding median and 16%/84% percentiles, which are shown as black solid lines and black dotted lines, respectively. distributions, we generated 61 slip models. They are (i) the original slip model, (ii) 5 models with different strike angles with respect to the strike of the original model, (iii) 5 models with different dip angles with respect to the dip of the original model, and (iv) 50 models that were generated using the preceding spectral analysis and synthesis procedure (see the 'Inversion-based source models' section). For (ii), strike angles were varied by −5°, −2.5°, 2.5°, 5°, and 7.5°w ith respect to the original strike, while for (iii), dip angle deviations of −2.5°, 2.5°, 5°, 7.5°, and 10°with respect to the original dip value were considered. These ranges were determined by inspecting various inversion-based models and the regional seismotectonic setting. We only changed one model component at a time. The full sensitivity analysis, which would vary all major parameters simultaneously and cover the parameter space comprehensively, will be the topic of a future study.

Tsunami simulation sensitivity: Satake et al. model
This section presents the tsunami sensitivity results for the Satake et al. model (i.e. model 2). The results based on different slip models will be discussed in the 'Tsunami simulation sensitivity: all inversion-based models' section. The main reasons for selecting this model as a base case are that the source inversion utilises tsunami data and that information on kinematic earthquake slip is available (thus, the influence of adopting instantaneous slip rather than kinematic slip can be assessed). Moreover, the tsunami simulation codes used in this study and by Satake et al. (2013) are similar in terms of the governing equations and numerical discretisation. This setup essentially closes the gap between forward analysis (i.e. tsunami prediction) and inverse analysis (i.e. source modelling), and results in consistent reproduction of the results presented by the original study. However, this is not always the case when simulation codes with different formulations are employed in forward and inverse analyses. This applies to the Yamazaki et al. model (model 6), where a nonhydrostatic model is used to account for wave breaking (Yamazaki et al. 2009). Figure 6 shows the maximum tsunami wave heights for instantaneous and kinematic earthquake slip cases. It also compares the observed and simulated tsunami wave profiles at four GPS buoys, i.e. Miyagi Central, Miyagi North, Iwate South, and Iwate Central (from South to North). The locations of the four GPS buoys are indicated in the maps. It is important to realise that maximum wave heights for instantaneous and kinematic ruptures are significantly different (Figure 6a). Large differences can be observed along the Japan Trench, where a large earthquake slip has occurred, and off the Miyagi Prefecture (near Miyagi Central buoy). The former is due to the concentration and timing of the large slip occurrence (i.e. earthquake slip for the kinematic case is more spread out), whereas the latter represents the combined effects due to the tsunami source and propagation.
The assumption of instantaneous slip leads to the generation of steep short-wavelength components of tsunami waves. This can be clearly seen in the tsunami wave profiles at the four GPS buoy locations (Figure 6b), and it is particularly evident in the large local peak of the tsunami wave at Miyagi Central. Overall, the tsunami predictions based on the original Satake et al. model agree well with the tsunami observations. The wave profiles for the kinematic slip case have relatively smooth temporal changes in the profiles, and they match well with the observed ones (as demonstrated by Satake et al. 2013). The results for the case of instantaneous slip have long-wavelength characteristics in regard to the tsunami wave profiles, and while these are broadly similar to the observations, they over-predict local peaks of the first major waves (around 2,000 s). The notable difference between the instantaneous and kinematic cases, which was observed for the Satake et al. model, cannot be generalised for other inversion-based slip distributions. For instance, kinematic models by Shao et al. (i.e. models 3 to 5) produce similar results as the instantaneous models, and the results contain large short-wavelength components, because a large amount of slip is released within a rather short time interval, thus exhibiting close-toinstantaneous rupture behaviour.
Next, we compare the tsunami simulation results at the four GPS buoys when considering different fault geometry parameters (i.e. strike and dip). The results for the strike and dip variations are shown in Figures 7 and 8, respectively. In the figures, in addition to the overall variability of the tsunami wave profiles, temporal close-ups show the changes of the profiles at the arrival of the first major waves. The results shown in Figure 7 indicate that strike variations mostly have little impact on the peak amplitude of the tsunami profile (note that this result should be cautiously interpreted because some focusing may occur owing to the strike variation); however, arrival times of the first major waves at Miyagi Central can be significantly altered by strike variations (where the effects due to large earthquake slip along the Trench are particularly remarkable). For the considered range of strike angles, increasing the strike with respect to the base case advances the arrival time of the first major waves (e.g. by about 250 s for an increase of 7.5°). This can be important for tsunami evacuation planning. Figure 8 shows the effects of steeper dip angles on tsunami wave heights (depending on the locations). Variation of the dip angle from −2.5°to 10°with respect to the base case results in differences in peak amplitudes of the tsunami wave profiles by a factor of 2. These changes were intuitively expected, as steeper dip angles increase the vertical seafloor deformation for a given earthquake slip.
However, a quantitative estimate of this variation, in particular when the fault is generally dipping gently (note that the precise dip angle over the shallow parts of the fault is not known or the dip changes along strike), is mandatory to assess the uncertainty of, and hence potential maximum of, wave heights that can be expected.
Next, we investigated the variability of the tsunami wave profiles due to uncertain earthquake slip using the simulation results based on 50 stochastic slip distributions (see the 'Spectral analysis and synthesis of earthquake slip' section). The results (Figure 9) demonstrate that tsunami wave profiles are sensitive to details of the slip distributions, and their variability due to stochastic earthquake slip is much greater than that due to variations in strike and dip angles (Figures 7 and 8). The results also suggest that wave amplitudes contained in the 16% and 84% curves (not drawn explicitly in Figure 9) vary at least by a factor of 2 (when peak values are attained at around 2,000 s), and these tend to bracket the simulated tsunami wave profiles based on the original slip model. This clearly indicates that the spatial distribution of earthquake slip has a major effect on the tsunami predictions, and it represents an important source of uncertainty for future scenarios. Figure 9 also displays four specific cases of earthquake slip distributions and their corresponding tsunami wave profiles. The four cases were selected such that the synthesised slip distributions have asperity areas in different parts of the fault plane; realisation 8 has an asperity close to the epicentre (i.e. more slip towards the downdip part of the fault plane in comparison with the original model), realisation 17 has a slip concentration in the northern part of the fault plane, realisation 29 has an asperity towards the southern part of the fault plane, and realisation 34 has a large slip concentration along the Trench, similar to the original slip model. By comparing the wave profile results for the four slip cases, we can therefore correlate the effects of asperity locations with the simulated tsunami wave profiles. The results for the specific cases show that realisation 17 generates large tsunami wave heights at Iwate Central and Iwate South (as expected), while realisation 29 results in less significant tsunami wave profiles at the four GPS buoys than the observations (however, much greater tsunami waves are generated in Fukushima and further South). Comparison of the results for realisations 8 and 34 indicates that when large slip is concentrated along the Trench (at shallow depths), local short-wavelength components tend to be dominant and produce large peak values. These observations are useful for understanding the influence of the spatial slip distribution on the tsunami predictions for different inversion-based models. Figure 10 compares the simulated tsunami inundation heights along the coast with strike/dip/slip variation to the Tohoku Tsunami Joint Survey (TTJS) data (Mori et al. 2011(Mori et al. , 2012. The results focus on the range of the tsunami predictions due to different source parameters (rather than individual results). To obtain the inundation height profiles, discrete cells containing zero elevation were set up along the coastline and peak values were extracted within individual cells. Figure 10a shows the tsunami inundation for variations in strike (−5.0°to 7.5°with respect to the original model). The degree of variability of the inundation heights is typically 5 to 10 m (depending on the location). The degree of variability of the inundation heights becomes greater (typically 5 to 15 m; Figure 10b) when the dip is varied (−2.5°to 10.0°with respect to the original model). The increased variability for dip was observed along the rias coast of the Miyagi and Iwate Prefectures (from 4.25 × 10 6 to 4.45 × 10 6 in terms of northing). This particular part of the coast was formed by submerged steep and narrow river valleys, and therefore, tsunami waves tend to be amplified by local bathymetry. A source model with steeper dip angles can match the observed inundation height better. The variability of the inundation heights due to stochastic slip distributions is much greater than that for the cases of variable strike or dip (Figure 10c). The range of the predicted inundation heights contains the TTJS observations comprehensively. The variability of the predicted inundation heights strongly depends on the location (i.e. northing). The southern and central parts of the Tohoku region (from 4.05 × 10 6 to 4.25 × 10 6 and from 4.25 × 10 6 to 4.38 × 10 6 for the southern and central parts, respectively) show relatively consistent trends. In contrast, the northern part of the Tohoku region (i.e. from 4.38 × 10 6 to 4.45 × 10 6 ) is subjected to much greater variability.
To quantify the variability of the tsunami simulation results due to slip variations, the mean and coefficient of variation (CoV) of the inundation heights for given locations are shown in Figure 11. The mean trend agrees with the TTJS results in the southern and central parts of the Tohoku coast, while a large discrepancy was observed in the northern part. The variability of the predictions (i.e. CoV) is spatially variable; in the northern and southern parts (i.e. Iwate and Fukushima Prefectures), the CoV is about 0.5 to 0.6, whereas in the central part (i.e. Miyagi Prefecture), it ranges from 0.2 to 0.3.
It is important to note that the spatial distribution of the CoV shows a coherent large-scale trend that is independent of the tsunami height and the source models. This result has important implications in tsunami hazard and risk assessments for coastal structures from a disaster risk reduction perspective. In current design manuals for coastal defence and evacuation structures (e.g. FEMA 2008), the tsunami safety margin is not specified as a spatially dependent quantity, and thus, the designed structures are likely to have non-uniform tsunami hazard/risk. Therefore, consideration of possible ranges of tsunami inundation heights can facilitate the implementation of more uniform and consistent tsunami risk reduction measures.

Tsunami simulation sensitivity: all inversion-based models
The analysis carried out for the Satake et al. model (see the 'Tsunami simulation sensitivity: Satake et al. model' section) has been extended for all 11 inversion-based slip models (Figure 1 and Table 1), and the results are discussed in this section. Figure 12 shows the tsunami wave profiles at the four GPS buoys for the 11 slip distributions. To enable the inspection of individual wave profiles, the results are displayed in two panels with different colours. The inversion-based slip distributions (Figure 1) can be categorised into several groups: (i) models 2, 4, 5, 6, and 11 have large slip amplitudes (more than 50 m); (ii) models 2, 6, and 11 have large slip concentrations along the Trench, whereas models 1, 3 to 5, and 7 to 10 have large asperity areas in the down-dip direction away from the upper edge of the fault plane; and (iii) models 3, 5, 7, and 10 have   Figure 12 Variability of tsunami wave profiles at four GPS buoys for different models. Models 1 to 6 (a) and models 7 to 11 (b). Note the overall agreement of all model predictions with the observations (black line) at long wavelengths, but the substantial variations in peak amplitude and timing (in particular at the Miyagi sites).
asperity areas that are more stretched along the strike direction (which can be observed in the estimated values of the correlation lengths shown in Table 2). A detailed inspection of Figure 12 suggests that slip models 2, 6, and 11 with large deformation near the upper edge of the fault plane generate distinct spikes in the tsunami profiles at the Miyagi Central and Miyagi North buoys, which are different from the observations at those sites. Note that the differences can be attributed to several reasons including instantaneous versus kinematic rupture cases and the use of hydrostatic versus hydrodynamic tsunami calculations when constructing the slip distribution from tsunami observations. The former is applicable to model 5 (as demonstrated in Figure 6), while the latter is applicable to model 6 whereby the consideration of hydrodynamic responses of the ocean filters out short-wavelength components (Løvholt et al. 2012). The slip models 3 to 5, which have relatively large slip over large areas, over-predict the observed tsunami wave profiles at the Miyagi buoys, whereas the slip models 1 and 7 to 10, which have moderate amounts of slip, produce tsunami wave profiles at the Miyagi buoys that are broadly similar to the observations. At the Iwate buoys, good agreement between the simulated and observed tsunami waves was achieved for models 1, 2, and 6 (including phase-arrival times) and for models 3 to 5 and 10 to 11 (mainly peak amplitudes). The predictions at the Iwate buoys based on models 7 to 9 lack local peaks because of the absence of large slip patches in the northern part of the fault plane, while the general long-wavelength fluctuations are captured well. The resemblance of the simulated tsunami wave profiles with respect to the observations can be understood from the slip distributions ( Figure 1). It is interesting to note that as an ensemble, the average predictions of the 11 slip models match well with the observations both in terms of amplitudes and timing (data not shown for brevity). Nevertheless, the variability of the predicted tsunami wave profiles can be as large as a factor of 3 (when the peak values are attained). Such results are useful for benchmarking epistemic (model) uncertainty associated with tsunami predictions. Next, we examine the sensitivity of wave profile predictions at the four GPS buoys while considering strike, dip, and slip variations together with different original slip distributions. The overall trends of the predicted profiles are similar to the results in Figures 7, 8, 9. Because the ensemble mean of predicted wave profiles using the 11 inferred slip distributions broadly captures the amplitude and timing of the observed tsunami waves, the predictions generally reproduce the observations. Remaining differences are due to additional variability contributed by the variation of the original slip distributions. To illustrate the variability of the tsunami predictions due to slip distribution, the results for slip variations are presented in Figure 13, whereby the data are displayed as individual curves (grey lines) for each rupture model, and the median and 16% and 84% percentiles based on the entire results in comparison with the observed wave profiles are also shown. Figure 13 indicates that the median agrees well with the observed profiles, particularly at the Miyagi buoys, whereas the 16% and 84% curves bracket the observed wave profiles at the Iwate buoys (note that the local peaks of the median predictions are lower than the observed profiles). Figure 14 shows the variability of tsunami inundation heights along the coast for the 11 slip distributions (separated into three panels for enhanced visibility). Their ability to match the observed inundation heights depends on the location of the coastal line of interest. In the southern part of the Tohoku coast, models 1 to 2, 6, and 10 generally produce predictions consistent with the TTJS results, whereas major over-/under-predictions occur for models 7 to 9. In the central part, models 1 to 4 and 10 to 11 agree well with the survey results, while larger inundation predictions are made for models 5 to 6 and 11, which tend to be associated with the models with large slip values off Miyagi Prefecture (Figure 1). In the northern part, differences of the predicted inundation heights are greater than in other parts of the Tohoku region. Models 2, 5 to 6, and 10 to 11 produce predictions that are slightly lower than the observed heights but broadly consistent with the survey results, whereas some models such as 7 to 9 generate significantly lower inundation heights than those seen in the observations and other predictions. This leads to wider variability of the tsunami predictions along the northern coast, when multiple original slip distributions are considered in addition to the strike/dip/slip variation ( Figure 15). Essentially, the conclusions that can be drawn from Figure 10 are applicable to the cases where the multiple inversion-based models are taken into account. The summary statistics of the inundation heights for different locations along the coast have already been presented in Figure 11. The results are similar to the cases where model 5 is considered for generating stochastic random-field slip distributions. The mean trend for the 11-model case is slightly decreased from the single-model case; at some locations, large increases of the mean inundation profile are observed because of the inclusion of slip models with large slip off Miyagi Prefecture. The CoV values for the 11-model case are similar to those for the single-model case. Therefore, the discussion regarding the implications for coastal structural design is generally applicable to the 11-model case.

Conclusions
This study investigated the sensitivity of tsunami hazards to different earthquake slip distributions and fault geometry for the 2011 Tohoku earthquake. The main objective was to quantify through numerical simulations the possible range of tsunami properties, focusing on nearshore tsunami wave profiles and on-shore tsunami inundation heights. We explored the sensitivity of the tsunami hazards to uncertain fault geometry by varying strike and dip, and to uncertain earthquake slip by considering stochastic slip distributions and inversion-based slip models. The method for stochastic slip was based on the spectral representation of slip heterogeneity and the spectral synthesis procedure for generating random fields, which is a method that was originally developed by Mai and Beroza (2002). The stochastic model was then employed to generate synthetic earthquake slip distributions. To ensure that the synthesised slip distributions were realistic and comparable to the original slip distribution, several pragmatic constraints such as edge tapering and nonlinear transformation were implemented. Subsequently, we used the synthesised earthquake slip distributions for tsunami simulations. The results were compared with actual observations of the 2011 Tohoku event such as GPS buoy recorded wave profiles and TTJS survey data on the land. The available extensive geophysical and tsunami datasets for the 2011 Tohoku event, compared to other major tsunami events, facilitated our detailed validation study. It also enabled us to quantify the epistemic uncertainty associated with tsunami predictions based on various source inversion models. Hence, the outcome of this study may be useful for benchmarking estimation uncertainty for future events empirically.
Comparisons of the tsunami simulation results with the observed data indicate that wave profiles at four GPS stations are well reproduced by the ensemble of the 11 inversion-based slip distributions, whereas the variability of the predictions is large at locations that face large slip/asperity areas (e.g. central/northern Miyagi Prefecture and southern Iwate Prefecture). Our results highlight strong sensitivity of tsunami wave heights to site location and slip characteristics, and also to variations in dip. While dip variations can be accounted for semideterministically (generally, the dip of a subduction zone is relatively well known and may only change gradually with depth or along strike), earthquake slip characteristics are the major source of uncertainty for predicting tsunami risks due to future mega-thrust events. Consequently, the spatial variability of tsunami hazards has important implications for the tsunami-resistant design of structures and evacuation plans because different safety margins may be required to achieve consistent reliability levels. Another relevant aspect that needs attention in tsunami sensitivity analyses is the consideration of new scaling relationships for fault dimensions of mega-thrust subduction earthquakes (e.g. Murotani et al. 2013). Therefore, future investigations should assess the influence of earthquake slip distributions and fault geometry when quantifying the variability of impact assessments for tsunami hazards.

Appendix
This section describes the spectral analysis and synthesis procedure in a workflow manner to further illustrate the method we used to generate stochastic slip realisations that resemble inverted slip distributions.
1. Read a cell-based slip model as a 2D matrix, convert it to a grid-based slip distribution, and bilinearly interpolate the slip distribution using a selected grid spacing. The new grid spacing should be related to the grid size of the input slip distribution. Preliminary analyses show that half of the original grid spacing may be a reasonable choice. 2. Taper the slip distribution by adding three rows and columns over which slip decreases to zero on each side of the fault plane. Note that tapering changes the power spectral characteristics. As part of preliminary analyses, different tapering cases (including no tapering and the use of mirror image at fault boundaries with large slip) were considered. When large slip values are concentrated at the along-strike edge, we find that the impact of no tapering is significant. 3. Calculate the 2D FFT of the slip distribution and obtain the 2D normalised power spectrum (computed as the squared amplitude spectrum). Define an applicable wavenumber range for spectral analysis by considering the characteristic dimension of the fault plane for the lower limit and the spatial resolution of the original slip model for the upper limit. 4. Extract the normalised power spectra in the downdip and along-strike directions, and fit the von Kármán spectral model (Equation 1) to estimate A z , A x , and H. To develop an anisotropic 2D von Kármán model, we first estimate the Hurst number using two methods: linear regression analysis of spectral decay for the circular average Anguilar et al. 1993) and a 2D search over the possible parameter space of A z , A x , and H by minimising the L 2 -norm (Mai and Beroza 2002). Based on the results for the two methods, the geometric mean of the estimated H values is adopted as a final value of H, and then A z and A x are obtained for each direction based on the onedimensional (1D) search. This pragmatic approach makes use of the fact that the circular averaged power spectrum is more stable for parameter estimation than the down-dip and along-strike spectra. Estimated H parameters usually fall between 0 and 1 in agreement with physical theories (Geist 2002;Mai and Beroza 2002). In the above approach, anisotropic features of slip distribution in the downdip and along-strike directions are captured solely by A z and A x . 5. To generate a realisation of a slip distribution with desired stochastic properties, compute the theoretical normalised power spectrum for given parameter values, generate the wavenumber matrix, calculate the power spectrum matrix, and obtain the amplitude spectrum. For the same wavenumber space, generate the random-phase matrix between 0 and 2π. Subsequently, the constructed matrix of complex Fourier coefficients is transformed into the spatial domain via a 2D inverse FFT. To ensure that the spatial random field generated by this approach is real, Hermitian symmetry is enforced when constructing the randomphase matrix. 6. The generated synthetic slip distribution is modified to have a marginal distribution of slip values with statistics (mean and standard deviation) similar to the original slip model. For this purpose, firstly, the location and extent of an asperity (rectangle) for the synthesised slip field are constrained to resemble those of the original model, whose asperity zone is defined as a set of sub-faults that have slip values greater than a specified threshold value. In this study, the threshold value was set to three times the average slip. The asperity rectangle of the synthesised distribution is specified as fractions of the fault length and width, whereas the extent of the slip concentration around the asperity is specified as a percentage of slip within the asperity rectangle with respect to the total sum of slip over the fault plane. The asperity parameters (i.e. dimensions and slip concentration) are determined based on the original slip distribution. Typically, fractions of fault length and width are 0.35 to 0.45 and 0.25 to 0.35, respectively, and the slip concentration percentage is 30% to 40%. These are broadly consistent with the results by Mai et al. (2005). To satisfy these requirements, multiple random fields are generated from the same set of parameters until the asperity of the synthesised field falls within the designated asperity zone and its spatial concentration satisfies the criteria. 7. Nonlinear scaling of the generated slip values is performed using the Box-Cox transformation. As this transformation potentially leads to very large slip, a realistic upper bound (i.e. maximum slip of the original slip model) is introduced by limiting the synthesised slip values. Then, the mean and standard deviation of the slip values are adjusted to ensure that the amount of seismic moment for the simulated slip field is similar to that for the original slip model. To allow for some variation of the slip statistics, random multiplication factors for the mean and standard deviation are generated from the uniform distribution with lower and upper limits (the range between the lower and upper limits is set to 0.1 in terms of either the target mean or standard deviation), and these are used to determine the slip statistics.