Shallow slip amplification and enhanced tsunami hazard unravelled by dynamic simulations of mega-thrust earthquakes

The 2011 Tohoku earthquake produced an unexpected large amount of shallow slip greatly contributing to the ensuing tsunami. How frequent are such events? How can they be efficiently modelled for tsunami hazard? Stochastic slip models, which can be computed rapidly, are used to explore the natural slip variability; however, they generally do not deal specifically with shallow slip features. We study the systematic depth-dependence of slip along a thrust fault with a number of 2D dynamic simulations using stochastic shear stress distributions and a geometry based on the cross section of the Tohoku fault. We obtain a probability density for the slip distribution, which varies both with depth, earthquake size and whether the rupture breaks the surface. We propose a method to modify stochastic slip distributions according to this dynamically-derived probability distribution. This method may be efficiently applied to produce large numbers of heterogeneous slip distributions for probabilistic tsunami hazard analysis. Using numerous M9 earthquake scenarios, we demonstrate that incorporating the dynamically-derived probability distribution does enhance the conditional probability of exceedance of maximum estimated tsunami wave heights along the Japanese coast. This technique for integrating dynamic features in stochastic models can be extended to any subduction zone and faulting style.

Probabilistic Tsunami Hazard Analysis (PTHA) requires the computation of numerous tsunami scenarios for each specific area [1][2][3][4][5] . In particular, for earthquake-generated tsunamis, that is for Seismic PTHA (SPTHA 2 ), spatial slip variation greatly influences tsunami impact on the coastlines located in the near field of the fault [6][7][8] . Producing probabilistic inundation maps starting from dynamic earthquake simulations is a computationally expensive problem. As a workaround, standard practice is to produce suites of likely earthquake slip using stochastic slip distributions which are based on general features observed across a wide range of geological and tectonic settings [9][10][11] . Generating stochastic models based solely on finite fault inversions from subduction zone earthquakes has shown to better account for inundation compared to standard models; however they still tend to underestimate inundation 12 . In general, the stochastic models are not site specific and therefore do not account for systematic variation in immediate environment of the fault, e.g. change in lithology or seismic wave/rupture interaction due to free surface/fault geometry which could influence the slip distribution over the fault plane 13 . The latter phenomena may explain features like the enhanced shallow slip observed in some recent large tsunamigenic earthquakes 14 , such as Tohoku 15 and Haida Gwaii 16 , or tsunami earthquakes like Mentawai 17 .
We introduce the concept of a slip probability density function (SPDF) to describe the spatial variability of slip in an ensemble of models. We use it to provide a description of the coverage of slip across the fault plane used in hazard analysis. This is particularly important for coastal areas located near subduction zones where variations in the location of peak slip has a large effect on the distribution of the tsunami wave height 6,12,18,19 . The SPDF is not to be confused with the probability function used in some stochastic models that are based on one-point statistics to describe the slip distributions 20 . The SPDF is based on an averaged stack of all the slip distributions in the ensemble (see Methods for a description of how the SPDF is calculated), and is related to the phase spectrum of the slip distributions at large wavelengths. Generally, the phase spectra of the slip distributions is mainly random in stochastic models 12 . Therefore, the stack of such distributions, in large enough ensembles should produce near uniform probability of slip across the fault plane 21,22 with the exception of fault border which is controlled by the tapering filter. This is demonstrated in Fig. 1a where the SPDF has been produced using 500 stochastic slip distributions on a 520 km by 250 km Tohoku-like fault plane.
Studies where the phase of the lowest wavelengths are coherently constrained are generally focusing on the reproduction of particular past events (e.g. Tohoku 19,23 ). This naturally leads to a concentration of the SPDF in the area of the designated asperity. An example of this is given in Fig. 1b where the slip was distributed based on Gaussian function centred in the middle of the fault (i.e. along strike = 269 km, down dip = 125 km, standard deviation = 240 km). This is desirable for recreating some variability around the estimated source for particular historical events. However, its application for describing future earthquakes 12 and consequently associated hazard analysis, should be considered only under the strong assumption that an individual past event has illuminated a persistent asperity (e.g. an highly coupled section of the megathrust) on which earthquakes tend to repeat in a broadly similar way over a long time period.
In this study our aim is to ascertain a SPDF that can be used in tsunami hazard that is specific to a particular fault/area but without the assumption that the next earthquakes will add-up similarly to the last one. Modulating factors deriving from possibly persistent features, such as distribution of highly coupled regions of the interface, Figure 1. Systematic slip probability density function based on different methods of distributing slip on a 520 km by 250 km fault. The x-axis is along fault strike, the y-axis is along dip with the fault surface along the top of the subfigures. 500 simulations were used to produce the SPDF in each subfigure. (a) Using a uniform PDF for the spatial distribution of slip with a boundary taper; (b) a stationary Gaussian function is used as the PDF, the slip taper along the surface has been removed (c) best attempt at producing uniform SPDF that extends to the free surface. Some examples of the slip distributions used to make each of SPDFs are provided below the SPDFs.
can be in principle added afterwards. Here, we use a large ensemble of 2D dynamic simulations to investigate if a Tohoku like fault exhibits depth-dependent systematic deviations from a spatially uniform slip probability, only as a consequence of dynamic effects during the rupture. We test this hypothesis for large tsunamigenic earthquakes by comparing systematic features observed in the dynamic models with the SPDF generated using a generic stochastic model. From this comparison, we derive a transfer function that corrects stochastic slip distributions to be used in SPTHA. As a test case, we analyse the effect of this correction on the conditional probability of exceedance of maximum tsunami wave heights along the Tohoku coastline; the probability of exceedance that we consider is conditional to the occurrence of a M 9 earthquake.

Results
As we are interested in tsunamigenic earthquakes we firstly aim to extend the near uniform SPDF in the stochastic models to the surface in order to allow large slip near the Earth's surface. This is shown in Fig. 1c and it is this SPDF, averaged along strike, with which the dynamic simulations are compared to later in this study.
In the dynamic rupture simulations, a stochastic initial stress distribution is generated by taking the spatial derivative of a 1D slip distribution obtained using the stochastic composite source model discussed earlier. The advantage of this technique is that it allows for the natural concentration of the pre-stress into high stress asperities. The location of nucleation is chosen randomly from locations of high initial stress. A linear slip weakening friction law with a regularisation of the normal stress evolution during rupture 24 was used and the material properties are homogeneous (see Methods and Fig. 2). Near the surface (depth < 3 km), a low shear stress is enforced as we assume that aseismic processes have lowered the shear stress in this zone 25 . There is a transition from aseismic zone to the seismic section of the fault (between 3-6 km depth) which occurs within the wedge. The effective normal stress, σ n , varies as a function of depth based on the difference between the hydrostatic and lithostatic pressure starting from the non-zero value of 0.5 MPa in the trench zone, in order to avoid rupture jumps at the surface due to unrealistic near zero strength. At 25 MPa we assume that the pore pressure tracks the increasing normal stress and the effective normal stress remains constant with depth 25,26 . This choice of frictional parameters produces a 7.5 MPa strength drop in the lower section of the fault, assuming a constant static value for the effective normal stress. Simulations were limited to a homogeneous 2D model with a 1D fault, due to the high computational cost of performing 3D simulations.
The large variability of the initial shear stress (see Fig. 3a,c,e and Supplementary Information Figures S5 and S6) and the location of the nucleation lead to large variations in the size and extent of slip in the 500 dynamic simulations due to the nonlinear behaviour of rupture 27 . The slip profiles on the 1D fault are converted to seismic moment by assuming that the effective along-strike length scales with the mean slip and width (see Methods). Using this scaling relationship the numerical slip distributions cover a range between M w 8.2-9.5 (see Figure S4 in Supplementary Information). Small events (i.e. M w < 5.5) have been omitted from further analysis as the nucleation patch predominantly controls their slip distribution. Ignoring these small events leaves us with 470 slip distributions.  Fig. 3 where the large variation between slip distributions collapses when viewed by magnitude. Figure 3 provides a sample of the shear stress and slip distributions for three magnitude bins while all 0.2 magnitude bins from 8.4 to 9.6 are plotted in Supplementary Information (i.e. Figures S4 and S5). In all magnitude bins there are instances of surface rupture, we discount the M 8.2-8.4 and M 8.4-8.6 slip distributions from this analysis as there are less than 20 events in these samples and are therefore not representative. The probability of an earthquake generating surface rupture increases with magnitude (e.g. 13% of earthquakes in the M 8.6-8.8 bin ruptured to the surface compared with 82.1% for the M 9.0-9.2, see Table S1 in the Supplementary Information). With increasing surface rupture the spatial distribution of slip becomes more asymmetrical compared with slip distributions that do not reach the surface. In the cases where rupture does reach the surface, the slip distribution is asymmetric with larger slip near the surface due seismic wave interaction with rupture which causes dynamic reduction of the normal stress and larger slip near the surface 13,28 .
In Fig. 4, the location of maximum slip in all cases occurs at or below a distance of 50 km down dip, this is due to the transition from the seismic (i.e. 59 km down dip, 6 km depth) to the aseismic section (27 km down dip, 3 km depth) the point at which the shear stress systematically decreases below the dynamic threshold. Changing the depth range of the aseismic zone, altering the level of shear stress in it, or using another frictional parameter to represent it (e.g. a rate strengthening zone or increasing d c ) could alter the location and size of peak slip. The maximum slip is distributed over the seismic section of the fault for lower magnitude events, i.e. 50 km-200 km down dip for the M 8.4-8.8 events. With increasing magnitude this spatial range decreases, for M 9.2-9.4 events the maximum slip is limited to occur between 50 km-120 km down dip.
In the simulations this increasing constraint of maximum slip with magnitude is due to the increasing influence of surface rupture; larger earthquakes are more likely to rupture to the surface which in turn produce asymmetric slip distributions.
Comparison between dynamic and stochastic models. Our aim is to compare the stochastic slip distributions with the slip distributions generated by the dynamic simulations. For the stochastic simulations we adopt the common working assumption, used in most stochastic models 9 , of a uniform slip density function that tapers to zero at the fault boundary for all magnitudes. However, Fig. 3 demonstrates that the dynamic slip distributions exhibit different systematic behaviour for different magnitudes and whether they reach the surface or not. The dynamic slip distributions are compared with the stochastic slip models in magnitude bins of width 0.2 between M w 8.4 to 9.4. The ranges 8.2-8.4 and 9.4-9.6 have been excluded due to the low number of events in these bins (i.e. less than 10). The slip distributions were further subdivided into those that reached the surface and those that did not. The SPDF for each bin was calculated using a 1D version of Eqn 2 in the Methods and are displayed in Fig. 5. In all bins the SPDF tends towards zero at the bottom of the fault due to the large slip weakening distance, d c , (see Methods and Figure S1 in the Supplementary Information) used at depth in the dynamic simulations.
The 2D SPDF in Fig. 1c was averaged along strike to produce a 1D depth dependent stochastic SPDF (see black line in Fig. 5) in order to compare with the SPDFs generated using the dynamic simulations. The shape of the stochastic SPDF is assumed to be similar over all magnitude ranges (i.e. near uniform across the fault). Comparing the different SPDFs (Fig. 5), the variation in amplitude is due to the spatial concentration of slip for a given bin as the integral of each SPDF across the fault plane is 1. Therefore, comparison of the SPDFs provides a means of comparing the concentration of slip between the different subgroups. As slip scales with magnitude, when comparing the amplitude of the SPDF between different magnitude bins, larger amplitude does not imply larger slip but rather relative concentration. In cases where comparing SPDF in the same magnitude bin, then amplitude of the SPDF can be viewed as a proxy for the relative slip location.  Comparing the dynamic SPDFs with the stochastic source model (black line) Fig. 5 demonstrates that the stochastic model systematically underrepresents the concentration of slip near the surface (i.e. < 20 km depth) for all magnitude bins. For M < 9.2 each magnitude bin contains two types of rupture: those that reached the surface and those that do not. These two classes of rupture produce different slip distributions: ruptures that reach the surface produce asymmetrical slip distributions with the maximum slip located closer to the surface; earthquakes that do not rupture all the way to the surface produce a more symmetric slip distribution. This bifurcation is related to the point at which the fault is producing ruptures that can penetrate the low shear stress zone near the surface 29 and dynamic-free surface effects (i.e. increased normal stress prior to slip due to the reflection of seismic waves onto the fault 13,28 ). The spatial segregation of the symmetry of the slip distribution compliments the concept of depth dependent failure domains where asperities dominate the fault plane at depth (i.e. > 20 km deep) while earthquakes that breach the surface contain strong effect from the free surface boundary as well as the aseismic zone 30 . In the ensemble, some of the M < 9 events that rupture to the surface may be similar to tsunami earthquakes 31 , where surface rupture occurs and slip is primarily located under or near the wedge with near-trench nucleation 27 . Above M 9.2 all earthquakes rupture to the surface producing a similar asymmetric slip distribution. The likelihood of surface rupture in the ensemble of simulations may be altered by the initial conditions used in the numerical model. For example changing the normal stress or the difference between static (μ s ) and dynamic (μ d ) friction coefficients which would lead to variations in the stress drop; shrinking the width of the aseismic zone, by decreasing the negative stress drop in it or by varying the means of reproducing an aseismic zone (e.g., increasing d c ) may all affect the extent of rupture. A test to evaluate the effect of introducing a more elastically compliant wedge is presented in the Supplementary Information (see Figure S7) where the final slip distributions between a model with homogeneous material (i.e. the same used in the ensemble described above) and that containing a wedge with a lower v p , v s and density (i.e. 4.7 km/s. 2.1 km/s. 2.5 kg/m 3 respectively). In the sample study, it was found that the wedge has a bigger effect on larger earthquakes than smaller ones in general causing larger amounts of slip. However, the general shape of the final slip distribution remained similar. A more comprehensive study on the effect of the wedge is beyond the scope of this study.
Additionally, features in 1D simulations may not be present in 2D simulations as rupture could propagate around potential barriers. Therefore, analysis of different dynamic conditions in the aseismic zone as well as comparisons between 1D and 2D simulations should be performed in future studies.
In order to produce a stochastic source model that better represents the systematic dynamic features depicted in Fig. 5 the stochastic methodology requires some modification. We introduce a depth dependent transfer function, Λ x ( ), representing the differences between the stochastic and dynamic SPDFs (displayed in Fig. 5 and described in the Methods section). Λ x ( ) is dependent on the magnitude and whether there is surface rupture, requiring the function to be changed based on the size of the stochastic earthquake and the probability of surface rupture for a particular magnitude size. An example of the application of the 1D Λ x ( ) as a depth dependent function to a 2D stochastic model is provided in Fig. 6. All the Λ x ( ) are provided in Figure S8 in the Supplementary Information.
Slip probability on the fault and exceedance probability for maximum tsunami wave height. The dynamic simulations have demonstrated that the slip probability density function is non-uniform and varies with magnitude; these are important features for hazard that are not generally considered.
The 'traditional' stochastic source models were produced using the same method that generated the SPDF displayed in Fig. 1c (i.e. not tapering the slip at the surface). The stochastic slip distributions are then multiplied by the transfer function Λ x ( ) generated using dynamic simulations in the 9.0-9.2 magnitude range. Two transfer functions were used: one in case where there is surface rupture and one where there is not (represented by the bright red curves (solid and dashed) in Figure S8). The choice of transfer function is taken based on the probability of surface rupture occurring in the numerical simulations (see Table S1). 500 slip distributions were generated using the 'traditional' stochastic source model and a further 500 using the modified stochastic source model, in both cases the slip distributions were magnitude 9.0 events. The importance of applying such a correction to the traditional slip distribution is shown by the SPDF on the Tohoku fault plane (Fig. 7a), constructed by stacking the modified slip distribution. This SPDF shows an increase of probability For gauging the effect this has on tsunami hazard, we compare how a tsunami hazard metric such as the maximum tsunami wave height (which is the peak of the tsunami waveform from unperturbed sea level at 1 m depth), H max, is affected by our proposed approach. As the dynamic simulations were based on a Tohoku like fault, the eastern Japanese coastline provides an appropriate location for examining variations in H max . For each earthquake, the corresponding slip distribution has been mapped on the 3D Tohoku fault plane subdivided into 398 subfaults, and the H max has been calculated by combining the slip with pre-computed tsunami Green's functions (further details in Methods section).
From each of the 500 magnitude 9.0 slip distributions the tsunami was simulated for both stochastic and corrected model types, thus producing a robust sampling of a wide variety in H max along the coastline in both cases. Figure 7b-c display the probability of exceedance of H max at each receiver for both ensembles. For assessing SPTHA, these probabilities should be combined with those of the earthquake occurrence 2,9 . The original stochastic model produced a median H max not greater than 20 m between 35°-41°N in comparison to the median of modified model reached 30 m in the same area (see Fig. 7b,c). In the ensemble of the modified stochastic slip models there is a subset of events producing very high H max values in several areas, namely between 39° and 40°N and around 37°N. In particular, the occurrence of + 60 m H max values are very infrequent being present in at most 5% of the ensemble.
The inset in Fig. 7a shows the H max hazard curves aggregated between 36°-40°N as we have only considered the Tohoku section of the fault system; outside of this latitude range other sections of the fault system not considered in this study may become more influential in determining the tsunami hazard. For the Tohoku section, the hazard curves in Fig. 7a differ if using traditional or modified stochastic models, with the former resulting in an underestimation of the hazard for large tsunami intensities (H max ). Therefore, the modification with Λ x ( ) produces the largest, H max values which are missed in generic stochastic source models.
The maximum wave heights and runups measurements for the 2011 Tohoku tsunami and for two historical tsunamis in the region are also plotted for comparison (Fig. 7d). Since very large earthquakes and tsunamis are infrequent, the available record of past events is likely too short to be complete 33 ; however, the events in such catalogues must be foreseen in the conditional probability of exceedance for a M9 earthquake assessed from numerical modelling for future earthquakes and tsunamis. Comparison of Fig. 7b-d indicates that some of the extreme values could not be reproduced with the traditional approach, while the present approach forecasts extremely low probability wave heights exceeding the historical observations at several places.
Most features between the two models are the same (i.e., a shear modulus of 30 GPa and stress drop per subevent of 1 MPa). Sensitivity analysis of the transfer function to these parameters and others (e.g., the type of friction law, heterogeneous media, strength of the normal stress, asperity size, variation in the length of the aseismic zones and fault geometry, etc.) is required in order to better sample the range of possible ruptures. The same holds for the scaling law used to calculate the rupture length (discussed in more detail in the Methods). However, a preliminary analysis presented in the Supplementary Information (Figures S12, S13 and S14) shows that using the M w 8.8-9.0 and M w 9.2-9.4 transfer functions produces a similar broad scale pattern to that presented in Fig. 7 using the M w 9.0-9.2 transfer function.
Therefore, while there is uncertainty in the calculation of the magnitude, due to the use of an empirical length scaling, the conditional probability of exceedance along the coastline produces robust features observed in all 3 ensembles (see Figure S14). However the 1D to 2D conversion is still an important feature that requires further investigation prior to real application in SPTHA, particularly in ascertaining the probability of surface rupture occurring for a given magnitude.

Discussion
We have introduced the concept of the SPDF and have highlighted the importance of accounting for its spatial variation when considering earthquake source models in SPTHA. We have proposed a new method based on the application of a transfer function for altering the stochastic source model according to systematic features observed in an ensemble of dynamic earthquake simulations. This allows the generation of rapidly computed slip models based on the "dynamic" SPDF, which accounts for systematic features, such as shallow slip amplification in mega-thrust earthquakes.
Taking the Tohoku fault as a case study, we computed 500 simulations using a simplified 2D dynamic model with an isotropic medium, a linear slip weakening friction law with a shallow low stress zone as an approximation for the aseismic zone. Uncertainty in the initial shear stress distribution and nucleation location was accounted for using a different stochastic shear stress distribution in each simulation and a randomly located nucleation patch. The simulations produced events that ranged in size from M w 8.4 to M w 9.5, with a range of diverse characteristics that collapse to distinct slip distributions when grouped in 0.2 magnitude bins and whether they produce surface rupture or not. Magnitude dependent, spatially heterogeneous dynamic SPDFs that are dependent on whether there is surface rupture or not were then generated and compared with the standard, near uniform stochastic SPDF. The comparison highlighted that the spatial variation and magnitude dependence of dynamic SPDFs in conjunction with the probability of rupture reaching the surface are not accounted for in generic stochastic source methods. For example stochastic slip distributions consistently underestimate the concentration of slip in a near surface zone that extends down to a depth of 20 km in cases when rupture reaches the surface.
Using a transfer function based on this difference, we corrected the stochastic source models and then we calculated the probability of exceedance for H Max along the eastern Japanese coastline for 500 M w 9 Tohoku-like events using both original and corrected slip distributions. We establish that the modified sources produce more extreme (and low-probability) H Max values. This demonstrates the importance of incorporating systematic features observed in dynamic simulations in SPTHA.
The extensive extreme impact of the Tohoku 2011 tsunami was considered somehow unforeseen, even if data from past events might have been perhaps more carefully considered 34 . Using our or other techniques in SPTHA shows the Tohoku event was a combination of two low probability events; an M9 earthquake (annual probability is roughly 10 −3 35,36 in the Tohoku area), and a large amount of shallow slip (Fig. 7a). This probability is conditional as only a M 9 earthquake is considered in the Tohoku region; the overall H max probability would be different were all other magnitude ranges considered and the fault plane extended along the whole subduction zone. The use of hybrid schemes such as the one we propose offers a computationally affordable means of including important dynamic features that may otherwise be overlooked in more generic earthquake source models.

Methods
Stochastic Modelling. We use a composite stochastic source model 37 to produce the stochastic kinematic slip distributions. This involves placing a hierarchical set of circular sub-events on the fault plane. The subevents are distributed based on a power law size-number relationship. The subevents are allowed to overlap with each other, when this occurs the slip from the different subevents at the particular location are added together. The number, n, of sub-events of radius R is: where D is the fractal dimension and p is a constant. In this study, D has been set to 2 in order to generate self-similar slip distributions 38 . R is the radius of the subevent and is bounded within the range [R min , R max ]; p is determined based on the seismic moment of the earthquake, the fractal dimension and the stress drop 39 . The slip distribution across the individual subevents is described by the Eshelby slip function 40,41 . The distribution of Scientific RepoRts | 6:35007 | DOI: 10.1038/srep35007 the subevents across the fault plane is based on a probability density function (PDF) which is usually a uniform function 39,41 . The SPDF in Fig. 1a was generated using a uniform PDF to describe the distribution of subevents coupled with the tapering of slip at the fault boundary. In Fig. 1b a fixed single Gaussian function centred at 269 km along strike and 125 km downdip is used to describe the distribution of subevents. For Fig. 1c the width of the fault plane was artificially doubled. The PDF that describes the distribution of subevents is composed by one or more Gaussian functions whose centres are randomly allocated to lie within one half of the fault plane (i.e. within the original width of the fault). After the subevents are distributed on the fault and then summed, the target fault is cut from the half of the plane where the Gaussian functions were centred. Slip is tapered on the other three fault boundaries. Slip is then normalised to the required moment.
Spatial slip probability density function (SPDF). The SPDF, Δ K , is based on the mean SPDF for an ensemble of models used, for example, in hazard analysis. Each slip distribution in the ensemble is converted into a probability density function by dividing the slip at each location on the fault plane by the total cumulative slip for that particular distribution (this is represented by the function in the brackets in Eqn 2). The SPDF is then calculated by taking the mean from the ensemble of PDFs across the whole fault plane. This can be expressed as: where N is number of models in the ensemble, δ x y ( , ) i K is the stochastic slip at position x,y on the fault plane in i th model. The superscript K stands for kinematic and is used to differentiate stochastic slip distributions and SPDFs from similar slip distributions and functions generated from dynamic simulations that contain the superscript D. The denominator in the Eqn 2 is the integral of slip for the i th model where A represents the area on the fault. Integrating the function Δ K over the whole fault gives a value of 1. The terms inside the brackets represent the slip probability density function for the individual stochastic slip models.
Dynamic Modelling. The 2D wave equation was solved using a spectral element method (SPEC2D) with a linear slip weakening friction law 42 with a temporal smoothening of the normal stress change on the fault plane 24 . A Tohoku like fault trace is used with the curved fault geometry based on Slab 1.0 fault trace that extends to the trench 43 (see Fig. 2). The material parameters are homogeneous in this cross-sectional model and are set from commonly observed seismic velocity and density values in the crust (v p = 5.98 km/s, v s = 3.2 km/s, ρ = 3000 kg/m 3 ); no water layer is present as it has been shown that wave propagation in this layer has negligible effects on the final slip 29 . For the friction law, the slip weakening distance is d c = 1 m; the static (μ s ) and dynamic (μ d ) coefficient of friction are 0.6 and 0.4 respectively; the reference velocity (ν*) was set to 0 m/s and reference slip (δ σ ) to 0.1 m for the temporal smoothening of σ n . At a depth of 51 km the d c linearly increases with down dip distance to a value of 300 m at the bottom of the fault, this is done to assist the rupture arrest with the assumption that the fault is becoming aseismic at depth (see Figure S1 in the Supplementary Information).
The shear stress distribution, τ(x), is taken from different stochastic shear stress distributions. These distributions are constructed using the composite source model discussed earlier where the subevents are distributed onto a 2D fault plane based on one Gaussian PDF with a standard deviation of 32 km with a centre that is randomly chosen over a depth range of 2-45 km. The subsequent stochastic slip distribution is converted into a stress distribution by taking the spatial derivative of the slip 10,38 . The 2D shear stress distribution is then averaged along strike in order to generate a 1D, depth dependent shear stress distribution which has a k −1 spectra if there is no surface rupture (see Supplementary Information). A k −1 stress spectra is used rather than a 1D self similar k −0.5 44 as the aim is to produce slip distributions that are representative of depth dependent features in 2D slip distributions. The spectra of the slip distributions generated in the dynamic simulations have a k −2 spectra (see Figure S2 in the Supplementary Information) which is consistent with the spectra of the stochastic models to which it is being compared with.
The element size and grid distance vary along the fault plane in the dynamic model requiring that the stochastic shear stress has been linearly interpolated onto the numerical grid points. The random location of the Gaussian PDF leads to a variation of location of the high stress asperity between simulations as shown in Fig. 2b- In the wedge zone (depth < 6 km) we assume that aseismic processes have lowered the shear stress 25 . Consequently the amplitude of the shear stress in the near surface zone (< 3 km) is scaled such that it is less than or equal to 0.5μ d σ n . For the rest of the fault the shear stress is scaled to μ s σ n with a linear transition over a depth of 3 km between the two scaling regimes. Nucleation on the fault is produced by lowering σ n smoothly using a Gaussian function with an amplitude of 3.3 MPa and a standard deviation of 4.2 km and example is provided in Fig. 2f. The size of patch was tested (see Supplementary Figure S3) and was found not to have a dominant effect on the final slip distribution. The location of this patch is randomly chosen (see Fig. 2c) with the only constrain being that the shear stress is high (i.e. τ(x) ≥ 0.9μ s σ n ). A 19-point moving average has been applied to the resulting slip in order to smoothen small jumps in the slip caused by small changes in the normal vector between adjacent cells along the fault plane.
Conversion of 1D slip profiles to 2D seismic moment. The 1D slip profiles are converted to 2D seismic moment by assuming that the effective along-strike length L scales with the mean slip and width based on the empirical relationship 45 Δσ is the average stress drop across the whole 2d slip distribution and was taken to be 7.5 MPa based on the choice of frictional parameters in the dynamic model and produces slip profiles that are comparable to inversion observations for Tohoku (see Fig. 3d,f). W is the width that is set to the rupture size in the individual simulations, and δ′ is the mean slip. Using this approximation, Figure S1 shows that the earthquakes magnitude varies between M w 7.8-9.6 with very small events (i.e. M < 5) being omitted as the nucleation patch predominantly controls their slip distribution. Varying the length/width scaling relationship would shift the corresponding magnitudes as well as their range. For example, assuming that rupture is square for all sizes produces a lower magnitude range of M w 8.2-9.2 (see Figure S9 Supplementary Information). Altering the average stress drop in Eqn 3 (or an equivalent scaling relationship) leads to a uniform shift in the magnitudes.
Modifying the stochastic slip to account for observations in simulations. The aim is to reproduce the magnitude dependent SPDF observed in the dynamic simulations in the stochastic slip distributions which can be produced much more rapidly. To do so we describe the SPDF generated by the dynamic simulations, Δ D (x), in terms of the stochastic SPDF, Δ K (x), using a transfer function, Λ x ( ): The transfer function describes the average difference between the dynamic and stochastic SPDFs, and is determined by dividing the two SPDF from each other at each point along the fault plane: This transfer function gives a value of < 1 when the stochastic SPDF consistently produces a larger amount of slip than the ensemble of dynamic slip distributions (an example of this is the bottom of the fault where a high d c is used to represent an aseismic zone). When Λ x ( ) is greater than one, the stochastic SPDF underestimates slip relative to the dynamic ensemble (e.g. near the surface for + M8.6 events). In Fig. 5b fifteen Λ x ( ) have been calculated based on grouping the slip distributions by 0.2 magnitude bins and whether they produced surface rupture or not.
In order to generate individual stochastic slip distributions that account for general features observed in the dynamic simulations, Eqn. 4 is rewritten where ∆ x ( ) K is replaced with a 1D version of Eqn 2, producing: inside of the brackets relates to individual slip distributions and by multiplying the individual stochastic slip distributions by Λ x ( ), the ensemble of slip distributions, will approach Δ D (x) as N increases. Therefore, to generate the modified stochastic source model we multiply the stochastic slip distributions by Λ x ( ), the slip distribution is then normalised in order to produce a predefined slip that correlates to a predefined magnitude. Additionally, Λ x ( ) is magnitude dependent requiring the function to be changed based on the size of the stochastic earthquake required. The SPDF and consequently the Λ x ( ) also varies depending on whether rupture reaches the surface or not. To incorporate this complexity in the stochastic slip models two transfer functions are generated: one based solely on events that reach the surface, the second contains all other events. The different transfer functions are given in Figure S8 in the Supplementary Information. In terms of practical application, the Λ x ( ) function is simply multiplied to the slip distribution produced using the technique described in Stochastic Modelling Section (see also refs 39, 41) leading to the generation of a 'modified stochastic source model' . The choice of transfer function is taken based on the probability of surface rupture occurring in a given magnitude bin in the dynamic simulations (see Table S1). In the case of 2D slip distributions we assume that Λ x ( ) is a depth dependent function this means the two dimensional slip distribution is multiplied by Λ x ( ) which remains constant along strike (i.e. δ x ( ) i K is replaced with δ x y ( , ) i K in Eqn 6).
Tsunami numerical simulations. The subduction interface geometry in the Tohoku region is constrained by Slab1.0 model 43 15 ). The vertical coseismic sea floor displacement associated to each slip distribution (i.e. the initial condition for the tsunami propagation) is obtained as a combination of the single displacement fields arising from each subfault, numerically computed using the commercial software Abaqus version 6.9 (www.simulia.com, date of access: 16/06/2016) and considering the medium as elastically homogenous; in particular, we also include the contribution of the horizontal deformation in the region of steep bathymetric slopes 46 and we apply a two-dimensional filter 47 to each field in order to take into account the attenuation of the sea floor deformation through the water column. Tsunami numerical modelling for each subfault (Green's functions) is performed by using the HySEA 48,49 code that solves the non-linear shallow water equations using a hybrid numerical scheme (Finite Difference two-step scheme similar to leap-frog for the propagation phase in open sea combined to a second-order Finite Volume TVD-Weighted Average Flux scheme for the inundation step). The bathymetric model used for the tsunami propagation is SRTM30+ (http://topex.ucsd.edu/WWW_html/srtm30_plus.html) and the spatial resolution of the computational grid is 30 arc-sec. We collect the waveforms at 2579 receivers located along the 50 m isobath off the eastern Japanese coast (black dots in Fig. 7a). For each slip distribution, the resulting tsunami waveform at each receiver is obtained by linearly combining the Green's functions; we extract at each receiver the maximum wave height to get H max 50 profiles which were then extrapolated in front to the coastline (1 m depth) using the Green's law: