3D Linked Subduction, Dynamic Rupture, Tsunami, and Inundation Modeling: Dynamic Effects of Supershear and Tsunami Earthquakes, Hypocenter Location, and Shallow Fault Slip

Physics-based dynamic rupture models capture the variability of earthquake slip in space and time and can account for the structural complexity inherent to subduction zones. Here we link tsunami generation, propagation, and coastal inundation with 3D earthquake dynamic rupture (DR) models initialized using a 2D seismo-thermo-mechanical geodynamic (SC) model simulating both subduction dynamics and seismic cycles. We analyze a total of 15 subduction-initialized 3D dynamic rupture-tsunami scenarios in which the tsunami source arises from the time-dependent co-seismic seafloor displacements with flat bathymetry and inundation on a linearly sloping beach. We first vary the location of the hypocenter to generate 12 distinct unilateral and bilateral propagating earthquake scenarios. Large-scale fault topography leads to localized up- or downdip propagating supershear rupture depending on hypocentral depth. Albeit dynamic earthquakes differ (rupture speed, peak slip-rate, fault slip, bimaterial effects), the effects of hypocentral depth (25–40 km) on tsunami dynamics are negligible. Lateral hypocenter variations lead to small effects such as delayed wave arrival of up to 100 s and differences in tsunami amplitude of up to 0.4 m at the coast. We next analyse inundation on a coastline with complex topo-bathymetry which increases tsunami wave amplitudes up to ≈1.5 m compared to a linearly sloping beach. Motivated by structural heterogeneity in subduction zones, we analyse a scenario with increased Poisson's ratio of ν = 0.3 which results in close to double the amount of shallow fault slip, ≈1.5 m higher vertical seafloor displacement, and a difference of up to ≈1.5 m in coastal tsunami amplitudes. Lastly, we model a dynamic rupture “tsunami earthquake” with low rupture velocity and low peak slip rates but twice as high tsunami potential energy. We triple fracture energy which again doubles the amount of shallow fault slip, but also causes a 2 m higher vertical seafloor uplift and the highest coastal tsunami amplitude (≈7.5 m) and inundation area compared to all other scenarios. Our mechanically consistent analysis for a generic megathrust setting can provide building blocks toward using physics-based dynamic rupture modeling in Probabilistic Tsunami Hazard Analysis.

Physics-based dynamic rupture models capture the variability of earthquake slip in space and time and can account for the structural complexity inherent to subduction zones. Here we link tsunami generation, propagation, and coastal inundation with 3D earthquake dynamic rupture (DR) models initialized using a 2D seismo-thermomechanical geodynamic (SC) model simulating both subduction dynamics and seismic cycles. We analyze a total of 15 subduction-initialized 3D dynamic rupture-tsunami scenarios in which the tsunami source arises from the time-dependent co-seismic seafloor displacements with flat bathymetry and inundation on a linearly sloping beach. We first vary the location of the hypocenter to generate 12 distinct unilateral and bilateral propagating earthquake scenarios. Large-scale fault topography leads to localized upor downdip propagating supershear rupture depending on hypocentral depth. Albeit dynamic earthquakes differ (rupture speed, peak slip-rate, fault slip, bimaterial effects), the effects of hypocentral depth (25-40 km) on tsunami dynamics are negligible. Lateral hypocenter variations lead to small effects such as delayed wave arrival of up to 100 s and differences in tsunami amplitude of up to 0.4 m at the coast. We next analyse inundation on a coastline with complex topo-bathymetry which increases tsunami wave amplitudes up to ≈1.5 m compared to a linearly sloping beach. Motivated by structural heterogeneity in subduction zones, we analyse a scenario with increased Poisson's ratio of ν = 0.3 which results in close to double the amount of shallow fault slip, ≈1.5 m higher vertical seafloor displacement, and a difference of up to ≈1.5 m in coastal tsunami amplitudes. Lastly, we model a dynamic rupture "tsunami earthquake" with low rupture velocity and low peak slip rates but twice as high tsunami potential energy. We triple fracture energy which again doubles the amount of shallow fault slip, but also causes a

INTRODUCTION
Earthquake sources are governed by highly non-linear multiphysics and multi-scale processes leading to large variability in dynamic and kinematic properties such as rupture speed, slip rate, energy radiation, and slip distribution (e.g., Oglesby et al., 2000;Kaneko et al., 2008;Gabriel et al., 2012;Bao et al., 2019;Ulrich et al., 2019a;Gabriel et al., 2020). Such variability may impact the generation, propagation, and inundation of earthquake-generated tsunami or secondary tsunami generation mechanisms such as triggered landslides (e.g., Sepúlveda et al., 2020). For example, unexpectedly large slip at shallow depths may generate large tsunami (Lay et al., 2011;Romano et al., 2014;Lorito et al., 2016).
To model earthquake-generated tsunami, sources can be approximated from earthquake generated uplift (Behrens and Dias, 2015, and references therein). Analytical solutions (e.g., Okada, 1985) describe seafloor displacements sourced by uniform rectangular dislocations within a homogeneous elastic half space. Models of tsunami generated by large earthquakes can routinely and quickly use kinematic finite fault models constrained by inversion of seismic, geodetic, and other geophysical data (Geist and Yoshioka, 1996;Ji et al., 2002;Babeyko et al., 2010;Maeda et al., 2013;Allgeyer and Cummins, 2014;Mai and Thingbaijam, 2014;Bletery et al., 2016;Jamelot et al., 2019), but are challenged by the inherent non-uniqueness of kinematic source models (Mai et al., 2016).
3D Dynamic earthquake rupture modeling can provide mechanically viable tsunami source descriptions on complex faults or fault systems on the scale of megathrust events (Galvez et al., 2014;Murphy et al., 2016;Uphoff et al., 2017;Murphy et al., 2018;Ma and Nie, 2019;Saito et al., 2019;Ulrich et al., 2020). Such simulations can exploit modern numerical methods and high-performance computing (HPC) to shed light on the dynamics and severity of earthquake behavior and potentially complement PTHA. For example, in dynamic rupture models shallow slip amplification can spontaneously emerge due to updip rupture facilitated by along-depth bi-material effects (Rubin and Ampuero, 2007;Ma and Beroza, 2008;Scala et al., 2017) and free-surface reflected waves within the accretionary wedge (Nielsen, 1998;Lotto et al., 2017a;van Zelst et al., 2019). Dynamic rupture earthquake models can yield stochastic slip distributions, too, under the assumption of stochastic loading stresses (Geist and Oglesby, 2014). Such physics-based models can be directly linked to tsunami models by using the time-independent or timedependent seafloor displacements (and potentially velocities) as the tsunami source (Kozdon and Dunham, 2013;Ryan et al., 2015;Lotto et al., 2017b;Saito et al., 2019;Madden et al., 2020). For instance, time-dependent 3D displacements from observational constrained dynamic rupture scenarios of the 2018 Palu, Sulawesi earthquake and the 2004 Sumatra-Andaman earthquake are linked to a hydrostatic shallow water tsunami model by Ulrich et al. (2019b). Bathymetry induced amplification of horizontal displacements are thereby accounted for by following Tanioka and Satake (1996).
Observational and numerical studies show that megathrust geometry and hypocenter location influence earthquake rupture characteristics. Ye et al. (2016) state that megathrust earthquakes across faults that are longer horizontally than they are deep vertically (with an aspect ratio of three or larger) tend to exhibit primarily unilateral behavior. Also, events with an asymmetric hypocenter location on the fault favor rupture propagation along strike to its far end (Harris et al., 1991;McGuire et al., 2002;Hirano, 2019). Weng and Ampuero (2019) show the energetics of elongated ruptures is radically different from that of conventional circular crack models. Bilek and Lay (2018) find that the complexity of slip as well as bi-or unilateral rupture preferences of large earthquakes highly depend on the depth location of the hypocenter.
Subduction zones worldwide are associated with tectonic, frictional, and structural heterogeneity along depth and alongarc impacting megathrust earthquake and tsunami dynamics (Kirkpatrick et al., 2020, e.g.,). Kanamori and Brodsky (2004) show that fracture energy varies between subduction zone earthquakes. A special case are so called tsunami earthquakes (Kanamori, 1972) that may require a large amount of fracture energy, low rupture velocity and low radiation efficiency. Structural heterogeneity in subduction zones includes variations in Poisson's ratio (V p /V s ) (e.g., Liu and Zhao, 2014;Niu et al., 2020)), while dynamic rupture and seismic wave propagation models often adopt an idealized Poisson's ratio of ν =0.25 governing seismic wave propagation (e.g., Kozdon and Dunham, 2013).
The initial conditions of dynamic rupture simulations that control earthquake rupture nucleation, propagation, and arrest include fault loading stresses, frictional strength, fault geometry, and subsurface material properties (e.g., Kame et al., 2003;Gabriel et al., 2013;Galis et al., 2015;Bai and Ampuero, 2017). These initial conditions may be observationally and empirically informed (e.g., Aochi and Fukuyama, 2002;Aagaard et al., 2004;Murphy et al., 2018;Ulrich et al., 2020) but remain difficult to constrain. Particularly in subduction zones where observational data are sparse, space and time scales vary over many orders of magnitude and both geometric and rheological megathrust complexities are likely to control rupture characteristics. Recently, initial conditions for 2D and 3D megathrust dynamic rupture earthquake models have been informed from 2D geodynamic long-term subduction and seismic cycle models (van Zelst et al., 2019;Madden et al., 2020;van Zelst et al., 2021). This approach provides self-consistent initial fault loading stresses and frictional strength, fault geometry and material properties on and surrounding the megathrust, as well as consistency with crustal, lithospheric, and mantle deformation and deformation in the subduction channel over geological time scales. Such subduction-initialized heterogeneous dynamic rupture models lead to complex earthquakes with multiple rupture styles (Gabriel et al., 2012), shallow slip accumulation and fault reactivation.
We apply the 2D geodynamical subduction and seismic cycle (SC) model from van Zelst et al. (2019) to inform realistic 3D dynamic rupture (DR) megathrust earthquake models within a complex, self-consistent subduction setup along with their consequent tsunami, following the subduction to tsunami run-up linking approach described in Madden et al. (2020). In this study we introduce a number of important differences to previous work. 2D linking including approximations to match SC and DR fracture energy during slip events leads to differences in slip magnitude between the SC and DR modeling and large magnitudes and high rupture speed in dynamic rupture scenarios (van Zelst et al., 2019). In contrast, we here constrain fracture energy independently from the long-term model. In Madden et al. (2020), a different long-term geodynamic and seismic cycle simulation was used, specifically, assuming different shear moduli. We here change the geometrical 3D extrapolation of the 2D fault geometry compared to the large blind dynamic earthquake scenario of M W 9.0 of Madden et al. (2020), to be consistent with empirical earthquake source scaling relations for M W 8.5 megathrust events (Strasser et al., 2010).
We use complex 3D dynamic rupture modeling to first study trade-offs and effects of along-strike unilateral vs. bilateral rupture and variations in hypocentral depth in subduction zone earthquakes (McGuire et al., 2002). By varying the hypocentral location along arc and along depth, we generate 12 distinct unilateral and bilateral earthquakes with depth-variable slip distribution, rupture direction, bimaterial, and geometrical effects in the dynamic slip evolution. We analyse the consequent time-dependent variations in seafloor uplift affecting tsunami propagation and inundation patterns. We define as reference model a bilateral, deeply nucleating earthquake. To this reference model we add a complex and more realistic coastline in the tsunami simulation and study the effects on tsunami arrival time and wave height at the coast.
The linkage from long-term geodynamic to co-seismic dynamic rupture modeling requires assumptions with respect to the incompressibility and visco-elasto-plastic, plane-strain conditions of the subduction model vs. the compressible, elastic conditions of the earthquake model. In two additional scenarios we analyse variations in the energy balance of the subductioninitialized dynamic rupture scenarios. We increase fracture energy in the reference model by changing the frictional critical slip distance within the dynamic rupture model and adapting nucleating energy accordingly. The increase in fracture energy leads to large uplift, low radiation efficiency and low rupture velocities, characterizing a tsunami earthquake (Kanamori, 1972). Lastly, we analyze the effect of a higher Poisson's ratio throughout the dynamic rupture reference model and the effect on tsunami genesis and inundation.
This leads to a total of 15 subduction-initialized 3D dynamic rupture-tsunami scenarios: 12 dynamic rupture models with varying hypocenters. For one "reference model" (model 3B) of these 12 we vary fracture energy or Poisson's ratio, or coastline bathymetry.

METHODS
Here, we summarize the computational methods used for simulating subduction-initialized dynamic earthquake rupture linked to tsunami generation, propagation, and inundation (Figure 1). For an in-depth description of the virtual laboratory for modeling tsunami sources arising from 3D co-seismic seafloor displacements generated by dynamic earthquake rupture models, we refer to Madden et al. (2020). We compute 3D dynamic earthquake rupture and seismic wave propagation with SeisSol (https:/seissol.org). Tsunami propagation and inundation uses sam(oa) 2 -flash, which is part of the open-source software sam(oa) 2 (https://gitlab.lrz.de/samoa/samoa). Both codes use highly optimized and parallel implementations of discontinuous Galerkin (DG) schemes. All simulations were performed on SuperMUC-NG at the Leibniz Supercomputing Centre Garching, Germany.
To link input and output data in massively parallel simulations, we use ASAGI (pArallel Server for Adaptive GeoInformation), an open source library with a simple interface to access Cartesian material and geographic datasets (Rettenberger et al., 2016, www.github.com/TUM-I5/ASAGI). ASAGI translates a snapshot of the 2D subduction model into 3D initial conditions for the earthquake model and bathymetry data and seafloor displacements from the earthquake model into initial conditions for the tsunami model. ASAGI organizes Cartesian data sets for dynamically adaptive simulations by automatically migrating the corresponding data tiles across compute nodes as required for efficient access.

3D Earthquake Dynamic Rupture Modeling With SeisSol
Physics-based 3D earthquake modeling captures how faults yield, slide and interact (e.g., Ulrich et al., 2019a;Palgunadi et al., 2020) and can provide mechanically viable tsunami-source descriptions (Ryan et al., 2015;Ma and Nie, 2019;Ulrich et al., 2019bUlrich et al., , 2020. We use SeisSol (de la Puente et al., 2009;Pelties et al., 2012Pelties et al., , 2014 to solve simultaneously for frictional failure across prescribed fault surfaces and high-order accurate seismic wave propagation in space and time (illustrated in Figure 1, right). SeisSol uses a discontinuous Galerkin (DG) scheme with Arbitrary high-order DERivative (ADER) time stepping on unstructured tetrahedral grids with static mesh adaptivity Käser and Dumbser, 2006). It is thereby particularly suited for modeling complex geometries such as those in the vicinity of subducting slabs. SeisSol is optimized for current petascale supercomputers Heinecke et al., 2014;Rettenberger et al., 2016;Uphoff and Bader, 2020;Dorozhinskii and Bader, 2021) and uses an efficient local time-stepping algorithm (Breuer et al., 2016;Uphoff et al., 2017;Wolf et al., 2020). Its accuracy is verified against a wide range of community benchmarks (Harris et al., 2011(Harris et al., , 2018, including dipping and branching faults with heterogeneous off-fault material and initial on-fault stresses Wollherr et al., 2018;Gabriel et al., 2020). We note that on-fault initial conditions such as frictional parameters or initial fault stresses are assigned with sub-elemental resolution (at each DG Gaussian integration point). However, within each off-fault element, all material properties are constant. We create a 3D complex structural model in GoCad (Holding, 2018) and discretize it with the meshing software Simmodeler by Simmetrix (Simmetrix Inc., 2017).
Within SeisSol, frictional failure is treated as an internal boundary condition for which the numerical solution of the elastodynamic wave equation is modified. In the dynamic rupture scenarios of our study, fault strength, i.e., its yielding and subsequent frictional weakening, is governed by the widely adopted linear slip weakening (LSW) friction law (Ida, 1972). Over a critical slip weakening distance D c the effective friction coefficient µ decreases linearly from static µ s until reaching dynamic µ d . We note that this is different to the rate-weakening friction used in the long-term geodynamic SC model. The process zone width is the inherent length scale defining the minimum resolution required on-fault, and is defined as the area behind the rupture tip in which shear stress decreases from its static to its dynamic value . Figure 2 depicts the inferred 3D initial conditions from the subduction seismic cycle model for all dynamic rupture scenarios. These include highly heterogeneous initial shear stress and strength as well as fault geometry and material structure that together govern earthquake nucleation, propagation, and arrest. The underlying 2D seismo-thermo-mechanical geodynamic seismic cycle (SC) model simulates subduction dynamics over millions of years and earthquake cycles over several hundreds of years (e.g., Van Dinther et al., 2013 (Figure 2A). The long-term phase of the simulation builds up stress as well as selfconsistent strength and fault geometries across the subduction interface and forearc. The geodynamic SC model includes 70 megathrust events that rupture almost the entire fault and nucleate near the downdip seismogenic zone limit. The same representative event as in van Zelst et al. (2019) is here chosen and linked to the newly designed dynamic rupture model applying the techniques described in Madden et al. (2020). We note, that at the chosen SC time step, the dominant deformation mechanism in the seismogenic zone is elastic behavior, which is consistent with the deformation mechanism in the dynamic rupture model. All material properties, stresses, and fault geometries are exported from the SC model at that timestep. Material density and fault strength can be adapted to be used in the DR model. Cohesion varies between 2.5 and 20 MPa and increases with deeper lithologies. At nucleation depth, cohesion is set to 5 MPa. The curved, blind megathrust interface evolves during this slip event and is characterized by large-scale fault roughness including characteristic topographic highs ("bumps") in terms of a distinct change in local gradients of the curved non-planar interface related to sediment intrusion on geodynamic time scales (zoomin box of Figure 2A). In contrast to van Zelst et al. (2019) and Madden et al. (2020), the critical slip weakening distance D c is here not inferred from the geodynamic SC model directly but is assigned to be a constant value (D c = 0.1 m for 13 out of 14 dynamic rupture models) along the entire fault. Increasing the critical slip weakening distance D c , over which effective friction is reduced to its dynamic value, results in a higher fracture energy G c according to G c = 1/2µ s P n − µ d P n D c with P n being the initial fault normal stress (Venkataraman and Kanamori, 2004). The fracture energy consumed within the frictional process zone indicates how much energy is necessary to initiate and sustain rupture propagation. A high fracture energy results in higher fault slip, higher stress drop and a higher moment magnitude for comparable dynamic rupture scenarios (given the nucleation energy is adopted accordingly). In one model, model 5, we triple the model-wide constant D c from 0.1 to 0.3 m while keeping all other parameters constant, thus, tripling fracture energy.

Subduction Seismic Cycle Modeling for Earthquake Initial Conditions
All material properties are extrapolated into the third dimension as constant along arc, for simplicity. We use a plane strain assumption, as in the 2D subduction model, to determine the out-of-plane normal and shear stress components in 3D. This implies that we omit eventual oblique subduction components by setting the out-of-plane shear stresses to zero and the out-ofplane normal stress component to be a function of the two inplane normal stresses and Poisson's ratio ν. In the SC subduction model a Poisson's ratio of ν = 0.5 is used, which is an appropriate assumption for large time scales. This Poisson's ratio needs to be reassigned within the DR model to represent compressible rocks and to solve the linear elastic wave equations. The choice of ν affects the material properties as they are transferred from the subduction model to the earthquake model since Lame's parameter is calculated from the re-assigned ν and the shear modulus G is taken from the subduction model. In all besides one dynamic rupture models we assume ν = 0.25 (Poisson solid, λ = G). In one 3D dynamic rupture (model 6, modified from reference model 3B), we use a larger Poisson's ratio of ν = 0.3, which is in the range observations of basaltic rocks (Gercek, 2007).
Using the imported on-fault conditions from the geodynamic SC model directly would lead to multiple locations of instantaneous failure on the fault ( Figure 2B). We thus use a pre-processing static relaxation step in which we relax the initial fault loading stresses to be just below fault strength before applying them to the DR model following Wollherr (2018). Spontaneous earthquake rupture is commonly initialized by assigning an over-stressed or frictionally weaker predefined hypocentral "nucleation" area. We here assign a locally lower static friction coefficient within a circular patch of an empirically determined minimum required size to initiate spontaneous rupture. We choose frictional coefficients µ s and µ d within this patch constrained from the geodynamic SC model. We prescribe all but the shallowest (model family 1) hypocentral nucleation depths at locations which are close to failure (low strength excess, Figure 2B) in the SC model. The reference geodynamic SC slip event nucleates at 40 km depth, corresponding to model family 3 in this study.

Tsunami and Inundation Modeling With Samoa
Tsunami are modeled with a limited second order Runge-Kutta discontinuous Galerkin solver (Cockburn and Shu, 1998;Giraldo and Warburton, 2008) for the two dimensional depthintegrated hydrostatic non-linear Shallow Water Equations (SWE, LeVeque et al., 2011) in the sam(oa) 2 -flash framework . The SWE are a simplification of the incompressible Navier-Stokes equations, under the assumption that vertical scales are negligible compared to horizontal scales. As a result, hydro-static pressure and the conservation of mass and momentum are part of the equations, while vertical velocity and turbulence are neglected. Manning friction is included in the model with an additional source term (Liang and Marche, 2009). To allow accurate and robust wetting and drying, the model uses a Barth-Jaspersen type limiter that guarantees the conservation of steady states (well-balancedness), mass conservation and positivity preservation of the water-depth (Vater and Behrens, 2014;Vater et al., 2015Vater et al., , 2019. The framework sam(oa) 2 -flash simulates hyperbolic PDEs on dynamically adaptive triangular meshes . It is based on the Sierpinsky Space filling curve, which allows cache efficient traversals of mesh-cells and -edges with a stack and stream approach. Cell-level adaptive mesh refinement in every time step and a water-depth based refinement criterion, allow for multiple levels of refinement for the tracking of wave fronts and other areas of interest (i.e., coasts). While tsunami are usually sourced by setting a perturbation to the seafloor and water surface as initial condition, sam(oa) 2 -flash can include the full spatio-temporal evolution of the seafloor displacement in the simulation. sam(oa) 2 -flash has been validated against a suite of benchmarks (Synolakis et al., 2008).

Dynamic Rupture Modeling for Tsunami Initial Conditions
The time-dependent seafloor displacement generated in the dynamic rupture model is used as input for the tsunami model. Since these displacements are written in form of a triangular unstructured grid by SeisSol, rasterization is required to obtain a regular grid that can be read by sam(oa) 2 -flash . The resulting grid is comprised of rectangular cells of size x × y = 500 m × 500 m. The geometric center of each cell is used to sample the triangular grid using a nearest-neighbor approach. Since all our examples share a relatively large source area and a short process-time, compared to the ocean depth (2 km), we here omit corrections required for landslide-induced tsunami (Kajiura, 1963) and neglect the effects of water flow from seafloor to sea-surface (Saito and Furumura, 2009;Wendt et al., 2009). We directly use the seafloor deformation as time-dependent (or static, for comparison) initial sea-surface perturbation.
The seafloor deformation data contains the seismic wavefield and the dominant static displacement, which remains unchanged after dynamic rupture propagation ceased since we do not account for post-seismic relaxation. Saito et al. (2019) shows that seismic waves have an effect on tsunami generation, but do not affect the solution in the far field. However, our earthquake scenarios also include reverberating seismic waves, e.g., trapped within the accretionary wedge. Simulating long enough for such complex seismic waves to stop imprinting on the time-dependent near-source seafloor deformation requires significantly longer simulation times. We here limit computational costs by stopping the DR earthquake simulation after the fault stopped slipping and the dominant static displacement remains constant. We apply the time-dependent seafloor deformation to source the tsunami model during that period and keep the seafloor elevation constant afterwards. In the presence of transient seismic waves, however, this approach may lead to artifacts such as spurious gravity-waves in the tsunami simulation. Thus, we here remove reverberating seismic waves using a filter-based approach on the seafloor perturbation.
We apply a Fourier filter approach that we base on an analytical test-bed in which we can separate the significant frequency-wavenumber coefficients of the permanent displacement from the ones of seismic waves . In the wavenumber representation for the displacement field we can identify the coefficients belonging to seismic waves clearly as radial symmetric lines in the wavenumber space (Supplementary Figure 1). Our analytical testbed confirms that these lines propagate in the frequency-wavenumber domain with the inverse velocity with which seismic waves propagate in time and space. Coefficients of the permanent displacement on the other hand are dominant close to the origin. To erase seismic waves from the seafloor perturbation we design a kernel to zeroout the radial symmetric waves in the frequency-wavenumber representation depending on their velocity. Close to the origin we leave the representation as it is, to keep the dominant coefficients of permanent displacement. To avoid rolling effects, the kernel is smoothed. As a result, seismic waves of chosen frequencies are effectively damped in the time-dependent displacement field (Supplementary Figure 2). The effects of seismic wave damping on the permanent displacement are negligible.

3D Heterogeneous Megathrust Dynamic Rupture Models
In DR modeling, rupture can only propagate across predefined fault interfaces. In distinction, fault geometry spontaneously arises during slip events in the SC model. The fault geometry for the chosen slip event at the coupling time step is shown in Figure 2A. The locations of the highest visco-plastic strain rate represent the fault. A moving average scheme is used to smooth this 2D fault geometry, which is then uniformly extruded alongarc to construct the 3D DR fault plane (van Zelst et al., 2019;Madden et al., 2020). The 3D dynamic rupture fault ( Figure 2C) does not intersect the surface but ends ≈5 km below the seafloor and extends to 93.5 km depth. The fault length in the alongarc y-direction is chosen to agree with the average scaling of the source dimensions of interface and intraslab subductionzone earthquakes with moment magnitude (Strasser et al., 2010). We here aim to model typical size tsunamigenic earthquakes of magnitudes close to M W =8.5, and thus choose a fault width of ≈313 km. The curved fault is characterized by large scale roughness. Noticeable are two topographic highs shown in terms of sharp increases of fault dip in Figure 2D. We use a constant element edge length of 1 km along the fault. We ensure that this element size along the fault is sufficient to capture the median process zone width following all error criteria of Day et al. (2005) in a series of models with different size elements following the analysis in Wollherr et al. (2018). We measure the median process zone width as 1,386 m for the reference model 3B and 1,224 m for model 1B. For a polynomial basis function of order p = 5  estimate the minimum required resolution to be 1.65 elements per median process zone width. A polynomial basis function of p = 5 leads to 6th order numerical accuracy in space and time in the seismic wave propagation solution. We note that each tetrahedral element fault interface is discretized by (p + 2) 2 Gaussian points. In Supplementary Figure 3 we detail two exemplary higher resolution versions (h=500 m, p=5) of the reference model 3B and of model 1B which has the shortest median process zone width. Compared to the coarser resolution of h = 1 km, we quantify the error to be <1% for the maximum vertical seafloor displacement, <6% for the rupture velocity, <5% for the average peak slip rate and <2.5% for the average fault slip which agree with the criteria for well-resolved models by Day et al. (2005). We note that in contrast to the simple depthdependent setup in Wollherr et al. (2018), the peak slip rate and not the average fault slip is most sensitive here.
The 3D model domain extends from x = −657 to x = 1,075 km and y = −1,023 to y = 700 km and to a depth of z = −700 km. The large model size prevents that waves reflected from non-perfect absorbing boundary conditions interfere with the rupture itself. To limit computational cost, the unstructured tetrahedral mesh is statically coarsened. The sea floor is flat in the DR model (but not in the tsunami simulations) and assigned a free surface boundary condition. The mesh size is 7.8 million elements. Each simulation took 1:10 h on 75 nodes, with 48 Intel Xeon Skylake cores each, of Supermuc-NG.
To analyze the effects of material directivity and complex initial conditions, as well as uni-and bilateral rupture behavior in the DR models, and the resulting tsunami, we vary the hypocenter location along arc and depth. We choose fault locations which are close to failure in the 2D SC model, at 30, 40, and 45 km depth ( Figure 2B). To analyze shallower earthquake nucleation locations, we additionally nucleate 3D DR scenarios at 25 km depth ( Figure 2C). The reference dynamic rupture model (Model 3B) is nucleated at the origin location of the slip event in the SC model, that is, at a depth of 40 km and at the center of the fault along strike (x = 267.25 km, y = −156.5 km). The static friction coefficient is reduced to µ s = 0.019 within a patch of radius 2 km, which represents the minimum value in the SC model within the nucleation area ( Table 1). For all earthquake model families, the assigned nucleation parameters, that is locally reduced static friction coefficients µ s and patch radii r, are listed in Table 1. To evaluate effects in lateral direction, we move the hypocenter from y = −156.5 km (fault width center) to y = −78.25 km (25% of the fault width) and y = −234.75 km (75% of the fault width), exploring observational and statistical inferences of large ruptures being predominantly unilateral or bilateral (McGuire et al., 2002;Mai et al., 2005) (see Figure 2C). We note that larger nucleation energy (larger radii, lower local strength) are required in fault regions further from failure than the reference model (hypocentral depth of 40 km).

Tsunami Model Setup
The tsunami modeling area extends from x = −600 to x = 600 km and from y = −750 to y = 450 km, the ocean depth being at a constant 2 km. A linearly sloping beach is placed with its toe at x = 500 km with an inclination of 5%, which results in the coastline being located at x = 540 km in most models (see Figure 3, left). We additionally analyze the inundation behavior along a more realistic coastline in one model (model 4D). To this end a nonlinear coastline is included in model 3B (Figure 3, right). We adapt the coastal topo-bathymetry of the Okushiri benchmark (Yeh, 1996;Honal and Rannabauer, 2020) and stretch it along the full y-direction model length. Every tsunami is generated using the time-dependent (dynamic) seafloor displacements generated in the 3D DR models from 0 to 200 s. Depending on the used DR model setups, different tsunami modeling refinement levels and output configurations are required. The minimum spatial resolution in sam(oa) 2 -flash is defined as x = domain-width · (1/2) d/2 . We use a minimum refinement level of d = 18 for all tsunami simulation runs, yielding a minimum spatial resolution of x = 2.34 km. To obtain detailed inundation patterns, we use a refinement of d = 30 ( x = 36.62 m) near the coast and a maximum refinement of d = 26 ( x = 146.5 m) in the remainder of the domain. On SuperMUC-NG, these simulations took 2:43 h on 100 nodes sourced by dynamic displacement. This corresponds to roughly 13,000 CPUh, respectively. Simulation outputs are in general written every 10 s of simulation time. If only sea surface height tracing measurements along a few axes are needed, a run with a maximum refinement of d = 24 ( x = 293.0 m) takes approximately 27 min across 100 nodes (2,187 CPUh) with dynamic displacement. Output of the full wavefield with a maximum of d = 24 took 1:21 h across 32 nodes (2,074 CPUh) for dynamic displacements, writing outputs every 100 s of simulation time between t = 0 and t = 3,000. We note that in all our tsunami models, including a slow "tsunami earthquake, " the ratio of tsunami source width/ (source time × tsunami wave speed), in water depth of 2 km is >>1, indicating that the tsunami does not propagate over the source duration (Abrahams et al., 2021). We find differences of up to 10% in the tsunami system energy balance of kinetic and potential energy when comparing static with dynamic sources and source all tsunami in the following from the dynamic seafloor displacement.

Dynamic Rupture Models
We first investigate the effects of varying earthquake hypocenter locations in complex subduction initialized dynamic rupture models. We vary the hypocentral depth between 25, 30, 40, and 45 km, which resembles one shallow earthquake and three low strength excess regions in the geodynamic subduction and seismic cycle model. It has been inferred that hypocenters of large earthquakes are not arbitrarily distributed across fault planes, but located close by regions of large slip (Mai et al., 2005). Additionally, large subduction zone events may propagate preferably unilateral (McGuire et al., 2002) or bilateral (Mai et al., 2005) and rupture dynamics of subduction earthquakes may be significantly affected by bimaterial contrasts (Ma and Beroza, 2008). In our 12 scenarios listed in Table 1, different hypocenter depths lead to pronounced differences in dynamic rupture propagation. The models differ in their propagation direction,  Table 2 summarizes all rupture characteristics.
In model family 1 (shallowest hypocenter) nucleation is initiated at a depth of 25 km which is located above topographic high 1 on the megathrust. We observe complex dynamic rupture behavior, including supershear transition at topographic high 2, localized high slip rates within the fault depression and reactivation of fault slip at a late stage. Rupture propagates at low slip rates before reaching the edge of the first topographic high, where the slip rate increases. As the main rupture front hits the second topographic high (10 s), supershear rupture initiates in downdip direction (Figure 4).
For model family 2, we observe complex downdip propagating dynamic rupture behavior with supershear rupture being triggered at the second topographic high. The hypocenter is located at the edge of topographic high 1 (Figure 2). As the rupture front passes topographic high 1, slip rate increases, similar to model family 1. Supershear rupture is triggered in downdip direction at 15 s simulation time (Figures 4, 5).
The nucleation location in model family 3 corresponds to spontaneous failure in the 2D SC model. It is located on the lowermost point of the geometric depression of the fault, at 40 km depth. Slip evolves circularly and propagates away from the hypocenter. After 8 s simulation time, supershear rupture initiates in updip direction at topographic high 1 (see Figures 4,  5).
For model family 4, the hypocenter was placed at a depth of 45 km, which lies at the edge of the second topographic high. A circular rupture front evolves during 9 s simulation time. Supershear rupture arises in updip direction after 13 s, when the rupture front hits the first topographic high (Figures 4, 5).
In all 12 models, rupture fronts reach the lower limit of the seismogenic zone at x = 282.25 km after passing the second topographic high. In the SC model, ductile behavior begins to dominate here which is expressed as strength increase in the DR model. Rupture propagation is spontaneously arrested at this depth. Thus, in all 12 models, slip stops at the same depth. At the lateral edges of the prescribed megathrust fault, rupture is geometrically stopped, and no tapering of initial stress or strength is applied. Close to the surface, in the shallowest part of the fault, the sedimentary region allows for small slip while smoothly stopping rupture.
Bilateral rupture evolution appears to be symmetrical for centered hypocenters, despite bimaterial contrasts above and beneath the fault potentially affecting strike-slip faulting contributions (Harris and Day, 2005). For asymmetric hypocenter locations (models A and C of each model family), unilateral dynamic rupture propagation is identical to bilateral models within the first seconds until stopping phases from the closer fault edge affect rupture dynamics. In both cases (A and C), rupture dynamics appear minor symmetric with slight bimaterial effects (see rupture characteristics summarized in Table 2). Larger slip is accumulated at the respective far side of the fault (Figure 6). The highest absolute slip is observed in the models with hypocenter locations at y = −234.75 km (C models) and the lowest absolute fault slip is observed for models with hypocenters placed in the center of the fault (B models).
In all 12 models we observe localized weak reactivation of slip after approximately 100 s simulation time due to dynamic triggering caused by trapped waves. Waves are trapped within the accretionary wedge between the uppermost part of the fault and the surface until the end of the simulation (200 s). They are reflected at the free surface boundary and propagate back to the fault, which leads to very small amounts of shallow slip in the order of centimeters. As noted before, in the tsunami linking step any artificial contribution of these waves to seafloor displacements will be filtered.
Across all 12 models, the highest peak slip rate (PSR) occurs on the lower part of the fault, inside the geometric depression at ≈40 km depth, and spread out in along-strike direction (Supplementary Figure 4, Table 2). Model family 2 produces the highest PSR (model 2A, 27.29 m/s). Deeper earthquakes produce lower peak slip rates, such that the earthquakes with the deepest hypocenter (model family 4) exhibits the lowest PSR (model 4A, 23.06 m/s). Along-arc, PSR first increases while rupture propagates away from the hypocenter. Then it decreases due to dynamic interaction with the free surface and other fault edges. At the hypocenter, the PSR is low.
The overall highest absolute slip of 34.31 m is observed for model 4C (see Table 2) and the lowest absolute fault slip is observed for model 1B (32.40 m). For models with the same hypocentral depth, the maximum accumulated fault slip is consistently observed when the hypocenter location is located laterally at y = −234.75 km. At the same time, these earthquakes show relatively low peak slip rates. The lowest maximum absolute fault slip for models with the same hypocenter depth is observed for a laterally centered hypocenter. Moment magnitudes vary from M W = 8.87 (model 3B) to M W = 8.89 (model 4B, Table 2).

Tsunami Simulations
At 200 s simulation time, the filtered vertical sea-surface uplift has a maximum of ≈4 m which is located above the buried dynamic rupture fault plane for all 12 models (see Supplementary Figure 5). The sea surface uplift and subsidence reflect the patterns of accumulated slip in Figure 6. Thus, despite the stark dynamic differences in rupture dynamics between the 12 models, including supershear rupture evolution in upor downdip direction, the static sea surface disturbance is nearly the same. For models of one family, we see lateral differences in the spatial extend of the sea surface uplift which are related to laterally varying hypocenter locations. For models with asymmetric on-fault slip distributions, we observe the same pattern in the sea surface uplift.
Supplementary Figure 6 illustrates the dynamically sourced tsunami propagation toward the simulation domain boundaries. After 2,300 s simulation time, the tsunami arrive at the coast with wave heights of up to ≈5.5 m. At the coast, we observe differences in the tsunami arrival times for laterally varying hypocenter locations of up to 100 s (Figure 7). The time delay between the tsunami waves caused by an earthquake with a centered hypocenter (B) and an earthquake with a hypocenter located at y = −78.25 km (A) are always some tens of seconds higher than the time delay between events with a centered hypocenter (B) and a hypocenter location of y = −234.75 km (C) (see Supplementary Figures 7-9). We observe only insignificant differences in the tsunami arrival time at the coast with varying hypocentral depths (Supplementary Figures 10-12). Tsunami that were generated by deeper earthquakes arrive few seconds later than those being generated by shallower ones. Figure 7 shows the sea surface height (ssh) of the tsunami when arriving at the coast. We observe non-symmetric differences in coastal ssh in dependence of earthquake alongstrike hypocentral location. A maximum wave height of ≈5.5 m can be observed in all 12 simulations. The difference in the tsunami height ( ssh) of models with a centered hypocentral location (B) and a hypocenter located at y = −78.25 km (A) present higher values of approximately 0.25-0.4 m than the ssh of earthquakes with a centered hypocenter location (B) and earthquakes located at y = −234.75 km (C) ( ssh is ≈ 0.1 m). This agrees with the differences in tsunami arrival at the coast. For larger time delays we observe larger differences in the tsunami height accordingly. In summary, comparing all model families 1A-C, 2A-C, 3A-C, and 4A-C, the largest difference of 6 cm in coastal sea-surface heights can be observed between the shallowest earthquakes (models 1A-C) and earthquakes nucleating at 40 km depth (models 3A-C) (see Supplementary Figures 13-18).
We calculate the potential energy transferred by the earthquake rupture to the sea surface (the "tsunami potential FIGURE 7 | (Top) Comparison of sea surface height (ssh) for tsunami fronts arriving at the coast for the models 2A-2C (upper row). The difference in arrival times between the models is shown as ssh (lower row). (Bottom) Inundation comparison for models 3A-3C. The green to white color scale shows the time delay for the tsunami fronts arriving at the coast (upper row). The difference in inundation between respective models is plotted as t (lower row). energy") (Melgar et al., 2019;Crempien et al., 2020): where η is the vertical, static sea-floor deformation of the DR model (corresponding to the sea surface heights after 200 s shown in Tables 2, 3), with a water density of ρ=1,000kg/m 3 and the gravitational acceleration being g=10m/s. The deepest earthquake (model 4B) produces the highest tsunami potential with ≈3,150 TJ which marginally differs from the tsunami potential of model 3B ( Table 2). The tsunami resulting from shallower hypocenter depths (models 1B and 2B) have slightly smaller tsunami potentials of 3,078 and 3,062 TJ. Within model families 1 and 2, a centered hypocenter location leads to the highest tsunami potential, while for model family 3 and 4, the centered hypocenter causes the smallest tsunami potential. Within each model family, tsunami simulations with a hypocenter location at y = −78.25 km (A) always present lower tsunami potentials than models with hypocenter locations at y = −234.75 km (C).

Tsunami Simulation With Complex Coastal Topo-Bathymetry
We perform an additional tsunami simulation, model 3D, which is adapting model 3B by replacing the linearly sloping beach with a complex coastline (see Figure 3, right). We adapt the coastal topo-bathymetry of the Okushiri benchmark (Honal and Rannabauer, 2020) and stretch it along the full y-direction model length (section 3.2). The resulting sea-surface height and inundation area are displayed in Figure 8. Tsunami inundation is observed along the entire length of the coast. While the coastal parts at the far-ends of the domain are hit by a shallow wave of approximately 1 or 2 m, its central part is hit by high wave amplitudes of up to 8 m. In this model, the tsunami amplitudes  are ≈2.5 m higher, compared to absolute sea-surface heights of ≈5.5 m for the linearly sloping beach.

Simulations With Increased Fracture Energy and Poisson's Ratio
We here analyse two additional dynamic rupture models based on reference model 3B varying on-fault or off-fault rheology. Firstly, we triple the critical slip weakening distance D c from 0.1 to 0.3 m (model 5) which triples fracture energy and generates a "tsunami earthquake." Secondly, Poisson's ratio is increased from ν=0.25 to ν=0.3 (model 6) everywhere in the domain. Both models result in shallow slip about twice as high as the reference model 3B (see Table 3). Table 1 shows the adapted nucleation characteristics that were necessary to initiate rupture on the fault. We do not further decrease the static friction coefficient, but increase the nucleation radius from 3.5 to 10.0 km (model 5). In model 6 a smaller nucleation area of only 1.8 km is sufficient. For the high fracture energy model 5, rupture dynamics evolve very differently to model 3B, specifically at much lower rupture velocities of max. 1,352 m/s. There is no supershear rupture triggered during the entire simulation time. In model 6 the dynamic rupture evolution is similar to model 3B and supershear rupture evolves in updip direction after 9 s. For both models, 5 and 6, trapped waves are observed until the end of the simulation (200 s) dynamically interacting with the shallow part of the fault. Figure 9 displays the rupture characteristics of model 3B, model 5, and model 6. Compared to the reference model 3B, both adapted models accumulate large shallow slip. The maximum and average fault slip of models 5 and 6 are about twice as high as in the reference model 3B (see Table 3), which reflects in increased earthquake magnitudes of M W =9.03 (model 5) and M W =9.04 (model 6). Also, their dynamic stress drops are higher and the maximum vertical dynamic seafloor displacement is increased by up to ≈2 m. Model 5 shows a much smaller PSR (12.59 m/s) and rupture velocity (1,352 m/s) then model 3B (24.12 and 2,124 m/s), while the PSR (24.80 m/s) and rupture velocity (1,936 m/s) of model 6 are similar to the reference model. The maximum PSR for all three models is always observed at the same depth which is located within the fault depression at ≈ 40 km depth. Figure 10, top, compares the sea-surface height of models 5 and 6 after 200 s, corresponding to the end time of the DR simulation. The overall tsunami waveforms in models 5 and 6 appear to be broader than in model 3B (Supplementary Figure 7) and the trajectories in Figure 10, bottom, show much higher tsunami amplitudes. After 1,400 s simulation time, the wave amplitudes of model 3B reach extrema of +2 and −3 m, whereas the tsunami in models 5 and 6 reach values of +2 and −5 m. After 2,300 s simulation time, the wavefronts hit the coast and result in maximum tsunami heights of over 7.5 m. This is ≈2.0 m higher than in the reference model 3B. An important difference between models 5 and 6 is the difference in rupture speed. While model 6 produces supershear rupture, the overall rupture velocity of model 5 is ≈1,352 m/s. Thus, although the model 5 earthquake scenario has a slightly smaller stress drop and magnitude than model 6, it produces the highest tsunami amplitudes.
As the waves hit the coast there is a time delay of 100 s and a difference in tsunami height of ≈0.5 m between models 6 and 5 (see Figure 11). Due to the overall higher tsunami waves of models 5 and 6, the water inundates further on-shore and reaches higher distances from the coast than for model 3B. Model 5 has the highest tsunami potential with roughly 6,950 TJ and model 6 has a tsunami potential of 6182.48 TJ ( Table 3). The tsunami potential of model 5 is twice as high as the one of reference model 3B. Figure 12, top, shows snapshots of the unfiltered DR seafloor displacements of models 1B (supershear rupture in downdip direction), 3B (supershear rupture in updip direction), and FIGURE 9 | Fault slip, peak slip rate, and rupture velocity for the dynamic rupture models 3, 5 (increased fracture energy), and 6 (higher Poisson's ratio) after 200 s at the end of the DR simulation. 5 (tsunami earthquake, no supershear rupture) after 100 s simulation time. Localized, sharp uplifting fronts are visible in the dynamic displacement off-set from the centrally located hypocenter overprinting the static deformation signal. The ocean response recorded within the source region during the tsunami generation process of all 3 models (Figure 12, bottom) reflect the seismic, and near-field displacements around a rupture front at 100 s simulation time. The time series shown is recorded at x = −100 km and y = −150 km, which is well inside the DR modeling domain.

Dynamic Effects During Tsunami Generation by Supershear and Tsunami Earthquakes
In our models, co-seismic ocean response phases appear for supershear earthquakes as well as for the "tsunami earthquake" propagating at sub-Rayleigh speed during the duration of earthquake slip and within the DR model, i.e., for the dynamic tsunami generation process. We here do not observe a faster, instantaneous supershear mach cone ocean response signature (e.g., identified in Elbanna et al. (2020) for strike-slip events). Additionally, due to our filtering approach, trapped seismic waves are effectively damped, also inside the near-fault region once rupture stops. We then do not observe dynamic phases propagating (Supplementary Figures 5, 6).

Simplifying Model Assumptions
In this study, we link 3D dynamic rupture initial conditions to a chosen slip event in a 2D long-term geodynamic subduction and seismic cycle. The linked initial conditions include a curved, blind fault geometry, spatially heterogeneous fault stresses, strength, and material properties. The SC 2D material properties, fault geometry, as well as stresses and strength are extruded into the third dimension without adding additional along-strike variability. While limiting complexity, we can in this manner isolate sensitivities, e.g., of hypocentral location, and their effects on rupture dynamics and tsunami generation, propagation and inundation.
In linking from the SC to the DR model, we adopt several simplifying assumptions to bridge the incompressibility and visco-elasto-plastic, plane-strain conditions of the subduction model to the compressible, elastic conditions of the earthquake model. The resulting 3D dynamic rupture is linked with the tsunami model through the time-dependent seafloor displacements, following the same methods as detailed in Madden et al. (2020). All resulting rupture models are characterized by uniformity along strike, rendering the dynamic earthquake source differences to be associated with along-dip variability in rupture dynamics controlled by the nucleation position relative to the two bands of topography on the megathrust. Future work may extend our analysis to include along-arc variations in stress, strength, rheology, and geometry across the megathrust as suggested by detailed imaging of fault geometry and shown to impact earthquake and tsunami dynamics (e.g., Galvez et al., 2019;Ulrich et al., 2020). The adopted "bumpy" fault geometry develops self-consistently with stress, strength and rheology in the SC model. Its along-depth dimensions aligns with conceptual asperity models governing megathrust fault slip. Future work may study the dynamic effects of subducted seamounts (e.g., Cloos, 1992) ensuring to adapt initial stresses consistently.
We use fully elastic material response in combination with linear slip weakening friction. The complexity of the DR model could be increased by including more complex physics, such as FIGURE 11 | (Top) Comparison of sea surface height for tsunami fronts arriving at the coast for the models 3B, 5, and 6. The difference between the models is shown as ssh. In contrast to Figure 11 the x-axis (distance from coast) indicates higher values. This is due to the greater inundation area of model 5 which exceeds 161 m distance. (Bottom) Inundation comparison for models 3B, 5, and 6. The green to white color scale shows the time delay for the wave fronts arriving at the coast. The difference between the models is plotted as t with a blue (negative values) to red (p color-scale, respectively. rapid velocity weakening rate-and state friction (Ulrich et al., 2019a), thermal pressurization of pore fluids , or off-fault plasticity Ma and Nie, 2019).
In the linking step from the DR model to the tsunami model, we filter the seafloor displacements using a spatialtemporal Fourier-transform (section 4.5). Detailed analysis of the effects of relatively small coseismic phases on tsunami genesis, propagation, and inundation is here challenging due to the filter we apply. Future studies may use the unfiltered seafloor displacement as input to the tsunami model to analyse the fully dynamic interaction of the seafloor movements with the tsunami and inundation dynamics. sam(oa) 2 -flash 's hydrostatic shallow water tsunami model enables the simulation of tsunami genesis, propagation, and inundation at the coast. The approach is limited by the assumption of long wavelengths. Additionally, it does not take the interaction of wind with the water interface into account. Overall, the usage of the shallow water equations might overestimate wave amplitudes.
Our DR models do not account for seafloor bathymetry. A more realistic bathymetry would translate the horizontal earthquake motion into vertical displacements, such that the tsunami amplitude might be amplified (Tanioka and Satake, 1996;Lotto et al., 2018;Saito et al., 2019;Ulrich et al., 2019a). Recent work in van Zelst et al. (2021) complement this paper by using the exact subduction bathymetry and topography arising from the 2D seismic cycle model and include it in 2D DR simulation. The resulting seafloor displacement is linked to a one-dimensional shallow water tsunami model. The resulting maximum tsunami amplitude is 6.5 m, larger than our comparable 3D reference model 3B. It is here left for future work to combine their findings with our 3D sensitivity analysis.
The computational costs of each of the presented 15 linked scenarios (see sections 3.1 and 3.2) is well within the scope of the allocation typically available to users of supercomputing centers. While hundreds of such simulations are readily possible, fully physics-based dynamic rupture models rather complement than replace cheaper (e.g., kinematic) source descriptions used for millions of PTHA forward models. Specifically, for narrowing down the high-dimensional and often non-unique source parameter space in conjunction with observational or longterm modeling constraints and for sensitivity analysis of other parameters influencing rupture behavior and tsunami generation and propagation.

Hypocentral Depth and Up-Dip vs. Down-Dip Supershear Rupture Propagation
We vary the hypocenter location across four depths (25, 30, 40, 45 km) to study the effects on rupture dynamics and tsunami evolution, propagation, and coastal inundation. Earthquakes with shallower hypocentral depths (25 and 30 km depth) generally generate lower slip than earthquakes with deeper hypocenters (40 and 45 km depth). The lower accumulated on-fault slip for events with shallower hypocenters leads to comparably lower vertical seafloor displacement. In our SC initialized models, however, these differences are relatively minor and have small impact on tsunami generation, propagation and inundation, in contrast to what is typically observed from observational data (Bilek and Lay, 2018).
For simulations with shallower hypocenter locations, we observe supershear rupture propagating in the downdip direction, while hypocenters located at 40 and 45 km depth lead to supershear rupture in the updip direction. In either case, the slab geometry (topographic highs) and rheology influences the nucleation and direction of supershear rupture propagation significantly. In all our simulations, supershear rupture initiates when the rupture front hits a topographic high on the fault plane. The first topographic high (1) is related to a pile up of subducted sediments. The weaker material at the depth of the topographic high 1 might facilitate supershear rupture. Supershear rupture is also triggered by the second topographic high highlighting the complex dynamic effects of the long-term self-consistently developing fault roughness, stress, and rheology heterogeneities. The direction of supershear rupture propagation is determined by whether the updip or downdip rupture front interacts with the rough fault geometry (Bruhat et al., 2016;Bao et al., 2019;Tadapansawut et al., 2021). Megathrust supershear rupture is challenging to identify observationally but has been suggested in back-projection analysis of the Tohoku-Oki megathrust earthquake (Meng et al., 2011) and was analyzed in Cascadia 2D dynamic rupture models by Ramos and Huang (2019) who observe down-dip supershear rupture propagation near the ETS region caused by normal stress gradients (similar to the conditions governing our topographic highs). Independent of where the earthquakes nucleate, the highest peak slip rate is consistently observed at the same location on the fault: inside the depression that separates the two local fault topographic highs, although the intensity of the slip rate decreases with increasing hypocentral depth. The calculated tsunami potential energy varies in the range of E T ≈78 TJ for earthquakes nucleating at different depths. This is caused by a difference in the maximum seafloor displacement of approx.

Bilateral vs. Unilateral Rupture on a Complex Bimaterial Megathrust
To study unilateral vs. bilateral rupture effects on rupture dynamics (Hirano, 2019;McGuire et al., 2002) and tsunami generation, we shifted rupture nucleation from a centered hypocenter to 25% (A) and 75% (C) of the fault width (y = −78,25 and y = −234,75 km). Overall, differences are small. However, we note that earthquakes with a centered hypocentral location (B) consistently produce the highest stress drops and lowest accumulated slip on the fault, which lead to the lowest vertical seafloor displacement. The simulations with a nucleation patch at 75% (C) of the fault width show the highest accumulated slip for earthquakes at the same depth, but relatively low PSR. Due to the high amount of accumulated slip, we would expect the seafloor displacement to be the highest as well, but this is not the case for all model families, due to bimaterial effects (Brietzke et al., 2009). These effects lead to (minor) differences in amplitudes of the resulting tsunami at the coast and the (short) time delay between the arriving wave fronts. The higher sea surface uplift is also seen in the trajectories of the different models (Supplementary Figure 7). We note again, that in difference to Madden et al. (2020), we here use the 2D geodynamic SC model developed in van Zelst et al. (2019) which has twice larger shear moduli. All our earthquakes show a higher stress drop of roughly 5 MPa compared to 2.2 MPa in Madden et al. (2020), whereas the rupture velocity remains the same (≈ 2, 100 m/s). They also show a comparably lower maximum slip of 32.8 m  observe a slip of 42.2 m) on the fault plane, which is visible in the lower seafloor displacement as well (4.60 m compared to 28.1 m). Overall, our earthquake and tsunami scenario agree with observational scaling, resulting in a typical tsunami generating earthquake magnitude of M W = 8.8 and a tsunami wave height of roughly 5.5 m at the coast. For all 12 models with hypocenter variations, we find that high earthquake magnitudes (or high fault slip) correspond to the highest seafloor displacements and result in greater tsunami heights, as expected.

Comparison of Tsunami Behavior for Linear and Complex Coastline
To analyze the effect of coastal complexity on inundation, we included a non-linear coastline in model 3B (see Figure 3). The results show that the complex and more realistic setup yields higher tsunami amplitudes (up to 8 m, Figure 8) than the model with a linear beach (sea-surface height of up to 5 m, Figures 7,  11). The overall distribution of sea-surface heights along a nonlinear coast is much more complex. Between y = 50 and y = −160 km, the distance between fault and coast increases. Beyond y = −160 km until y = −350 km this distance decreases. While the part between y = 50 and y = −160 km is hit by tsunami heights of up to 8 m, the sea-surface heights at the coast between y = −160 and y = −350 km are ≈2 m lower. This effect may be enhanced when combining a complex coastline with lateral varying earthquake source characteristics. Even though the waves of both, model 3B and the scenario with the complex beach arrive nearly simultaneously at the coast, they need ca. 1,000 seconds longer to reach the farthest onshore point in the non-linear case.

Large Shallow Slip
Most earthquakes of high magnitudes tend to have a large stress drop, accompanied with a high radiation energy (Festa et al., 2008). Venkataraman and Kanamori (2004) state that this is different for a special class of megathrust events, "tsunami earthquakes, " which have a comparably small radiation efficiency. We here model this dynamically by increasing the critical slip distance, the amount of slip over which the static friction coefficient in the dynamic rupture model decreases to the dynamic friction coefficient. Polet and Kanamori (2000) state that tsunami earthquakes tend to rupture the shallow portions of a fault, which results in a large amount of shallow slip. The increased D c value in model 5 leads to a lower rupture velocity of 1,352 m/s compared to the rupture speed in the initial model 3B of 2,124 m/s. The mean velocity of 1,352 m/s is lower than 80% of the S-wave speed (van Zelst et al., 2019), which specifies this earthquake by definition as a tsunami earthquake (Kanamori, 1972;Bilek and Lay, 1999). As the rupture propagates updip, we observe high slip on the shallow fault. The slip is twice as high as in the initial model 3B, which leads to a higher max. vertical seafloor displacement (1.5 times the displacement of 3B) and a higher tsunami height (see Figure 11). This behavior is reflected in the tsunami potential of model 5, which is double the tsunami potential in model 3B.
In model 6, we increase the Poisson's ratio from ν = 0.25 to ν = 0.3. The increase in Poisson's ratio results in a reduction of the critical maximal shear stress on the fault (Xie et al., 2009). This means that rupture initiation is facilitated. Whereas, van Zelst et al. (2019) observe a decrease in slip with increasing Poisson's ratio, we find higher slip on the fault, together with a higher stress drop; model 6 has also double the amount of slip of model 3B, just as the "tsunami earthquake" model 5. This causes a higher vertical seafloor displacement and thus a higher tsunami. At the same time the rupture velocity remains nearly the same as in model 3B. Due to the high magnitude and accumulated shallow slip, the tsunami potential of model 6 is nearly twice as high as in model 3B. The differences in rupture dynamics between van Zelst et al. (2019) and our model 6 from varying the Poisson's ratio is due to the different choices in linking frictional parameters (we here assume constant D c ) which leads to differences in stress drop. While the stress drop in van Zelst et al. (2019) remains the same with increasing ν, we observe an increase of 1.08 MPa.
The main difference between model 5 and 6 are summarized in Tables 2, 3. In model 6, the peak slip rate and stress drop get amplified, leading to a similar rupture behavior than in model 3B with a higher magnitude and absolute slip, resulting in higher tsunami amplitudes. In model 5, the rupture velocity gets reduced and the rupture characteristics change. Although model 5 produces an earthquake with a slightly lower magnitude and stress drop, it produces a seafloor displacement that is 46 cm higher than for model 6. The tsunami earthquake (model 5) generates the highest tsunami amplitude of all these models and consequently the greatest inundation area at the coast, while the waves arrive significantly later due to the lower rupture velocity. The effect of a Poisson's ratio increase in model 6 is not quite as large as the change of the rupture dynamics and tsunami generation in model 5. While in model 5 rupture propagates at sub-shear speeds, supershear rupture still evolves in model 6. Nevertheless, an increasing critical slip weakening distance D c just as a change in the material properties in the earthquake rupture model can drastically change the rupture dynamics and influence tsunami generation and propagation (see Figure 11). We note that for future analysis of the effects of enhanced shallow slip such as occurring in both models 5 and 6, it will be crucial to combine our analysis with 3D non-constant water depth in the source region, since a realistic subduction zone geometry (van Zelst et al., 2021) will place shallow slip in deeper water and cause additional effects due to horizontal motions and steep topography contrasts.

CONCLUSION
We investigate the influence of hypocentral depth, rupture propagation direction and bimaterial effects, as well as the influence of fracture energy and Poisson's ratio on rupture behavior and tsunami generation and propagation. We analyse 15 subduction-initialized 3D dynamic earthquake rupture tsunami propagation and tsunami run-up scenarios. We vary the hypocentral depth between 25, 30, 40, and 45 km, which resembles four low strength excess regions in the geodynamic subduction and seismic cycle model. In all models, supershear rupture is triggered once the earthquake rupture front crosses one of two distinct topographic highs in the megathrust geometry, which are related to sediment subduction on geodynamic time scales. Earthquakes from shallow hypocenters exhibit supershear rupture in the downdip direction, whereas supershear rupture propagates updip for earthquakes that nucleate deeper. Albeit dynamic earthquakes differ (rupture speed, peak slip-rate, fault slip, bimaterial effects), the effects of hypocentral depth on tsunami dynamics are negligible. Earthquakes with deeper hypocenters accumulate higher slip during up-dip rupture compared to shallow hypocenters, in which rupture mainly propagates downdip. Larger fault slip correlates with larger vertical seafloor displacement by up to 13 cm, which is reflected in the tsunami potentials. These tendencies barely affect the tsunami run-up behavior at the coast, where the maximum difference in tsunami height is only a few centimeters and the wave arrival times vary by few seconds.
Lateral hypocenter variations lead to small effects such as delayed wave arrival of up to 100 s and differences in tsunami amplitude of up to 0.4 m at the coast. To study unilateral vs. bilateral directivity effects on dynamic megathrust rupture, tsunami generation, propagation, and inundation, we varied the hypocenter location along-strike at all of above depth locations. We find that the highest fault slip is always observed for unilateral rupture with hypocenters located at 75% of the fault width (at y = −234.75 km), whereas a centered rupture initiation leads to purely bilateral rupture including the lowest dynamically accumulated slip. In between models of one model family, fault slip varies up to ≈1.5 m. We find only minor bimaterial effects; models with hypocenters located at 25% of the fault width mostly mirror those with hypocenters at 75% of the fault width.
We dynamically generate a "tsunami earthquake" by increasing the critical slip distance, and thus increasing the amount of fracture energy and decreasing radiation efficiency of the bilateral, 40 km deep dynamic earthquake rupture model. This results in lower rupture velocities (average rupture velocities in model 5 are 64% of those in model 3B) and doubles the amount of on-fault slip which is then, in contrast to the initial model, concentrated on the shallow part of the fault. This leads to a ≈2 m higher vertical seafloor displacement and a ≈2 m higher tsunami amplitude at the coast. Increasing Poisson's ratio has a similarly large effect on shallow fault slip, but less on tsunami height and run-up. Increasing ν from 0.25 to 0.3 doubles the amount of fault slip and favors shallow slip, leading to a vertical seafloor uplift of ≈6 m, which is an increase of 1.5 m and a difference of up to ≈1.5 m in coastal tsunami amplitudes.
Our sensitivity analysis based on 15 physics-based linked earthquake and tsunami and inundation models for a generic megathrust setting can provide building blocks toward dynamic rupture modeling complementing Probabilistic Tsunami Hazard Analysis (PTHA). Virtual laboratories, such as we use here, using computationally efficient and open source earthquake and tsunami computational models enable hypothesis testing and physics-based plausibility assessment of linked tsunami and earthquake models of varying complexity.

AUTHOR CONTRIBUTIONS
SW further post-processed the geodynamic seismic cycle (SC) data, designed the DR models, analyzed the DR data and visualized it, and wrote the initial draft of the study. A-AG initiated the study, revised the draft, supervised SW, and acquired the financial support for the ChEESE project leading to this publication. MS performed tsunami simulations and contributed substantially to the manuscript. EM provided tools (used in Madden et al., 2020) to analyze the data of the SC model and to format the data in a way such that they can be used as input for the dynamic rupture model (DR), helped designing the DR models, and critically reviewed the manuscript. IvZ provided firstly postprocessed data output from SC model, contributed to and critically reviewed the manuscript. LK co-supervised MS and contributed a conversion tool from dynamic rupture to tsunami models and provided critical feedback on the manuscript. YvD provided critical comments on the manuscript. LR designed the tsunami studies together with SW, A-AG, and MS analyzed and visualized the tsunami model output, analyzed and postprocessed the DR models and co-supervised MS. All authors contributed to the article and approved the submitted version.

FUNDING
This research has been supported by the European Union's Horizon 2020 Research and Innovation Programme under the projects ChEESE, grant no. 823844 and TEAR, grant no. 852992. IvZ was funded by the Royal Society (UK) through Research Fellows Enhancement Award RGF\EA\181084. Computing resources were provided by the Institute of Geophysics of LMU Munich and the Leibniz Supercomputing Centre (projects no. pr63qo, pr45fi, and pn68fi).