Earthquake Rupture on Multiple Splay Faults and Its Effect on Tsunamis

Detailed imaging of accretionary wedges reveals splay fault networks that could pose a significant tsunami hazard. However, the dynamics of multiple splay fault activation during megathrust earthquakes and the consequent effects on tsunami generation are not well understood. We use a 2‐D dynamic rupture model with complex topo‐bathymetry and six curved splay fault geometries constrained from realistic tectonic loading modeled by a geodynamic seismic cycle model with consistent initial stress and strength conditions. We find that all splay faults rupture coseismically. While the largest splay fault slips due to a complex rupture branching process from the megathrust, all other splay faults are activated either top down or bottom up by dynamic stress transfer induced by trapped seismic waves. We ascribe these differences to local non‐optimal fault orientations and variable along‐dip strength excess. Generally, rupture on splay faults is facilitated by their favorable stress orientations and low strength excess as a result of high pore‐fluid pressures. The ensuing tsunami modeled with non‐linear 1‐D shallow water equations consists of one high‐amplitude crest related to rupture on the longest splay fault and a second broader wave packet resulting from slip on the other faults. This results in two episodes of flooding and a larger run‐up distance than the single long‐wavelength (300 km) tsunami sourced by the megathrust‐only rupture. Since splay fault activation is determined by both variable stress and strength conditions and dynamic activation, considering both tectonic and earthquake processes is relevant for understanding tsunamigenesis.

Dynamic rupture modeling is a useful tool to understand the role of splay faults in rupture dynamics (e.g., Aslam et al., 2021;DeDontney & Hubbard, 2012;DeDontney et al., 2011;Kame et al., 2003;Lotto et al., 2019;Tamura & Ide, 2011;Wendt et al., 2009). These studies show that parameters such as the initial stress, branching angle, frictional properties, strength of the accretionary wedge, and material contrasts along the megathrust affect splay fault rupture. Several coupled models have been employed to solve for splay fault rupture dynamics and tsunamis sequentially or simultaneously Lotto et al., 2019;Li et al., 2014;Wendt et al., 2009;Ulrich et al., 2022).
Dynamic rupture models of branching faults typically use simple, planar fault geometries, even if observed splay fault geometries are more complicated (e.g., Collot et al., 2008;Moore et al., 2007;Park et al., 2002). Besides that, most dynamic rupture studies include only a single splay fault, which is partly necessitated by the difficulty of modeling fault junctions with numerical methods (e.g., Aochi et al., 2002;Pelties et al., 2014). Another reason for using predominantly simple fault geometries in dynamic rupture modeling up to now is the difficulty in constraining consistent initial stress and strength conditions on complex fault geometries. However, recent studies (Madden et al., 2020;Van Zelst et al., 2019;Wirp et al., 2021) have shown that initial conditions for 2-D and 3-D megathrust dynamic rupture earthquake simulations can be constrained from 2-D geodynamic long-term subduction and seismic cycle models. Indeed, 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 over geological time scales.
To understand the effect of multiple splay fault rupture with non-planar geometries and subduction-initialized stress and strength on the free surface displacements and the ensuing tsunami, we model dynamic rupture constrained by a geodynamic model of long-term subduction and the subsequent tsunami propagation and inundation.

Modeling Approach
We use the modeling approach presented in Van Zelst et al. (2019), where a geodynamic seismic cycle (SC) model is used to constrain the initial conditions of a dynamic rupture (DR) model. We extend this approach by using the resulting surface displacements of the DR model as input for a tsunami propagation and inundation (TS) model. Our modeling framework accounts for the varying temporal and spatial scales from geodynamics to tsunami inundation (see also Madden et al., 2020). We apply this framework to understand the dynamics of splay fault rupture by including six splay fault geometries constrained by the SC model within the DR model setup.

Geodynamic Seismic Cycle Model
We use the same SC model as Van Zelst et al. (2019) which is based on the Southern Chilean subduction zone. We use the output of the SC model as input for the DR model for one event. The SC model solves for the conservation of mass, momentum, and energy with a visco-elasto-plastic rheology (Gerya & Yuen, 2007). It models 4 million years of subduction followed by a seismic cycle phase with a 5-year time step with spontaneous slip events driven by a strongly rate-dependent friction (van Dinther, Gerya, Dalguer, Corbi, et al., 2013) using the seismo-thermo-mechanical (STM) modeling approach (van Dinther, . For a full description and discussion of the methods, we refer the reader to Van Zelst et al. (2019).
We observe widespread visco-plastic shear bands in the sedimentary wedge in the SC model forming during megathrust slip events, which we interpret as faults (Figure 1b). Both in-and out-of-sequence thrusting fault geometries that are typically observed in nature (e.g., Kimura et al., 2007) are present.
For one slip event, we use the output of the SC model as input for the DR model according to Van Zelst et al. (2019). We pick six splay fault geometries according to the highest accumulated visco-plastic strain during and initial sea surface height (blue). The coastline is located at x = 282.25 km. Note that the x-axis differs for each panel depending on the model setup size (trench indicated by the yellow triangle). Also note that the SC model has positive z-axis down, whereas the other two models have positive z-axis up.
Section 3, therein). We model mode II along-dip rupture propagation (e.g., Ramos & Huang, 2019). SeisSol is based on an Arbitrary high-order accurate DERivative Discontinuous Galerkin method (ADER-DG, Dumbser & Käser, 2006) and uses unstructured tetrahedral meshes enabling geometrically complex models, such as branching and intersecting faults (de la Puente et al., 2009;Pelties et al., 2014). The on-fault element edge length is 200 m, which, combined with using basis functions of polynomial degree p = 5 (spatio-temporal order 6 numerical accuracy for wave propagation) results in an effective resolution of 28.6 m through (p + 2) Gaussian integration points on the fault, which is sufficient to resolve the cohesive zone size (Day et al., 2005;Wollherr et al., 2018). At the top of the DR model setup, we employ a free surface boundary condition with topography derived from a third order polynomial approximation of the rock-sticky air (Crameri et al., 2012) interface in the SC model from x = −72.8 to 499.6 km, beyond which we assign constant topography values (Figure 1e). We run the model for 180 s, which ensures smooth coupling to the TS model, as the surface displacements do not vary significantly after that time. To obtain the surface displacements of the DR model, we place 601 virtual seismometers from −100 to 500 km at 5 m below the free surface with a spacing of 1 km to record the velocity field. To optimally capture the surface displacements, we place the seismometers within elements that have a free-surface boundary edge.

Tsunami Propagation and Inundation Model
To model tsunami propagation, we solve the one-dimensional shallow water equations, which consist of the conservation of mass and momentum and consider the hydro-static pressure caused by gravitational acceleration. Recently the more advanced Boussinesq equations have grown in popularity to model tsunami propagation (Spiegel & Veronis, 1960). However for models of the type that we simulate in this work (i.e., a large domain compared to a small wave amplitude) the shallow water equations have been validated and proven to be an accurate model (Carrier & Greenspan, 1958). To solve the non-linear shallow water equations, we employ a first order finite volume scheme (LeVeque, 2002) and we use a well-tested augmented Riemann solver to solve for inundation (George, 2008).
To incorporate dynamic surface displacements, we consider the bathymetry as a constant, defined by the unperturbed topography from the SC model, plus a time-dependent deformation from the DR model that incorporates all effects. Following Abrahams et al. (2020), this approach is sufficient to capture all components of the deformation that contribute to the tsunami. The constant topography from the SC model has an average beach angle of 7.2 ⋅ 10 −6 ( Figure 1e). To compute the seafloor deformation from the DR model, we use the method by Tanioka and Satake (1996), which adds the vertical displacement to a linear approximation of the vertical contribution of the horizontal displacement. We then add the computed seafloor deformation displacements Δb(x, t) from the DR model to the SC model topography. The resulting displacement field contains fast traveling seismic waves, which are radiating from the earthquake source during the DR simulation. Waves are trapped within the sedimentary wedge between the uppermost part of the fault and the surface until the end of the simulation. To avoid imprinting of wave signals on the near-source seafloor deformation, we remove the seismic waves from all displacements used as tsunami sources. To this end, we apply a Fourier filter (Wirp et al., 2021) to the seafloor displacements which removes transient displacements resulting from waves with a ratio of frequency over wave number higher than 300 m/s (Figures S17-S18 in Supporting Information S1).
At this point other approaches additionally account for the energy transfer from the seafloor to the water surface and apply a low pass filter to the horizontal displacements (Kajiura, 1970;Wendt et al., 2009). In our models the source size is relatively large compared to the source duration, so we follow Saito (2013) in ignoring this energy transfer and instead directly adapt the change of the seafloor to the sea surface. We consider a model domain from x = −300 to 500 km, with the initial bathymetry from the SC model ( Figure 1e). We set the coastline at x = 282.25 km to coincide with the downdip limit of the seismogenic zone (Klingelhoefer et al., 2010). This results in a maximum water depth of 4,117 m. To discretize the model, we use 20,000 points, which translates to a uniform spacing of 40 m. We use adaptive time stepping and run the model for a total simulation time of 2 hr with maximum time steps of 0.5 s and minimum time steps of 0.08 s. The time step size is adapted according to the maximum wave speed in the model, which depends on the water column. Close to the coast, the size of the water column reduces to values close to zero, which increases the wave speeds and reduces the maximum admissible time step size according to the Courant-Friedrichs-Lewy condition (Courant et al., 1928). To avoid numerical instabilities, we consider cells with a water column of less than 10 −6 m as dry.

Stress Field and Splay Fault Geometries
The chosen six splay fault geometries that are activated during a representative slip event (Appendix A) result from realistic tectonic loading during retreating subduction on geodynamic time scales (Figure 1). The four shallowest splay faults (SFs) 1-4 are located within the sediments scrapped off from the ocean floor and SF5 follows the contrast in shear modulus between the incoming sediments that make up the accretionary wedge and the sediments of the pre-existing sedimentary wedge (Figure 1d). At the branching point with the megathrust, the largest splay fault (SF6) is initially situated in the weaker incoming sediments, but then travels through the stronger basalt and into the sedimentary wedge sediments. The dips of the splay faults average 24.0° and the branch angles between the splay faults and the megathrust average 14.4° (Table S1 in Supporting Information S1), which is in line with observations (Park et al., 2002) and Mohr-Coulomb theory.
At nucleation, the sedimentary and accretionary wedge are largely under compression with the principal stress direction approximately 22° from horizontal (red bars in Figure 2). This agrees with dynamic Coulomb wedge theory (Wang & Hu, 2006) Figure S7-S10 in Supporting Information S1). However, larger strength excess of 1-6 MPa exists across the large splay fault SF6 and the deeper parts of SF4 and SF5 ( Figure 3, Figures S7-S10 in Supporting Information S1).
The likelihood of fault activation through earthquake rupture can be analyzed through a comparison to theoretical fault growth angles (e.g., Kame et al., 2003). Faults form at an angle to the local stress field, which is generally believed to obey the Mohr-Coulomb failure criterion (e.g., Anderson, 1905;Heidbach et al., 2018;Sibson, 1994). We calculate the Coulomb angle α (white bars in Figure 2b) at which faults theoretically form with respect to the principal stress direction in a compressional stress regime according to (e.g., Choi & Petersen, 2015;Kaus, 2010;Wang & Hu, 2006;Zang & Stephansson, 2010):  in Supporting Information S1 for the failure analysis on the megathrust and the other splay faults. The shallowest part of the fault is at 0 km; the splay fault connects to the megathrust on the right hand side of the figure. Initial shear stress τ, fault yield stress (strength) dr yield , and strength excess dr yield − are shown for the DR model in the fault coordinate system. Frictional regimes dependent on temperature are indicated with corresponding isotherms (solid black lines). Background colors represent the material through which the fault is going: incoming sediments (orange), pre-existing wedge sediments (brown), and basalt (green). Fault strength and therefore stress depends on lithology, which results in stress and strength variations along SF6 which cuts through multiple lithologies. The sharp, localized increase in strength on SF6 is due to a local lack of fluids in the host rock.
where ϕ = tan −1 (μ d ) is the angle of friction with μ d being the dynamic friction coefficient of the sediments. We use the dynamic friction coefficient μ d = 0.105 to calculate the Coulomb angle instead of the static friction coefficient angle μ s , since strain localization forming shear bands in the SC model typically occurs during a slip event. Slip events are characterized by increased slip velocity and therefore reduced effective friction coefficient . This results in a Coulomb angle of −42° with respect to the principal stresses. Throughout the sedimentary wedge, this leads to a Coulomb angle of approximately −20° with respect to the horizontal. The splay fault geometries generally align very well with the theoretical optimal faulting angles ( Figure 2b), indicating that they are favorably orientated for activation during earthquake rupture. Interestingly, the deepest sections of SF1-3 and SF5, where they branch off the megathrust, are not aligned with the theoretical optimal Coulomb angles. Instead, the megathrust aligns with the Coulomb angle near the branching junctions.
During each slip event in the SC model, the entire accretionary wedge experiences large strains (Figure 1b, Figure  S1 in Supporting Information S1), resulting in repeated strain localization on the same splay fault geometries.
During the slip event, the amount of stress drop on the different splay faults is highly variable in the SC model ( Figure 4, Figure S11-S14 in Supporting Information S1). The largest splay fault SF6 shows stress drops up to 7.5 MPa in the basalt and sediments and an isolated large stress drop of 14.2 MPa in the accretionary wedge sediments (Figure 4c). SF5 generally exhibits stress drops of 1-2 MPa, with the deepest part of the fault featuring stress drops up to 3.3 MPa ( Figure S14 in Supporting Information S1). SF4 shows stress drops of 1-2 MPa near the branching point. The shallow, small splay faults (SF1-3) do not experience any significant stress drop.

Dynamic Earthquake Rupture
We compare a model in which only the megathrust is allowed to rupture (Figures 5a and 5c; Van Zelst et al., 2019) to the model in which the megathrust and the six splay faults are theoretically allowed to slip. The ruptures show similar rupture speeds, but different rupture duration with the model including splay faults rupturing for longer (89 s instead of 82 s). Slip differs significantly between the two ruptures with the model including splay faults exhibiting lower slip and slip velocities ( Figure 5).
Earthquake rupture initiation is non-prescribed and solely driven by the initial conditions from the geodynamic seismic cycle model. After a two-stage nucleation at very low slip rates (a 4 s period of low rupture speed, followed by a 2 s high speed phase), spontaneous rupture emerges on the megathrust ((1) in Figure 5b). Subsequently, rupture propagates both updip and downdip, where it is spontaneously arrested at the brittle-ductile transition (2) in both models (Figures 5a and 5b). In the updip direction, the main rupture front in the splay fault model encounters SF6 after 14.1 s. While the dynamic activation of SF6 appears to resemble rupture branching (DeDontney et al., 2011; Movie S1 and S2), we observe a high degree of complexity on smaller scales. The passing megathrust rupture dynamically unclamps SF6, that is, there is a decrease in the normal stress σ n (Oglesby et al., 2008), which results in negligible slip over 1 km of the splay fault close to the fault junction without spontaneously propagating rupture. Subsequently the rupture jumps from the megathrust to SF6 due to dynamic triggering, omitting the deepest 3 km of the splay fault that had a higher initial strength excess (Figure 3), which only ruptures in a down-dip direction after 18 s ( (3) in Figure 5j). Unilateral dynamic rupture then propagates updip on the splay fault with slip velocities of 4.7 m/s. Simultaneously, ahead of this rupture front, secondary ruptures are dynamically triggered by the main megathrust rupture (4) leading to an apparently very high updip splay rupture speed. Behind this first, apparently fast splay rupture front, we observe fault reactivation due to multiple passing rupture fronts on the megathrust and free surface reflected seismic waves (5), resulting in a static slip maximum of 13.8 m. Due to the splay fault rupture, the slip velocities on the megathrust updip of the splay fault are sharply reduced compared to a model which only ruptures the megathrust. This leads to a slip discontinuity on the megathrust (Figure 5d).
The main rupture front on the megathrust passes SF5 without activating it (6), that is, neither by branching nor dynamic triggering (Figure 5i). This difference in splay fault activation dynamics can be attributed to the local non-optimal orientation of SF5 near the branching junction, which forms an effective barrier for dynamic rupture propagation. Instead, SF5 is activated at shallow depths of ∼5 km at 32.8 s due to waves reflecting from the free surface (7). Multiple rupture fronts then propagate downdip (i.e., hosting reverse slip) on SF5, but the unfavorably oriented deepest 2.5 km of SF5 never fully ruptures (8). Since the passing of the primary megathrust rupture front does not trigger slip on SF5, there is no decrease in slip rate on the megathrust after it passes SF5.
Although the passing of the main rupture front induces small slip rates on SF1-4 on the order of ∼0.02 m/s due to unclamping, they only rupture afterward in a self-sustained manner at slip rates larger than 1 m/s due to static and dynamic stress changes. These are induced by secondary rupture front complexity on the megathrust as well as on SF5 and SF6 and multiple reflected (trapped) waves within the sedimentary wedge. The long rupture duration on these shallow splay faults leads to a maximum slip of 12.6 m for SF4 and 10.0 m, 8.1 m, and 8.0 m for SF1-3, respectively, barring some numerical outliers. Since slip occurs on the splay faults and the slip velocity on the megathrust is reduced when the rupture interacts with a splay fault, the maximum slip on the megathrust in the model including splay fault rupture (48.9 m) is lower than in the model without splay fault rupture (57.6 m). The decrease in slip on the megathrust at the branching points with the splay faults is visible as sharp discontinuities in Figure 5d, which contribute to the slip observed on the splay faults (Figures 5k-5p). This indicates that slip is transferred from the megathrust onto the splay faults for SF1-4 and 6.
In summary, we find that all splay faults rupture coseismically, albeit in three different fashions. SF6 is unclamped by the primary rupture front, SF1-4 are dynamically triggered by the reactivated megathrust, and SF5 slips near the surface due to dynamic stress transfer from wave reflections from the free surface. We ascribe these differences in fault activation to variable along-fault strength excess and fault orientation with respect to the prevailing stress field. However, we generally observe that the splay faults are favorably orientated with respect to the stress field and have low strength excess resultant from high pore-fluid pressures.
The maximum stress drop, computed on-fault, on the megathrust on the order of ∼17 MPa is comparable in the models with and without splay faults ( Figure S11, S15a and S15b in Supporting Information S1). SF 6 shows the largest stress drop of all splay faults on the order of ∼19 MPa (Figure 4c and Figure S15 in Supporting Information S1). The other splay faults show maximum stress drops of 2.5-6.5 MPa, with the deeper splay faults exhibiting larger stress drops than the shallow splay faults (Figure 4, Figures S12-S15 in Supporting Information S1). In general, the stress drop is relatively constant along the splay faults, with the exception of the branching point which typically shows a larger stress drop than the rest of the splay fault. Splay fault 6 is the only splay fault which shows varying stress drop along the fault with higher stress drops in the basalt and incoming sediments directly below the basalt.
The model without splay faults has relatively uniform static vertical surface displacements of ∼5 m and a smooth profile of horizontal displacements of 47.8 m seawards ( Figure 6). In contrast, the model with splay faults shows clear vertical surface displacement peaks corresponding to the shallow tips of the splay faults near the surface ( Figures 6b and 6c). The wavelengths of these peaks are ∼80%-95% smaller than the wavelengths of the vertical surface displacements due to rupture purely on the megathrust. The largest peak of 9.3 m at 180 s is associated with SF6, whereas the other peaks with amplitudes ranging from 4.7 to 6.5 m are associated with SF1-5. Hence, rupture on splay faults increases the amplitude of the vertical displacements up to 86%. The amounts of vertical displacement and slip are not linearly correlated ( Figure S16 in Supporting Information S1) as other factors, such as the dip angle and slip distribution on the fault also play a role. The effect of splay fault rupture is less pronounced in the horizontal displacements with a 17% lower amplitude of the horizontal displacements compared to the model without splay faults (Figure 6d).

Tsunami Propagation and Inundation
The tsunami resulting from the model without splay faults consists of a single wave with a wavelength of 300 km and a maximum sea surface height of 6.5 m (Figure 7a). It arrives at the beach after 11 min and it takes a total of 74.5 min for the whole wave to arrive at the coast. There is one episode of flooding at the coast with a run-up distance of 1,250 m. Here, we define run-up distance as a measure of how far inland the tsunami reaches horizontally compared to the original coastline (Satake, 2015).
In the model including six splay fault ruptures, the tsunami consists of one high wave crest corresponding to slip on SF6 ((7) in Figure 7b) and a broad wave packet resulting from slip on the other splay faults and shallow part of the megathrust ((1-6) in Figure 7b). Similar to the tsunami of the model without splay faults, the waves span a region of 300 km, but have smaller individual wavelengths. The tsunami first reaches the coast after 11 min and impacts the coast until 71.3 min. It reaches a maximum sea surface height of 12.2 m, which is almost double the height of the model without splay faults. Besides that, the flooding at the coast occurs in two episodes (Figure 8) in contrast to one flooding episode for the model without splay faults. The first episode is related to the large wave resultant from rupture on SF6, whereas the second episode relates to a wave originating from the interference of the smaller waves related to the other splay faults and shallow megathrust. The run-up distance of the tsunami is 2,210 m, which is 77% larger than that of the tsunami sourced by a rupture without splay faults.

Discussion
Observational studies of accretionary wedges image multiple splay faults which pose a tsunami hazard (Kopp, 2013). It is difficult to asses if multiple splay faults rupture during a single earthquake and how that affects the ensuing tsunami. It is often assumed that only one splay fault at the time is seismically active in conjunction with the megathrust (e.g., DeDontney & Hubbard, 2012; Park et al., 2002). However, the uncertainty in tsunami source location (Sibuet et al., 2007;Waldhauser et al., 2012) and the exact locations of the ruptured fault planes could allow for multiple, closely spaced (partially) ruptured splay faults during a single earthquake. Numerical models can shed light on the process of rupture on multiple splay faults, but initial fault stresses are difficult to constrain (e.g., Van Zelst et al., 2019) and the choice of numerical discretization method can hamper the geometric complexity of dynamic rupture models (e.g., DeDontney & Hubbard, 2012). Here, we explicitly account for self-consistent initial fault stresses, complex topo-bathymetry, and a shallowly dipping megathrust intersecting with six different splay fault geometries. This linked framework including geodynamics, seismic cycles, dynamic rupture and tsunami propagation and inundation allows us to address questions about the viability and likelihood of splay fault ruptures and their impact on tsunamigenesis.
It is currently unknown under what circumstances earthquakes will produce large offsets of the seafloor, which is one cause of unexpectedly large tsunamis (e.g., Brodsky et al., 2021;Dunham et al., 2020). Slip on the megathrust propagating onto splays through dynamic or static stress changes has been inferred for past and recent tsunamigenic earthquakes (e.g., Cummins & Kaneda, 2000;Fan et al., 2017;Obana et al., 2017). Our results highlight that studying compound rupture of megathrusts and multiple or segmented splay faults is important for the assessment of future hazardous events and to better understand the details of near-trench rupture processes that control seafloor uplift and hence tsunami generation (Madden et al., 2020;Saito et al., 2019;Satake, 2015;Tanioka & Satake, 1996;Ulrich et al., 2022;Wirp et al., 2021). Future efforts could aim to include region-specific observations, such as high-resolution seismic imaging and geological data in modeling workflows that link earthquake source models to tsunami models to improve our understanding of tsunamigenesis.

Fault Geometries
One of the choices in our coupled modeling framework is the choice of fault geometries in the SC model as input for the DR model. The chosen fault geometries determine which stresses and strengths are ultimately used as input for the DR model, where the initial stresses on the faults are crucial for the ensuing dynamic rupture.
The chosen megathrust geometry for this slip event is picked from the highest visco-plastic strain rate during the event . For all slip events in the SC model, the megathrust is blind and follows the lithology contrast between the basaltic oceanic crust and the incoming sediments below the sedimentary wedge. The megathrust is consistently located at that location, because it is the location of the largest differential strain build-up and thus largest interseismic stressing rates. The incoming sediments do not completely subduct together with the slab, as parts are also accreted to the accretionary wedge (Clift & Vannucchi, 2004;Cloos & Shreve, 1988;Von Huene & Scholl, 1991). A blind, or buried, megathrust is thought to be less common in nature, but has been inferred for, e.g., the Cascadia subduction zone where no evidence of the megathrust breaching the seafloor has been found (e.g., Flueh et al., 1998;Lotto et al., 2019). Coupled earthquake-tsunami models by Lotto et al. (2019) show that the tsunami profile resulting from a buried megathrust rupture with simple loading and strength properties is complex. Many small peaks and troughs are caused by the effect of enhanced shallow slip and the vertical seafloor displacement, which we also observe in our model ( Figure 6).
Similar to the megathrust geometry, all six splay fault geometries considered here are blind with the tip of the splay faults located at 2 km depth on average. This results in more gradual, and hence less discontinuous surface displacements compared to studies where splay faults breach the seafloor (e.g., Li et al., 2014;Ulrich et al., 2022). In addition, the surface displacements resulting from rupture on blind faults could also have different, and specifically smaller, amplitudes that might affect tsunami height. However, the surface displacements associated with interactions between the rupture and the free surface typically have lower wavelengths (Nielsen, 1998) that are not thought to have a strong effect on the tsunami (Saito et al., 2019). Therefore, we hypothesize that the use of blind faults in this work does not affect our main conclusions.
The splay fault geometries do not significantly change between slip events in the SC model and we observe no stress drop on SF1-3 during events in the SC model (Figure 4, Figures S12-S14 in Supporting Information S1), although strain localizes and slip occurs on them. We do observe a stress drop in the SC model on the larger splay faults of up to 7 MPa. However, in the DR models, we observe significantly more stress drop on each of the splay faults ( Figure 4, Figures S12-S14 in Supporting Information S1). This indicates the importance of the coupled modeling framework, where the DR model fully resolves the ruptures, resulting in stress drops on the ruptured faults and incorporating dynamic waves effects. These latter effects in particular have been shown to be important for SC models (Thomas et al., 2014;Van Zelst et al., 2019). We speculate that the repeated reactivation of the same splay fault geometries in the SC model over each seismic cycle might not occur if two-way coupling of the codes were to be employed such that the resulting stress state after the dynamic rupture would be fed back into the SC model. The reason for this is that the stresses on the splay faults would be reduced further after rupture due to larger dynamic stress drops, such that subsequent interseismic periods start with lower stresses that then would need to be increased slowly and steadily before the faults rupture again. Since the build-up of stress mostly occurs near the downdip limit of the seismogenic zone below which ductile creep leads to differential displacements, stresses would first need to be transferred updip to the splay fault locations (see e.g., Dal Zilio et al., Herrendörfer et al., 2015;Kammer et al., 2015). Reloading the splay faults for activation may thus take a significant amount of time.

Rupture on Splay Faults
Our models show that all six splay faults rupture when we use the self-consistent initial conditions from the SC model. This is partially due to the orientation of the splay faults with respect to the local stress field, which is generally favorable for rupture according to Coulomb theory (Figure 2b; e.g., Kaus, 2010;Wang & Hu, 2006). The splay faults also exhibit low strength excess, particularly at shallow depths (Figure 3, Figures S8-S10 in Supporting Information S1), indicating that they are close to failure at the start of the rupture (e.g., Li et al., 2014). Here, we define strength excess as dr yield − , where dr yield is the fault yield stress and τ is the initial shear stress. The low strength excess of the shallow splay faults can largely be explained by the low strength of sediments in the sedimentary wedge due to the presence of fluids and prevalent high pore-fluid pressures (van Dinther et al., 2014).
Here, we assume a high pore-fluid pressure ratio of 0.95. This results in reasonable recurrence intervals on seismic cycle time scales, while allowing for subduction along a shallow megathrust on geodynamic time scales (van Dinther, . In addition, 3-D dynamic rupture simulations  support the presence of high coseismic pore fluid pressure at megathrusts (Audet et al., 2009;Saffer & Tobin, 2011;Tobin & Saffer, 2009). The deeper splay faults SF4-6 are not as close to failure as the shallower splay faults, but still rupture due to the overall energetic rupture and wave reflections and the resulting stress changes. Splay fault 5 in particular does not rupture at the branching point due to the large strength excess and a high branching angle (21.8°) that misaligns SF5 with respect to the local stress field. Instead, it is activated at shallow depths due to reflecting waves from the free surface where the strength excess on the fault is small. Hence, our results suggest that multiple splay faults rupture during an energetic event with reflecting waves when (a) they are favorably orientated with respect to the local stress field for rupture, that is, they are strong faults according to Andersonian faulting theory, and (b) they have a low strength excess, that is, they are close to failure.
Slip on our simulated splay faults is on the order of ∼10 m. In nature, slip on splay faults is generally hard to observe, but Chapman et al. (2014) report an estimate of 3 m slip on a splay fault during the 1964 Alaska earthquake. Similarly, Cummins and Kaneda (2000) infer up to 3.5 m splay fault slip during the 1946 Nankai earthquake. In recent 3D dynamic rupture models of the 2004 Sumatra-Andaman earthquake, three co-seismically activated large-scale splay faults host larger slip on the order of 6-8 m . Additionally, another, shorter splay fault near the trench with a steeper dip hosts up to tens of meters of slip in these models. While observational estimates of splay fault slip are lower than our simulated slip, the modeled vertical surface displacements do align with those typically associated with splay faults. For example, Suito and Freymueller (2009) and Suleimani and Freymueller (2020) report 10 m of vertical surface uplift for the splay fault rupture of the 1964 Alaska earthquake.
The relatively high geodynamically constrained 2D dynamic rupture fault slip can be reduced  by simulating rupture in three dimensions, as illustrated in Li et al. (2021) for a strike-slip setting and Madden et al. (2020) and Wirp et al. (2021) for a subduction zone setting, where the stresses and strengths are adapted to avoid unilateral nucleation along the 2D fault and thereby reduce fault slip. The 3D megathrust dynamic rupture models in Wirp et al. (2021) additionally use a constant, not geodynamically informed, characteristic slip distance, which further reduces the amount of slip. We expect that considering off-fault plasticity would further decrease shallow fault slip on the megathrust and splay faults . Future studies could use these findings to expand on the approach presented in this article to obtain more realistic slip values on splay faults. Since our models are restricted to two dimensions and show a relatively large amount of slip on the splay faults, we caution that our models may overestimate the absolute tsunami heights and run-up distance.

Tsunamis Resulting From Rupture on Splay Faults
In the tsunami models, the effect of slip on splay faults is visible in the propagating wave and the inundation pattern at the coast (Figures 7 and 8; Goda et al., 2014). The tsunami model without splay fault rupture also shows localized crests (Figure 7a), although to a lesser extent. This indicates that crests in the tsunami data cannot exclusively be contributed to splay fault rupture. Similarly, the absence of complexity in the tsunami data, particularly with regards to the second wave packet, does not necessarily mean that rupture only occurred on one splay fault. Indeed, the effect of rupture on other, smaller splay faults might not be distinguishable based on tsunami data alone. To relate our findings directly to tsunami data, the here found splay fault effects should be analyzed with more complex bathymetry and 3-D complexity in future studies (Bletery et al., 2015;Matsuyama et al., 1999;Tonini et al., 2020;Ulrich et al., 2019). Recent studies using a similar methodology to the one presented here have already attempted this for megathrust-only events (Madden et al., 2020;Wirp et al., 2021). However, one of the major limitations in these 3-D studies is the uncertainties in how to accurately account for any lateral variation in the initial stresses and strengths on the megathrust since the considered geodynamic seismic cycle model is two-dimensional. This limitation is enhanced when complex splay fault geometries are considered in addition to the megathrust. Lastly, the here used hydrostatic shallow-water-based tsunami modeling approach does not fully capture smaller-scale complexity during tsunamigenesis nor dispersive effects during tsunami wave propagation, and future studies may extend our approach to account for a nonlinear hydrodynamic response (Kim et al., 2017;Saito et al., 2019), corrections for dispersive Earth elasticity, and non-dispersive water compressibility (Tsai et al., 2013) or fully coupled seismic, acoustic, and gravity modeling Lotto et al., 2019).

Conclusions
We develop and use a novel modeling framework that combines geodynamics, seismic cycles, dynamic rupture, and tsunami propagation, and inundation to understand the rupture dynamics and tsunamigenesis of multiple splay faults. This linked framework constrains the geometry, stress, and strength of the megathrust, six splay faults and the surrounding rocks in a physically self-consistent manner. To first order the on-fault stresses of splay faults increase linearly with depth. However, deviations occur due to variations in lithology, pore fluid pressure and deep stress release. In our geodynamic seismic cycle model, we analyze theoretical fault growth angles with respect to the principal stress direction assuming a Coulomb angle in a compressional stress regime. We find that large portions of most splay faults are favorably orientated, aiding activation during megathrust earthquake rupture. In addition, the splay faults generally have low strength excess due to high pore fluid pressures in large parts of the sedimentary wedge, indicating that they are close to failure.
We find that all splay faults are dynamically activated by various mechanisms in the dynamic earthquake rupture model, such as the passing of the megathrust rupture front and stress changes from reflected waves in the sedimentary wedge. We observe rupture branching from the megathrust to the largest splay fault, and detail the small-scale dynamic fault interactions of unclamping and rupture jumping. The main rupture front on the megathrust passes all other splay faults without activating them by branching. We attribute this difference in splay fault activation dynamics to local variations in strength excess and non-optimal orientations of all shorter splays near the branching junction, which form an effective barrier for dynamic rupture propagation. The second largest splay, SF5, is slipping only partially and in a down-dip reverse manner due to waves reflecting from the free surface. While the passing of the main rupture front unclamps the four shorter splays SF1-4, they rupture delayed due to static and dynamic stress changes from megathrust rupture complexity and slip on the respectively larger splays. Eventually, all splay faults experience slip reactivation during the same earthquake simulation due to stress changes induced by multiple reflected (trapped) waves within the sedimentary wedge.
Rupture on the largest splay fault results in a local, short-wavelength increase in tsunami height. A second, broad wave packet in the tsunami is due to slip on multiple splay faults and the shallow megathrust. This wave packet is similar in width to the one produced in the model with a pure megathrust rupture, albeit with larger amplitude and shorter wavelength wave crests relating to the activation of each splay fault. However, at the coast, the multiple wave crests can no longer be distinguished, making it difficult to determine if multiple splay faults ruptured from tsunami data alone.
Our multi-physics models imply that simultaneous rupture on multiple splay faults is mechanically viable and is facilitated by the low strength and favorable stress orientation of the faults resulting from long-term tectonics and the strong dynamic (re-)activation potential of splay faults. It is therefore important to take the possibility of rupture on multiple splay faults into consideration in tsunami hazard.

Appendix A: Defining Splay Fault Geometries From the Geodynamic Seismic Cycle Model
To provide well-defined fault geometries as input to the DR model, we approximate the splay geometries in the sedimentary wedge by analyzing the visco-plastic strain ɛ vp visualized as the accumulated visco-plastic slip d = 2Δx ⋅ ɛ vp in Figure 1b with Δx = 500 m representing the fault width (van Dinther, . We calculate ɛ vp from the second invariant of the visco-plastic strain rate = √ 2 + 2 according to where t = 0 is the coupling time step for which the output of the SC model is used as input for the DR model according to Van Zelst et al. (2019). Δt is the time step (5 years) of the SC model, and t max is the final time step of the coupled SC event. We verified that the viscous component in the visco-plastic strain rate is negligible , such that shows the effect of plastic rock behavior.
To pick discrete splay fault geometries from the visco-plastic strain distribution in the SC model ( Figure S1 in Supporting Information S1), we only consider regions where the minimum slip is 0.16 m (Figure 1b). This corresponds to a strain rate of 10 −12 s −1 . This threshold highlights the regions of strain during the event, and hence the splay fault geometries. We then pick six representative splay fault geometries (Figure 1b). We show the complete procedure for picking each splay fault in Figure A1 for splay fault 6 (see Figures S2-S6 in Supporting Information S1 for the other splay faults). For each splay fault, we manually determine the x-extent of the fault.
x min is initially determined by visual inspection of Figure 1b, which is then iteratively adjusted based on the highest strain. We choose an arbitrarily large value of x max , which is later adjusted based on meshing requirements at the branching point between the splay fault and the megathrust. Then, for each nodal x-coordinate, we pick the z-coordinate with the highest strain in the sedimentary wedge, that is, disregarding the megathrust at which the largest strain is accumulated (red dots in Figure A1a). We manually reposition any outliers that clearly belong to adjacent faults to align with the observed strain localization (red dots with cyan borders in panel (a) of Figure A1; see Figures S2-S6 in Supporting Information S1).
We then smooth the fault geometry with a moving average low-pass filter scheme with a span of 25 points (red dots in Figure A1b; Van Zelst et al., 2019). To ensure that the splay faults connect to the megathrust in the most efficient manner for the mesh, we limit the x-extent of the splay faults (red dots with yellow borders in Figure A1b). The final geometries of the splay faults are then shown in Figure A1c and Figures S2-S6 in Supporting Information S1. Details of the splay fault geometries are listed in Table S1 in Supporting Information S1 and the full geometry of the splay faults can be found in Data Sets S1 to S6.
We do not connect the splay faults to the surface, because there is no indication that they reach the surface in the geodynamic seismic cycle model ( Figure A1). This is due to the predefined decreased pore-fluid pressure ratio of 0.4 in the top kilometer of the SC model. Hence, we only consider blind splay faults here. There are also fault geometries other than splay faults present in the yielding sedimentary wedge of the geodynamic seismic cycle model, such as antithetic fault planes (Figures 1b, Figure A1). However, here we focus solely on the more conventional splay fault geometries and do not include any antithetic fault geometries to limit the complexity of our model.