Stress, rigidity and sediment strength control megathrust earthquake and tsunami dynamics

Megathrust faults host the largest earthquakes on Earth, which can trigger cascading hazards such as devastating tsunamis. Determining characteristics that control subduction-zone earthquake and tsunami dynamics is critical to mitigate megathrust hazards but is impeded by structural complexity, large spatio-temporal scales and scarce or asymmetric instrumental coverage. Here we use high-performance computing multi-physics simulations to show that tsunami genesis and earthquake dynamics are controlled by along-arc variability in regional tectonic stresses together with depth-dependent variations in rigidity and yield strength of near-fault sediments. We aim to identify dominant regional factors controlling megathrust hazards. To this end, we demonstrate how to unify and verify the required initial conditions for geometrically complex, multi-physics earthquake–tsunami modelling from interdisciplinary geophysical observations. We present large-scale computational models of the 2004 Sumatra–Andaman earthquake and Indian Ocean tsunami that reconcile near- and far-field seismic, geodetic, geological, and tsunami observations and reveal tsunamigenic trade-offs between slip to the trench, splay faulting and bulk yielding of the accretionary wedge. Our computational capabilities render possible the incorporation of present and emerging high-resolution observations into dynamic-rupture-tsunami models and will be applicable to other large megathrust earthquakes. Our findings highlight the importance of regional-scale structural heterogeneity to decipher megathrust hazards. Tsunamis generated by megathrust earthquakes are controlled by regional-scale structural heterogeneity, according to numerical modelling based on the 2004 Indian Ocean earthquake and tsunami.

V ariations in megathrust earthquake rupture behaviour are associated with tectonic, mechanical and structural factors, highlighting the importance of depth-dependent and along-arc subduction-zone heterogeneity [1][2][3][4][5] . Large tsunamis may be caused by various coseismic mechanisms, including large slip to the trench, as observed during the 2011 Tohoku earthquake 6 and inferred from tsunami deposits along the Kuril trench 7 , slip on splay faults 8 and inelastic coseismic deformation generating substantial uplift near the trench 9,10 .
Here, we hypothesize that a few regional-scale subduction-zone characteristics control earthquake kinematics, earthquake dynamics, and tsunami genesis. These characteristics include rock rigidity and yield strength, plate convergence rate, and regional tectonic stresses. Such features and their heterogeneities can be readily constrained by seismic, geodetic, and geological observations (for example, ref. 11 ).
The 2004 moment magnitude (M w ) 9.1-9.3 Sumatra-Andaman earthquake and subsequent Indian Ocean tsunami initiated a new era of geophysical studies focused on deciphering the complex source characteristics of large subduction-zone earthquakes and how these relate to tsunami generation. The earthquake's unexpectedly large moment magnitude challenged the idea that convergence rate and age of the oceanic lithosphere control the largest earthquakes on subduction interfaces [12][13][14] . The Andaman trench involves a middle-aged oceanic lithosphere, aged 50-70 million years 15 , converging at relatively slow rate as low as 20 mm yr -1 in the northern Andaman section 16 , and historically it hosted only thrust earthquakes of less than M 7.9 (Fig. 1a). Then the 2004, dynamically complex rupture ripped 1,300-1,500 km of faults along the trench for 8-10 minutes and challenged data-driven analyses. Inversion of seismic and geodetic data [17][18][19] is complemented by less conventional data inferences, including normal modes 20 , multiple centroid-moment-tensor inversions 21 , gravity 22 , tsunami 23 or hydroacoustic waves 24 , and seismic arrays 25 . Nonetheless, data-driven source models include large variations 26 and produce synthetic tsunamis in varying agreement with satellite or inundation records 27 . Furthermore, the earthquake generated a tsunami up to 30 m high in northern Sumatra 28 , which seems at odds with slip on the shallowly dipping subduction interface. Slip on several splay faults during or soon after the megathrust event was inferred from seafloor mapping [29][30][31] , tsunami modelling 32 and aftershock relocation 33 . Alternatively, near-trench slip in lithified sediments at the southern end of the rupture may be responsible for the tsunami [34][35][36][37][38] .
Rising computational resources and efficiency now allow data-integrated dynamic-rupture modelling of crustal earthquakes that complements data-driven analyses. Dynamic-rupture models provide physically self-consistent descriptions of how faults yield and slide, while their nonlinearity and complexity limit the total number of feasible numerical experiments. These methods have unravelled complex and/or poorly instrumented events in various tectonic contexts [39][40][41][42][43][44] . However, constraining the dynamics and extracting hazard implications from three-dimensional (3D) dynamic models of megathrust earthquakes requires new methods that assimilate distinct data and inferences.
In this article, we introduce a stringent forward modelling approach to evaluate the key factors controlling earthquake-tsunami behaviour. Large-scale dynamic-rupture modelling provides a bridge between scales harnessing the ability to uncover the physical mechanisms and parameters relevant to the coseismic rupture of naturally complex fault networks comprising megathrusts. We connect subduction mechanics with a scenario of the 2004 Sumatra-Andaman earthquake and tsunami produced by a physics-based computational model. Evaluation against independent observations in the near and far fields ensures that the scenario captures the key Stress, rigidity and sediment strength control megathrust earthquake and tsunami dynamics observed earthquake characteristics (magnitude, moment release, focal mechanism(s), rupture duration, teleseismic waveforms and ground displacements) and main characteristics of the Indian Ocean tsunami. The model reveals that the interplay of a curved megathrust geometry with the along-arc variability of regional tectonic loading and plate convergence rates, and strong rigidity variations with depth near the megathrust interface, are dominant factors controlling the earthquake. We then analyse a major conundrum 45 of earthquake-tsunami interaction that is ill constrained by observational data: the interplay of near-trench earthquake dynamics and sedimentary yielding controlling tsunami generation and propagation. We present two additional scenarios with alternative near-trench sedimentary rock strength. These demonstrate that the strength of trench sediments and the associated changes to on-fault megathrust and splay fault slip versus off-fault plastic yielding control the seafloor displacements and, in turn, the tsunami.
To model earthquake dynamics, choices about the preexisting state of stress and the friction law governing fault strength and sliding, as well as the lithological structure and fault geometries, are required. We constrain all required initial conditions on subduction-zone scale and specific to the Sumatra-Andaman trench. Regional principal stress orientations from stress inversion are aligned with geodetic observations and combined with seismo-tectonic observations and the Mohr-Coulomb failure theory to constrain principal prestress magnitudes and fault friction parameters ( Fig. 1a and 'Principal stress orientations and magnitudes' in Methods). We account for along-arc variations in convergence rate and assume over-pressurized fault fluids (Fig. 1b). Without regional prestress variation, rupture fails to stop dynamically at the northernmost megathrust, and slip larger than 40 m accumulates there (Extended Data Fig. 1a and 'The importance of convergence-rate-modulated regional driving stresses' in Methods). Horizontal ground displacements then exceed observations by factors of 2 to 5 (Extended Data Fig. 1b). Inferences of principal stress rotations after megathrust earthquakes provide direct information about the ratio of stress drop to prestress 4 . We propose to constrain the coseismic frictional strength drop on the megathrust on the basis of stress rotation observations ( Fig. 1c and 'Coseismic stress rotation to constrain fault friction' in Methods). The 3D structural model ( Fig. 1d and 'Regional lithological structure' and 'Fault geometry' in Methods) accounts for strong rigidity variations with depth in the near-fault region ( Fig. 1e and 'Regional lithological structure' in Methods), bathymetry and topography, and a complex, embedded fault network consisting of the curved megathrust and three splay faults, a backthrust and two forethrusts splaying off of it. All faults intersect each other and the ocean floor, with megathrustbathymetry intersection angles as low as 30° (ref. 46 ). We find that the following initial conditions are required to produce a realistic earthquake-tsunami scenario of the 2004 event: a non-Andersonian prestress state with pronounced along-arc and depth-dependent variations, near-lithostatic fluid pressure that results in low effective normal stress near the megathrust, strong frictional weakening, a  68 . Arrows compare the trends of the maximum horizontal stress (SH max ) used in this scenario (blue) with those by ref. 68 (black). Dashed black lines bound the region where the trend of σ 1 linearly increases from 309° to the south to 330° to the north in this scenario. Cross-sections show the plunge of σ 1 (blue) and the minimum principal stress (σ 3 , red) used in this scenario, compared with those in ref. 68 (black). red (radius 10 km) and black (radius 60 km) circles, centred at the hypocentre, bound the region where the slip weakening distance (D c ) increases linearly from 1.0 to 2.5 m. b, Modulation of the relative prestress ratio R 0 by the magnitude of the convergence rate inferred from rigid plate tectonics considerations, using the Euler pole inferred by ref. 69 (grey dot). c, rotation of the maximum, σ 1 (grey to blue), and minimum, σ 3 (grey to red), principal stress axes due to an earthquake (data from ref. 68 ). d, Schematic view of the here adapted regional 3D structural model. e, Strong rigidity variations with depth near the megathrust interface compared with global averages L2012 9 and Sr2019 5 . For the latter comparison, we assume 4 km water depth.
3D velocity structure specifically accounting for the low rigidities at shallow depth within the subduction channel ('The importance of strongly depth-dependent rigidity' in Methods and Extended Data Fig. 2), and sediments of high yield strength.
Regional variation of tectonic driving stresses modulated by convergence rates and the geometry of the slab and splay faults explains first-order static earthquake characteristics, such as geodetic displacements at regional distance. Along-depth variation in rock rigidity and yield strength of near-trench sediments determine rupture dynamics, such as rupture speed, as well as coseismic uplift and, thus, tsunami generation. Figure 2a shows the rupture evolution of the observationally constrained dynamic earthquake scenario of the 2004 Sumatra event. Complex rupture dynamics evolve across 1,300 km of the fault   Boomerang-shaped megathrust slip pulse consisting of multiple rupture fronts. Snapshots of absolute fault particle velocity (slip rate). b, Moment-tensor representation of the dynamic-rupture scenario (blue, large) compared with the moment-tensor representation inferred by the USGS (black, large) and comparison of a five-point source model derived from the dynamic scenario (blue, small) with a respective solution from teleseismic inversion 21 (black, small). c, Synthetic moment rate release compared with observational inferences from teleseismics by ref. 17 (II from intermediate-period surface wave and long-period seismograms, III from teleseismic body waves, intermediate-period three-component regional waves and long-period teleseismic waves) and ref. 47 . d, Along-arc and along-depth variation of rupture speed, compared with inference from rayleigh waves 48 and from acoustic observations 49 . e,f, Comparison of synthetic ground displacements (blue) and geodetic observations (orange). e, Horizontal ground displacements. f, Vertical ground displacements. The complete uplift and subsidence resulting from the dynamic-rupture scenario (in m) is shown in red to blue, including noticeable splay faulting signatures.
system, including dynamic-rupture transfers to both backthrusts and the forethrust. Along-depth variation in rock rigidity leads to faster rupture at depth and slower rupture at the trench, yielding the boomerang-shaped slip pulse apparent in Fig. 2a. This pulse consists of multiple rupture fronts, which are caused by reflected waves and head waves generated at structural interfaces and the complex free surface.
Along-arc variations in convergence rates do substantially modulate earthquake evolution. These translate into modulations of stress drop and fault slip. In the hypocentral region, the subduction interface is nearly optimally oriented. While rupture is forced to initiate (Supplementary Information, section Earthquake nucleation), locally decreased prestress reproduces spontaneous arrest immediately to the south of the hypocentre. The rupture also arrests spontaneously to the north, where lower convergence rates yield a critically lower prestress.
Our modelling insights are strengthened by the scenario's agreement with the main source characteristics of the 2004 earthquake. Its M w 9.2 falls in the range of inferred estimates (M w 9.1-9.3 (ref. 26 )). The synthetic focal mechanism, with strike/dip/rake angles of 343°/10°/114° (Fig. 2b) agrees with the respective observational inferences of 336°/7°/114° by the US Geological Survey (USGS). The synthetic moment rate release resembles available source time functions inferred from teleseismic data 17,47 (Fig. 2c). The scenario is also in excellent agreement with all five focal mechanisms of the multi-centroid-moment-tensor teleseismic inversion by ref. 21 (Fig. 2b). In Supplementary Information, section Multiple point source analysis: source time functions, we detail the derivation of equivalent multiple moment tensors and source time functions from dynamic-rupture models. The scenario also reproduces the timing, amplitudes and envelopes of recorded long-period teleseismic waveforms (Extended Data Fig. 3a at nine stations (Extended Data Fig. 3b) at periods of 150 s to 500 s with an average relative root-mean-square misfit (rRMS, as defined in ref. 44 , equation 6) of 0.7 ('Teleseismic model validation' in Methods)). Figure 3 shows the rupture speed and slip characteristics along the megathrust from the dynamic-rupture scenario. The speed of the modelled rupture front is slower than the S-wave speed of the host rocks ( Fig. 3a) everywhere. Shallow supershear rupture is prevented by off-fault plastic yielding. Rupture speed strongly varies along depth due to the incorporated rigidity variations (Fig.  1e). The average speed resembles observational Rayleigh wave 48 and acoustic 49 inferences including along-arc variations (Fig. 2d). The scenario captures changes in rake along megathrust strike, from thrusting in the south to increasingly oblique faulting farther north (Fig. 3b). This variation follows the change in the strike of the non-planar megathrust relative to the laterally varying regional loading.
The dynamic-rupture scenario yields a smooth fault slip distribution with along-arc variations (Fig. 3c). Regions of high slip at the trench alternate with flatter portions of the shallow interface that experience less slip. Synthetic fault slip agrees in terms of along-arc and along-depth distribution to first order with data-driven kinematic models (Fig. 3d,e).
Near-field and far-field horizontal and vertical ground displacements reproduce observations (Fig. 2e,f and Extended Data Fig. 4) in terms of amplitude and orientation. We denote deviations to spatial and temporal model limitations. The dynamic scenario does not account for Earth curvature, small-scale 3D crustal structure in the far field or the contribution of after-slip. Near-field synthetics may be affected by the simplistic treatment of earthquake nucleation and arrest in our model.
All three splay faults incorporated in the model are dynamically activated (Fig. 2a) and cause a noticeable signature in uplift and subsidence (Fig. 2f). The splays host substantially less slip (6-8 m) than the subduction interface (up to 25 m); however, slip efficiently transfers into vertical uplift due to their steeper dip. The total contribution of splay fault slip to vertical uplift remains limited in spatial extent and amplitude.
While outer-rise aftershock sequences help to constrain the extent of slip to the trench during megathrust ruptures 50 , the balance between shallow slip and off-fault yielding remains elusive from geodetic and seismic observations. However, both effects have been shown to impact vertical uplift 45,51 and therefore may influence tsunami behaviour. To address this, we model tsunami genesis and propagation sourced by dynamic earthquake rupture as detailed in 'Tsunami modelling' in Methods. Synthetic tsunamis are generated from the time-dependent horizontal and vertical coseismic seafloor displacements.
Relative to the dynamic-rupture scenario described previously, we vary the yield strength of near-trench sediments in two alternative dynamic scenarios with weaker and stronger near-trench sediments, yielding, respectively, more or less off-fault deformation (Fig. 4a). Cohesion and bulk friction depend on local This study, dynamic rupture scenario Rhie et al. 19 Chlieh et al. 18  mineralogy and lithology, depth, fluid pressure, stress orientation and magnitude, and the respective damage level of the host rock within and beyond the subduction channel. In all models, the bulk cohesion is a combination of a depth-dependent term, C 1 (z)σ ′ c with z the vertical coordinate, positive upwards, σ ′ c the effective confining stress, which accounts for rock strengthening with depth, and a constant term C 0 , which controls localization of off-fault yielding at shallow depth, as detailed in Off-fault plasticity in Methods. Mineralogical and lithological specific laboratory measurements of cohesion and its dependence on confining stress are challenging and may reveal strong local heterogeneity 52 . For simplicity, in the base dynamic-rupture-tsunami scenario, we assign C 0 = 1 MPa and C 1 (z) = 1, appropriate for partially consolidated sediments 35 . Slip to the trench (Fig. 4b) is limited by asymmetric, widespread plastic yielding in the shallow part of the accretionary wedge (Fig.  4a). The alternative scenarios result in distinctly different distributions of shallow slip in the upper 15 km (Fig. 4b). Stronger sediments (C 0 = 10 MPa and C 1 (z) = 1 representing lithified sediments entering the trench) cause localized and weak, yet effective (see Extended Data Fig. 5 for a comparison with a fully elastic model) off-fault yielding, and fault slip locally exceeds 30 m at the trench. By contrast, weaker sediments (C 0 = 0.3 MPa, C 1 (z) = max(1, (z/18,000) 2 ) and at shallow depth decreased bulk internal friction, representing unconsolidated sediments) prevent near-trench slip but cause widespread plastic deformation off the megathrust and off the splay faults (Extended Data Fig. 6d). In addition, the shallow rupture speed in the uppermost 10 km of this scenario (Extended Data Fig. 7a,b) is reduced, averaging only 0.5 km s -1 , which is close to the characteristic slow speed of tsunami earthquakes 53 .
While all three scenarios vary near the trench, they yield indistinguishable teleseismic (Extended Data Fig. 3) and comparable geodetic synthetics (Extended Data Fig. 6). They also are equivalent in terms of moment magnitude (Extended Data Table 1) and macroscopic earthquake kinematics (Extended Data Fig. 7c); however, rupture duration is 100 s shorter for the strong sediment scenario compared with the base scenario. The radiation efficiency of the flat, shallow parts of the megathrust is low in all scenarios.
We evaluate the dynamic-rupture-tsunami scenarios against geodetic tsunami observations. All three dynamic-rupture scenarios source a large-scale tsunami that propagates across the entire simulated domain (Fig. 4c). The amplitudes and shorter period characteristics of the dynamically sourced tsunami synthetics are sensitive to both fault slip at the trench and off-fault plastic deformation. The tsunami sourced by the base scenario and the stronger-sediment scenario match timing and amplitude of water height anomalies captured by the Jason-1 satellite about 2 hr after the mainshock (Fig. 4d). The tsunami sourced by the weaker-sediment earthquake scenario featuring pronounced off-fault yielding and no slip to the trench produces tsunami waves amplified by a factor of 2.
All three scenarios capture the pronounced two-peak characteristic of the first tsunami wave, a feature that is recovered by few of the kinematic models in ref. 27 . However, the overall better agreement of the base and stronger-sediment scenarios reconciles previous stochastic analysis 38 , implying that surface faulting in the southern rupture zone best reproduces the double peak in the altimetry data and is consistent with the inference of lithified sediment in this region 37,54 .
We identify key factors controlling megathrust earthquake dynamics and tsunamigenic potential using multi-physics numerical modelling capable of incorporating sophisticated imaging of complex lithological and fault structure in combination with a stringent framework to determine earthquake initial conditions. The such-found conditions governing a geometrically complex model of the 2004 Sumatra-Andaman earthquake and Indian Ocean tsunami include elevated fluid pressure and reveal a high friction drop on a statically strong and well-oriented megathrust interface. A simple linear-slip weakening friction law combined with off-fault yielding captures key features of the seismic, geodetic and tsunami observations in the near and far fields.
Natural megathrusts vary in additional ways due to local complexities in slab morphology and convergence rates, as well as sediment cover. Our approach yields smoothly varying rupture dynamics and kinematic characteristics, such as fault slip and rupture speed, across the megathrust. While spatial variations in geometric and rheological complexities are found to correlate on length scales of hundreds of metres in 3D imaging 11 , smaller-scale geometrical complexity of the fault interfaces may produce more heterogeneity in these features, as inferred in some data-driven kinematic models. The subduction interface remains nearly optimally oriented in this model (Extended Data Fig. 9), and the variation of absolute prestress magnitudes is small. The variations in regional-scale megathrust geometry are counteracted by the observationally constrained variations in the regional driving stresses. This approach can be extended to include further sources of spatio-temporal heterogeneities, which may modulate prestress and thus megathrust rupture dynamics. This can include variations in seismic coupling 55,56 and subducted structural features 57 , as well as stress changes due to aseismic creep and previous earthquakes 58, 59 .
Both sediment load entering the subduction channel and lithification of those sediments have been connected to the near-surface behaviour and tsunamigenic potential of subduction zones. While high sedimentation has been linked statistically to large earthquakes 60-63 , the scenarios we present here emphasize that the strength of these sediments is key. Indeed, we find that depth-dependent rigidity and the amount of off-fault plastic yielding of near-trench rocks are dominant factors controlling earthquake dynamics and tsunami genesis. Specifically, the Bengal fan to the north is associated with increased sediment thickness 64,65 , and a thick layer of lithified sediment has been evidenced in the southern part of the Andaman trench 54 . However, the nature and amount of sediments certainly vary along the megathrust. We expect that taking along-arc variations of rigidity and sediment yield strength into account will lead to along-arc variations in shallow rupture dynamics and tsunami genesis 45 .
The physics-based numerical models presented here reveal relationships between megathrust fault geometries, tectonic loading, rigidity variations with depth and sediment strength with megathrust hazard and emphasize the importance of regional-scale structural heterogeneity for earthquake dynamics and tsunami genesis. Novel seafloor observatories 66 have the potential to provide detailed images of near-trench regions worldwide. Our computational capabilities render possible the incorporation of present and emerging high-resolution observations into dynamic-rupture-tsunami models. We have demonstrated for the Great Sumatra-Andaman earthquake and tsunami that dominant factors controlling megathrust hazards can be identified, which, in principle, can be applied to subduction zones globally. Incorporating smaller-scale heterogeneities and observational uncertainties into our models will enhance both scenario-focused modelling and physics-based probabilistic hazard assessment of the future.

online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41561-021-00863-5.

Methods
Regional lithological structure. We construct a 3D subsurface computational model of the Sumatra-Andaman region to capture spatial variations of lithology and seismic rock properties. The model incorporates a complex subduction interface, large-scale layers varying in thickness 70 , and topography and bathymetry data 71 ( Supplementary Information, section 3D velocity model). We add a 3-km-wide layer to the velocity structure enclosing the subduction interface. Within this approximated subduction channel, strong along-depth rigidity variations are incorporated (Fig. 1e), which we constrain from recent global V p velocity models 5 in agreement with earlier global inferences from earthquake rupture duration 1, 9 We increase rigidity piecewise linearly with depth 5 from 5 to 43 MPa, thus accounting for softer rocks at shallow depth, including materials scraped off the oceanic crust and younger deposits (Fig. 1e). Ref. 5 empirically derives density and shear wave velocity from V p , which includes the assumption that Poisson ratio increases from 0.25 to 0.35 towards shallower depths. Dynamic-rupture modelling increasing Poisson ratio has been associated with enhanced shallow fault slip and tsunami potential 72 . In addition, the global average in ref. 5 is defined with respect to an absolute interplate boundary depth. In our subduction channel velocity model, we approximate the seafloor depth with respect to the water surface as ranging from 5 km in the south to 3 km in the northern model region, accounting for along-arc variation 71 .
We demonstrate the importance of the assumed 3D velocity structure including the strong rigidity variations in the subduction channel in 'The importance of strongly depth-dependent rigidity' by comparing with a 1D preliminary reference Earth model (PREM 73 ) velocity structure scenario.

Fault geometry.
We construct a non-planar model of the subduction interface and possible splay faults that may have ruptured during the 2004 earthquake (Figs.1c and 2a). The slab interface is based on the geometry of slab 2.0 74 . Evidence that splay faults slipped coseismically during the 2004 event stems from large tsunami run-ups observed in Aceh province in the near-source region 75 , high reflectivity in shallow seismic reflection data 31 , and aftershocks 33,76 , the double peak in the tsunami waves recorded by the Jason-1 satellite 32 , alongside other observations 29 . We add one long forethrust dipping 65° landwards and two backthrust faults dipping 65° trenchwards. Each extends from its inferred trace to the megathrust interface. The forethrust unifies the upper splay fault mapped by ref. 29 and the splay fault suggested by ref. 32 . The northern backthrust is mapped from seismic reflection data by refs. 29,30 , and the southern backthrust is identified by ref. 31 from seismic reflection data. The true extents of these faults are not known; mapped lengths are restricted to the area of data coverage.

Stress state and fault friction.
The stresses acting on faults and their strength, which are key initial conditions of dynamic-rupture models, are poorly known. Laboratory experiments offer insight on frictional strength parameters, but extrapolating these results to a natural scale is complicated. We introduce new procedures to constrain these parameters from stress inversion results 68 , including the inference of coseismic rotation of the principal stress axes 4,68 .
Principal stress orientations and magnitudes. In line with the results of ref. 68 , we apply a non-Andersonian initial stress field that features a shallowly plunging maximum compressive stress axis (σ 1 ), a near-horizontal intermediate stress axis (σ 2 ) and a steeply plunging minimum principal stress (σ 3 ). σ 1 plunges at 22°, which optimally stresses a fault dipping 8°.
We align stress inversion and geodetic observations in choosing the σ 1 orientation. The model accounts for the inferred clockwise rotation of σ 1 from 309° in the south to 330° in the north 68 , with a linear transition over 2° near 5° latitude (Fig. 1a). This leads to synthetic displacements with orientations in overall agreement with geodetic observations. The principal stress magnitudes are systematically constrained on the basis of seismo-tectonic observations, fault fluid pressurization and the Mohr-Coulomb theory of frictional failure.
Rupture dynamics are in their general sense governed by the stress drop and the fault strength. These are constrained by the effective confining stress, the fault friction drop μ s − μ d (where μ s and μ d are the static and dynamic friction coefficients, respectively) and the relative prestress ratio R 0 . The relative prestress ratio is defined as the ratio of prestress to strength drop 77 : where τ 0 and σ n are the initial shear and normal tractions on a fault plane; R 0 is the maximum value of R, associated with a virtual, optimally oriented fault plane. We assume lithostatic pressure and a near-lithostatic fluid pressure, P f = 0.974σ c , to ensure a realistic stress drop. The resulting effective confining stress increases with depth with a gradient of ~1 MPa km -1 . From 25 to 45 km depth, we taper the deviatoric stresses to zero to represent the transition from a brittle to a ductile regime.
We account for along-arc variations in convergence rate by varying the initial stress along strike of the megathrust. Observations by ref. 68 suggest that the strain accumulated since the last great earthquake, inferred using along-arc variations in convergence rate, is consistent with the along-arc variation in slip inferred for the 2004 earthquake, which is characterized by more slip in the south than in the north.
We modulate R 0 by the convergence rate inferred from rigid plate tectonics considerations, using the Euler pole inferred by ref. 69 (Fig. 1b). This modulation is applied over the full rupture length.
We demonstrate the importance of along-arc prestress modulation in 'The importance of convergence-rate-modulated regional driving stresses' by comparing with a scenario in which no prestress modulations are considered.
Coseismic stress rotation to constrain fault friction. The inference of coseismic rotation of the principal stress axes (Fig. 1c) during the 2004 earthquake 4,68 constrains the ratio of the stress drop over the prestress. A substantial principal stress rotation indicates that the stress drop is large enough to change the deviatoric background stress 78 ; that is, the prestress and the stress drop are of the same order. This can lead to a sudden switch of orientation of principal axes of stress, as suggested by the occurrence of normal faulting following megathrust earthquakes 4,68 .
We estimate the average stress drop dτ in the dynamic-rupture model as We define ξ as the ratio of residual shear traction to prestress to get which can be written as Equation (5) relates μ d to μ s and R 0 . Assuming a static friction value of μ s = 0.6 (ref. 79 ), a ratio of residual shear traction to prestress of ξ = 0.4, as inferred in the southern part of the rupture 68 , and a ratio of prestress to strength drop of R 0 = 0.65, yields μ d = 0.18. On the basis of these considerations, we use a linear-slip weakening friction law setting μ s = 0.6, μ d = 0.2 and D c = 2.5 m. The relatively large D c allows recovering the northern rupture arrest, as well as an average rupture speed comparable to published inferences (Fig. 2d). It leads to a wide process zone in our model, which may reflect the larger slip pulse width associated with large megathrust earthquakes 80 . At shallow depth above 10 km, we linearly increase μ d . Fault friction is then slip neutral at 5 km depth and slip strengthening above. In the southern rupture region, transition to slip strengthening occurs 5 km deeper.
Earthquake computational model. Computational earthquake models that are able to account for the curved thrust interface and splay fault networks intersecting with bathymetry, the narrow accretionary wedge and the complex lithology of subduction zones are challenging. Off-fault yielding processes, which can substantially modulate fault slip near the trench 45,81 , pose additional computational demands. These challenges are addressed using the SeisSol software package [82][83][84] , which solves the nonlinear coupled problem of spontaneous frictional failure across prescribed fault surfaces, seismic wave propagation and (visco-)plastic Drucker-Prager off-fault damage 42 .
SeisSol is based on an arbitrary high-order accurate derivative discontinuous Galerkin method 82 and solves the seismic wave equations with high-order accuracy in space and time. SeisSol uses fully non-uniform, adaptive, unstructured tetrahedral meshes enabling geometrically complex models, including curved and mutually intersecting faults. Mesh resolution can be adapted to ensure fine sampling of the faults while satisfying the requirements regarding numerical dispersion of pure wave propagation away from the fault. End-to-end computational optimization [85][86][87][88] , including an efficient local time-stepping algorithm 84 , allows for high efficiency on high-performance computing infrastructure. SeisSol is verified in a wide range of community benchmarks 89 by the Southern California Earthquake Center/USGS Dynamic Rupture Code Verification project 90,91 . SeisSol is freely available (https://github.com/SeisSol/ SeisSol).
Tsunami modelling. To model the tsunami, we use GeoClaw 92,93 , a highly scalable software solving the shallow-water equations using adaptive mesh refinement. GeoClaw implements a shock-capturing finite volume method on logically rectangular latitude-longitude grids on the sphere. Its adaptive mesh refinement approach allows for stable and accurate simulation of large-scale wave propagation in the deep sea and small-scale wave shoaling and inundation at the shore.
We use the GEBCO 2019 topography and bathymetry dataset 71 with a horizontal resolution of 15 arcsec (approximately 450 m). This allows for a sufficiently accurate representation of bathymetric features in deeper sea regions, where tsunamis have large wavelengths, but is relatively inaccurate at shallow depths. The tsunami is sourced by time-dependent bathymetry perturbations resembling the coseismic dynamic-rupture seafloor displacements 94 . In addition to the vertical displacement, we account for the effects of east-west and north-south horizontal displacements by the linear approximation proposed by ref. 95 .
The domain of the computational tsunami model (Fig. 4c) encompasses the source region and a large part of the Indian Ocean. The computational grid incorporates two levels of dynamically adaptive mesh refinement, with respective resolutions of 45 and 90 arcsec (~1,300 m and 2,600 m). The simulation is run for 2 h 30 min (simulation time).
Our hydrostatic shallow-water-based tsunami modelling approach may be extended to capture also smaller-scale complexity during tsunami genesis 96 and dispersive effects during tsunami wave propagation 97 by allowing for a nonlinear hydrodynamic response 98,99 , corrections for dispersive Earth elasticity and non-dispersive water compressibility 100 or fully coupled seismic, acoustic and gravity modelling 101,102 .
Teleseismic model validation. Following ref. 43 , we translate the time histories of dynamic fault slip of the dynamic-rupture scenario into a subset of 50 double couple point sources (25 along strike by 2 along depth). To that purpose, we subdivide the fault in 50 subregions. We compute the moment rate as G dS dt where S is the slip and G the local rigidity in the dynamic-rupture model at the centre of each fault element face. The equivalent moment tensor of each fault element face is calculated from the rake of the accumulated slip and each fault face orientation in terms of strike and dip angle. The moment tensor of each subregion is computed by summing the moment tensors of each fault element face within a subregion. The respective location of the equivalent moment tensor is calculated as the mean of the locations of the barycentres of each slipping fault element face weighted by its moment magnitude. We generate broadband synthetics from a Green's function database using Instaseis 103 , which is based on global 2.5D axisymmetric spectral element simulations, and the PREM model 73 for a maximum period of 10 s and including anisotropic effects.
The in-this-manner rapidly computed, yet simplified, teleseismic synthetics reproduce the recorded long-period teleseismic waveforms with an average rRMS of 0.7 for the preferred dynamic-rupture scenario at periods of 150 s to 500 s (Extended Data Fig. 3). The arrival times and amplitudes of the surface waves at the frequency considered (150 s to 500 s) are well recovered at most stations. The envelopes of the recorded long-period teleseismic waveforms are also overall well captured, with noticeable discrepancies in rupture forward direction (specifically at seismic stations HYB, LSA and ENH).
The teleseismic waveform synthetics combine InstaSeis Green's functions based on a PREM model with dynamic-rupture equivalent moment-tensor point sources and thus do not account for the off-fault 3D velocity structure we use in the dynamic model. Accounting accurately for the seismic potency contribution of the low-rigidity shallow parts of the megathrust would require fully coupled dynamic-rupture and seismic wave propagation simulations across teleseismic distances. However, this is currently computationally not feasible in the scope of thorough sensitivity analysis.
Off-fault plasticity. We account for the possibility of off-fault inelastic energy dissipation by adopting a Drucker-Prager elasto-viscoplastic rheology. The implementation is described in ref. 88 . While purely elastic dynamic-rupture models require defining initial stresses acting on the fault(s) only, models with off-fault plasticity require domain-wide parameterization of initial stresses, cohesion and friction of the bulk rock.
Off-fault initial stresses resemble the depth-dependent initial stresses described in Stress state and fault friction. The viscoplastic relaxation time T v is set to 0.03 s.
Cohesion and bulk friction are challenging to constrain. For example, the strength of drill core samples, which are often free of preexisting cracks, may not be representative of the fractured rock mass in a fault zone 104 . High-resolution geophysical imaging is required to resolve the geometrical complexities and characterize the damage state of megathrust fault zones 11 . In line with our focus on regional driving factors, we here use a simple cohesion model averaging local effects with a limited number of parameters. We assume that effective off-fault deformation shall be limited to shallow depths (above15 to 20 km depth) by adapting a linear dependence of bulk cohesion on the effective confining pressure 105 . We introduce an additional constant term to allow for modulation of the severity of plastic deformation at shallow depths in between scenarios. Then, the bulk cohesion C(z) is set as (Extended Data Fig. 8a): C 1 (z) accounts for the hardening of the rock structure with depth, while C 0 controls localization and magnitude of off-fault yielding at shallow depth and σ ′ c = σc − P f is the effective confining stress. We set C 0 to 1.0 MPa in the base dynamic-rupture scenario to represent partially consolidated sediments 35 , to 0.3 MPa for the alternative scenario featuring weaker sediments, and to 10.0 MPa for the alternative model featuring stronger sediments. C 1 (z) = 1 for both the base model and the alternative model featuring stronger sediments and C 1 (z) = max(1, (z/18,000) 2 ) for the alternative model featuring weaker sediments.
To enhance the shallow effectiveness of off-fault plasticity without causing widespread yielding at depths greater than 15-20 km, we decrease the internal frictional shear strength of the bulk rock in the weak sediment model. The bulk material's internal friction coefficient ν is set equal to the static fault friction coefficient, μ s = 0.6. In the weaker sediments scenario only, we linearly decrease it at shallow depth, from 0.6 at 18 km to 0.2 at the sea surface. Exploring additional hypotheses or observationally driven cohesion models is left to future work. For example, ref. 105 proposes a procedure based on the Geological Strength Index and the Hoek-Brown model that may be adapted to constrain a locally heterogeneous cohesion model.
Closeness to failure (Extended Data Fig. 8a) is quantified by the CF ratio 81 : where I 2 is the second invariant of the deviatoric stresses, and τ c is the Drucker-Prager yield criterion, given by: with Φ = arctan(ν) the internal angle of friction and σm = ∑ 3 n=1 σii/3 the mean stress.
The such-constrained CF ratio is depth dependent and varies along arc for all three scenarios presented, aligned with the assumed distribution of initial stresses. In the base and stronger-sediment scenarios, cohesion is dominated by the constant term C 0 at shallow depth. Closeness to failure is there low, while plastic yielding is not prevented since absolute rock strength is low as well. In the weak sediment model, a quadratic reduction of cohesion C(z) down to a depth of 18 km leads to pronounced higher closeness to failure peaking at 4 km depth.
The total seismic moment M 0,t is the sum of the moment due to slip on the fault, M 0,e , and the moment due to off-fault plastic strain, M 0,p . We adapt previous 2D inferences 106 to our 3D non-associated Drucker-Prager plasticity implementation 88 and compute: with μ being the rigidity, V the volume of each tetrahedral element i and η a scalar quantity measuring the accumulated material damage at the end of the dynamic-rupture simulation. Following ref. 107 , it is with ε p ij being the inelastic strain rate. The contribution of plastic strain to the total moment is small for the base and stronger sediments scenarios but non-negligible (M 0,p /M 0,t ≈ 7.5%) in the weaker-sediment scenario (Extended Data Table 1).
Ratios of M 0,p /M 0,e , which are on the order of a few percent, are consistent with 2D pulse-like dynamic-rupture simulations 106 at comparable relative fault strength S, with S = (1 − R)/R. However, the contribution of plastic moment can become dominant under certain conditions, such as large S and large angles between the maximum compressive stress and the fault. This may be explored in future work focusing specifically on relatively slow tsunami earthquakes 53 .
Off-fault plastic deformation affects all here presented scenarios compared with a purely elastic model (Extended Data Fig. 5a). Specifically, earthquake rupture dynamics of a fully elastic scenario, with no off-fault yielding, differ from the stronger-sediment scenario as illustrated in Extended Data Fig. 5a: specifically near the trench, fully elastic fault slip is considerably increased to a maximum of 52 m. Fully elastically sourced tsunami synthetics (Extended Data Fig. 5b) are overall similar to the base and stronger-sediment scenarios and also reproduce the early two tsunami peaks.
The importance of strongly depth-dependent rigidity. We demonstrate the importance of strongly depth-dependent rigidity in an alternative model, which adopts a 1D PREM 73 velocity structure. To generate a comparable dynamic-rupture earthquake scenario requires few changes to the model. We increase the nucleation radius from 10 to 15 km and the maximum normalized prestress ratio R 0 from 0.65 to 0.75 to restore the rupture potential of the megathrust. The PREM alternative earthquake scenario yields a comparable slip distribution to the base model (Extended Data Fig. 2a) and therefore fits GPS observations well (Extended Data Fig. 2d,e). However, its M w = 9.32 is higher than most kinematic models (based on body and surface waves, normal models, GPS, corral uplift or Rayleigh waves 26 ). Rupture speed of the PREM model is faster, which is shortening rupture duration to 450 s. This is 50-100 s shorter than most inferences from normal, body and surface wave inversion as well as back-projection 17,21,25,108 . In contrast to the along-depth variation of rupture speed and the boomerang-shaped slip pulse of the base scenario, rupture speed and slip rates in the PREM scenario (Extended Data Fig. 2b and Supplementary Video 2) are rather homogeneous along depth and along margin. In addition, the PREM slip pulse is smoother and wider. The Fig. 3 | Teleseismic waveforms. (a) Comparison of synthetic (red and blue) and observed (black) teleseismic waveforms. red waveforms are computed using the base scenario, while blue waveforms correspond to the model featuring less off-fault yielding and more slip to the trench, yet, indistinguishable synthetics, shown in Fig. 4. A 150-500 s band-pass filter is applied to all traces. Synthetics are generated using Instaseis 103 and the PrEM model including anisotropic effects and a maximum period of 10 s. (b) Locations at which synthetic data are compared with observed records. Fig. 4 | comparison of synthetic ground displacements in the base scenario and geodetic observations. Synthetics are shown in blue, geodetic observations in orange and magenta. In the far-field a different scaling is applied. (a) Horizontal ground displacements. (b) Vertical ground displacements. The complete uplift and subsidence resulting from the dynamic rupture scenario (in m) is shown in red to blue including noticeable splay faulting signatures. Fig. 8 | off fault plasticity. (a) Depth dependence of bulk cohesion C(z) (Eq. 6, blue) and of the failure criterion CF (Eq. 7, warm colors) at four locations along the trench (from south to north, the curves are colored in gold, orange, brown and black). Different line styles refer to different scenarios. (b) and (c) Off-fault plastic strain (quantified as η, Eq. 10) accumulated in the base and 2 alternative earthquake scenarios with differing visco-elastic characteristics yielding stronger or weaker sediments. We use a linear scale, cut off at η > 10 −5 . The red line is the vertical seafloor uplift (not to scale, amplified × 200). The locations of the southern (b) and northern (c) cross-sections are illustrated beneath. The distribution of η across the central rupture region is shown in the main text Fig. 4a.