The State of Pore Fluid Pressure and 3-D Megathrust Earthquake Dynamics

We study the effects of pore fluid pressure ( P f ) on the pre-earthquake, near-fault stress state, and 3-D earthquake rupture dynamics through six scenarios utilizing a structural model based on the 2004 M w 9.1 Sumatra-Andaman earthquake. As pre-earthquake P f magnitude increases, effective normal stress and fault shear strength decrease. As a result, magnitude, slip, peak slip rate, stress drop, and rupture velocity of the scenario earthquakes decrease. Comparison of results with observations of the 2004 earthquake support that pre-earthquake P f averages near 97% of lithostatic pressure, leading to pre-earthquake average shear and effective normal tractions of 4–5 and 22 MPa. The megathrust in these scenarios is weak, in terms of low mean shear traction at static failure and low dynamic

In seismogenic regions of subduction zones, lower P f conditions have been proposed as a mechanism for locking (Saffer & Tobin, 2011). Heise et al. (2017) co-locate a geodetically identified locked region with a patch of high electrical resistivity attributed to lack of fluid or low P f on the Hikurangi subduction interface, while shallow creep occurs in a region of conductivity that can be explained by high fluid production or high P f (Heise et al., 2013). However, heat flow studies (Gao & Wang, 2014) and force-balance inversions (Lamb, 2006) find shear to normal stress ratios that indicate high P f near the megathrust. Lamb (2006) finds evidence for P f at 95% of the lithostatic pressure at 7 of 9 subduction zones, including Sumatra. Two exceptions to this are Northern Chile and Tonga, with P f at 81% of the lithostatic pressure.
Temporal variation in P f is central to the fault-valve model of Sibson (1992Sibson ( , 1994, which attributes earthquakes to both tectonic loading (shear stress building until an earthquake occurs) and fluid-pressure cycling (P f building and effective normal stress falling over time until an earthquake occurs). Petrini et al. (2020) show that fluid pressure variations in time can control subduction zone seismic cycling. Analyses of borehole fluids suggest cycles of 10,000-100,000 yr (Saffer & Tobin, 2011), which may correlate with fault formation, while shorter period variations correlate with slow slip events in Costa Rica. In addition, observed increases in V p /V s following the 1995 M 8 Antofagasta earthquake (Husen & Kissling, 2002) suggest the rapid movement of fluid during or directly after megathrust earthquakes. Eberhart-Phillips et al. (1989) note that such changes can occur only when P f is near-lithostatic.
This variety of observations and inferences about P f in subduction zones is reflected in the variety of ways that P f is considered in faulting and earthquake models. Quasistatic models of fault slip may not incorporate P f explicitly, but set stress gradients that produce reasonable fault slip distributions (e.g., Madden et al., 2013;Madden & Pollard, 2012). Models of earthquake sequences and rupture dynamics commonly prescribe normal stress following effective stress theory as σ n − P f , where σ n is the normal stress (Brace & Kohlstedt, 1980;Hubbert & Rubey, 1959). P f typically increases with depth and is chosen ad-hoc to help reconcile realistic earthquake characteristics with friction and fault shear strength (Kozdon & Dunham, 2013;Liu & Rice, 2005;Ulrich, Gabriel, et al., 2019;Wollherr et al., 2019). Ulrich et al. (2022) and Uphoff et al. (2017) incorporate near-lithostatic P f following depth-dependent gradients into large-scale, three-dimensional (3-D) dynamic rupture models. Others initialize dynamic rupture models with conditions, including initial P f , from geodynamic and seismic cycling modeling that capture long-term subduction zone deformation and fluid flow Wirp et al., 2021;van Zelst et al., 2019). Rice (1992) shows that fluid at elevated pressures within a fault zone may follow the same gradient with depth as the lithostatic stress, causing constant effective normal stress with depth. Data from crustal sedimentary rocks support this theory (Suppe, 2014). This condition is assumed in some dynamic rupture models (e.g., Ramos & Huang, 2019;Ramos et al., 2021), but not others (e.g., Kozdon & Dunham, 2013;Lotto et al., 2019;Ulrich et al., 2022;Ulrich, Vater, et al., 2019). Other models consider the coupled, dynamic effects of fluids, such as dilatancy (e.g., Aochi et al., 2014;Segall & Rice, 1995) and thermal pressurization (e.g., Garagash, 2012;Rice, 2006;Schmitt et al., 2011;Segall & Bradley, 2012). Recent two-dimensional (2-D) antiplane earthquake sequence modeling by Zhu et al. (2020) couples earthquake and pore fluid dynamics by incorporating fluid migration and periodic P f variations over earthquake cycles. These models produce fluid-driven aseismic slip at the base of the seismic zone, large earthquakes, and earthquake swarms. 2-D seismo-hydro-mechanical modeling of subduction zone earthquake cycling shows high P f moving progressively updip due to compaction inside an evolving fault, eventually leading to a seismic event (Petrini et al., 2020). P f prior to an earthquake can be constrained by these observations and inferences with simultaneous consideration of the normal stress and static frictional strength of a megathrust, but it has not been measured directly and little data is available, particularly deep along subduction zones. Few studies integrate knowledge about megathrust mechanics with megathrust earthquake rupture dynamics to study P f at the time of rupture. Specifically, 3-D dynamic simulations at the megathrust scale that take realistic slab geometries into account remain challenging. To supplement this gap, we explore the dynamic effects of different hypotheses about P f magnitude and gradient in megathrust systems using a 3-D dynamic earthquake rupture and seismic wave propagation model that matches near-and far-field seismic, geodetic, geologic, and tsunami observations of the 2004 Sumatra-Andaman earthquake and Indian Ocean tsunami (Ulrich et al., 2022;Uphoff et al., 2017).
Our focus is to highlight the effects of pre-earthquake P f conditions on earthquake behavior within a structurally complex megathrust system. We analyze how various hypotheses on P f magnitude and depth gradient affect the pre-earthquake stress state near a megathrust, the subsequent earthquake rupture characteristics, and the postseismic stress field. Specifically, we generate six scenario earthquakes with P f magnitudes at 31%, 62%, 93%, and 97% of the lithostatic pressure and under two different depth gradients that cause either increasing or constant normal stress near the megathrust. We compare results against observations of the 2004 earthquake as well as general observational inferences about subduction zone earthquakes.
We note that the range of pre-earthquake conditions captured by our six scenarios may reflect the variety of conditions present along a single megathrust at the same time, due to spatial variations in P f magnitude and/or gradient. In addition, hydromechanical processes likely vary in space and time as a consequence of rock deformation processes that modulate the permeability of both fault and host rocks, in turn affecting fluid diffusion. Coupling these processes during the full seismic cycle to determine realistic fluid conditions at the start of earthquake rupture is a clear future step. However, modeling these processes in 3-D is beyond the state of the art, despite the recent progress of 2-D numerical models reviewed above. Our results provide key advances regarding the influence of P f on earthquake behavior and provide opportunity for future calibration with site-specific friction and pore-fluid measurements to constrain dynamically plausible megathrust strength and P f gradients.

Computational Model
The earthquake models are performed with SeisSol (www.seissol.org), a software package that solves for dynamic fault rupture and seismic wave propagation with high-order accuracy in space and time. SeisSol solves the seismic wave equation in velocity-stress formulation using an Arbitrary high-order DERivate Discontinuous Galerkin (ADER-DG) scheme (Dumbser & Käser, 2006). Computational optimizations target supercomputers with many-core CPUs Heinecke et al., 2014;Krenz et al., 2021;Rettenberger et al., 2016). SeisSol uses local time stepping, which increases runtime efficiency by decreasing dependence of the time-step on the element with the smallest radius (Breuer et al., 2016;Uphoff et al., 2017;Wolf et al., 2020). Following the SCEC/USGS Dynamic Rupture Code Verification exercises (Harris et al., 2018(Harris et al., , 2009), SeisSol has been validated against several community benchmarks (De La Puente et al., 2009;Pelties et al., 2012Pelties et al., , 2014Wollherr et al., 2018).

Structural Model
The structural model and computational mesh are shown in Figure 1. Use of an unstructured tetrahedral mesh allows for a realistic representation of the non-planar slab interface, splay faults, curved oceanic crust, and high-resolution bathymetry. The megathrust geometry follows Slab1.0 (Hayes et al., 2012). The splay faults, one longer back thrust, and two shorter forethrusts, are interpreted from aftershocks (Waldhauser et al., 2012), seafloor observations (Chauhan et al., 2009;Sibuet et al., 2007;Singh et al., 2008), and tsunami modeling (DeDontney & Rice, 2011). The mesh for this model has elements with edge lengths of 1 km along the faults, 4 km at the surface, and 100 km in the volume far from the fault; mesh resolution varies gradually between these conditions. We ensure that the element size along the megathrust and splay faults is sufficient to capture the cohesive zone following the analysis in Wollherr et al. (2018) and detailed in Appendix A.
The regional rock properties are adapted from Laske et al. (2013) and include four layers of oceanic crust and four layers of continental crust with the properties outlined in Table 1. As shown in Figure 1, the layers of oceanic

Fluid Pressure, the Regional Stress Field and Fault Tractions
We assume a laterally homogeneous regional stress tensor. Its orientation is from an inversion of focal mechanisms near the hypocenter of the 2004 Sumatra-Andaman earthquake by Karagianni et al. (2015;region 7.1.22). Taking a compression negative sign convention, the maximum compressive stress (σ 3 ) has an azimuth of 225° and plunges 7°. The intermediate principal stress (σ 2 ) has an azimuth of 315° and plunges 7°. The least compressive stress (σ 1 ) has an azimuth of 90° and plunges 80°. In all scenarios, the absolute stresses are proportional to the lithostatic stress (σ v = ρgz, where ρ is the density of rock, g is gravitational acceleration and z is depth) as σ 1 = 0.98σ v , σ 2 = 1.5σ v , and σ 3 = 2σ v . Below 23 km depth, we taper the differential stress to zero at 50 km depth to approximate the transition from brittle to ductile deformation.
We present six scenarios with different P f magnitudes and depth gradients applied to this absolute stress state (Table 2). Following the effective stress principle (Brace & Kohlstedt, 1980;Hubbert & Rubey, 1959), the effective principal stresses ( for each scenario are determined relative to  However, Rice (1992) shows that fluid at elevated pressures within a fault zone may follow the same gradient with depth as σ v , which causes a constant effective normal stress with depth. We follow this assumption in Scenarios 5 and 6, where high and very high P f follow the gradient in σ v , but are offset by constant values (K) of 42 MPa in Scenario 5 and 20 MPa in Scenario 6: We refer to this as a lithostatic P f gradient and it is applied below 5 km depth.
To resemble borehole stress and fluid-pressure measurements in continental margins (e.g., Suppe, 2014), we apply a lithostatic gradient above 5 km in both scenarios. On average over the rupture area, P f in Scenarios 5 and 6 is 93% and 97% of σ v , respectively, mirroring values in Scenarios 3 and 4. The lithostatic P f gradient, the absolute principal stresses and the effective principal stresses are shown for Scenario 6 in Figures 2d-2f.
In all scenarios, stresses and P f vary only with depth and do not vary with horizontal location. As P f increases in these scenarios, the magnitudes of ′ , ′ 3 , ′ 2 and ′ 1 all decrease. In addition, effective normal stress magnitudes (through the effective mean stress) and shear stress magnitudes (through the effective deviatoric stress) decrease as well. Figure 3a shows the relatively low stress magnitudes present at all orientations when a very high P f magnitude is applied in Scenario 4, while also demonstrating how these stress magnitudes increase with depth in Scenarios 1-4. Figure 3b shows the relatively low stress magnitudes present at all orientations when a very high P f magnitude is applied in Scenario 6, while also demonstrating how these stress magnitudes are constant with depth in Scenarios 5 and 6.
The initial shear and effective normal tractions, τ s and ′ , are determined by projecting the local effective stress tensors onto the non-planar megathrust and splay faults. As for the shear and effective normal stresses, both τ s and ′ magnitudes decrease overall as P f increases from scenario to scenario. In Scenarios 1-4, τ s and ′ magnitudes increase with depth, while in Scenarios 5 and 6, both are relatively constant with depth. The pre-earthquake tractions are shown for each scenario in Figure 4 and mean values are summarized in Table 2. Setting the effective stress magnitudes relative to ′ as we do maintains the same ∕ ′ distribution on the megathrust across all scenarios ( Figure B1), which isolates the influence of P f on earthquake behavior, as desired in this study.
While the on-fault tractions mirror the near-fault stresses in many ways, our 3-D, geometrically complex fault structure comprised of a non-planar megathrust and splay faults modulates the fault traction distributions. As a result, they depart in certain locations from the linear stress gradients and feature additional spatial variations and heterogeneity, as both τ s and ′ vary with fault geometry in all scenarios. Figure B1 illustrates how this distribution varies due to the non-planar megathrust geometry. In Scenarios 5 and 6, where the P f gradient is lithostatic and τ s and ′ are relatively constant with depth, the variation due to the megathrust geometry is ≈5 MPa.

Failure and Spontaneous Propagation
In all scenarios, dynamic earthquake rupture starts by forced nucleation in the southeastern corner of the megathrust at 30 km depth. Failure occurs when τ s exceeds the static fault strength, T fs , which is determined from the on-fault frictional cohesion, c, and the product of the coefficient of static friction, μ s , and ′ as (compression is negative):

Table 2
Initial Conditions for All Scenarios c is the frictional strength of a fault in the absence of ′ . In situ, c depends on local mineralogy and lithology, but here c is used as a standard proxy for near-surface behavior that we do not model explicitly, mainly the constitutive behavior of shallow sediments in the near-trench region (e.g., Harris et al., 2018;Kaneko et al., 2008). We set c = 0.4 MPa along the megathrust and splay faults below 10 km depth, but c linearly increases to 15 MPa from 10 to 0 km depth. Due to topography, the intersection of the fault and the seafloor ranges between 3 and 5 km depth, so maximum c values on the faults at the seafloor range from 8 to 11 MPa. For further discussion of c, please see Section 5.1 and Appendix B. . Whether the P f gradient is sublithostatic or lithostatic changes the effective principal stress depth profiles; magnitudes scale inversely with P f magnitude (c and f). Stresses and P f vary only with depth, not with horizontal location.
We assign μ s = 0.4 to all faults in all scenarios. Borehole estimates of stress in upper crustal rocks suggest that rocks follow Byerlee's law with μ s = 0.6-1.0 (Suppe, 2014;Townend & Zoback, 2000. Our choice of μ s = 0.4 is motivated by the lithology of the shallow megathrust characterized by high, clay-rich sediment input that is progressively strengthened by dehydration and compaction near the megathrust (Hüpers et al., 2017). Our choice to keep μ s constant across all faults and all scenarios allows us to here focus on the effects of P f magnitude and depth gradient.
We apply a linear slip-weakening friction law (e.g., Andrews, 1976) to represent dynamic weakening of a fault after failure. μ s decreases to the coefficient of dynamic friction, μ d , over the slip-weakening distance, D c . After weakening, the dynamic strength of the fault during slip, T fd , is given by: We assign μ d = 0.1 and use a constant value of D c = 0.8 m. The rupture continues to propagate as long as τ s locally exceeds T fs and a fault continues to slip as long as sufficient strain energy is available. Note that τ s at the rupture front is typically higher than the initial τ s , so statically stronger parts of a fault may fail after the rupture initiates elsewhere. Table 3 summarizes average characteristics of the earthquakes in each scenario. As pore fluid pressure (P f ) increases from low to very high, the moment magnitude (M w ) decreases, as do mean cumulative slip, peak slip rate (PSR), mean dynamic stress drop (Δτ s ), and rupture velocity (Vr). This reflects our here chosen set-up, in which both shear and effective normal tractions scale inversely with P f . M w of the earthquakes in Scenarios 1 and 2 are unrealistically large, which supports the conjecture by Saffer and Tobin (2011) that pore fluid is likely overpressured everywhere along the seismogenic megathrust. Further details about Scenarios 1 and 2 are given in Appendix C. M w for the earthquakes in Scenarios 3-6 are reasonable for a rupture area the size of the Sumatra earthquake (Strasser et al., 2010), thus, we focus on the results for these four scenarios. Videos of the slip rate evolving along the megathrust during each of these scenarios are available by link from Appendix D.

Earthquake Source Characteristics
In all four scenarios, an initially crack-like rupture develops into sharp, boomerang-shaped rupture pulses propagating along-arc on the megathrust. Each pulse consists of multiple rupture fronts, which are caused by reflected waves and head waves generated at structural interfaces and the complex free surface (Huang et al., 2014). We note that pulse-like rupture is here not caused by self-healing due to the dynamics of fault strength (Gabriel et al., 2012), but due to geometric constraints (Weng & Ampuero, 2019). Figure 5 compares slip, PSR, Δτ s , and Vr on the megathrust at the end of the earthquakes in Scenarios 3-6. All three splay faults incorporated into the base model are dynamically activated in all scenarios. In general, they slip an order of magnitude less than the megathrust.
The magnitude of P f inversely affects average cumulative slip, while its gradient (sublithostatic or lithostatic) influences the slip distribution on the megathrust ( Figure 5). As P f increases from high in Scenario 3 to very high in Scenario 4, mean slip decreases from 26 to 8 m. This is reflected in the decrease in earthquake moment magnitude from M w 9.3 in Scenario 3 to M w 9.0 in Scenario 4. The slip is similarly distributed in both scenarios, with maximum slip in the middle of the fault in the downdip direction. Slip is highest in the center of the fault along strike. Likewise, as P f increases from high in Scenario 5 to very high in Scenario 6, mean slip decreases from 36 to 10 m and moment magnitude decreases from M w 9.4 to M w 9.1. Scenario 6. As shown for Scenario 4 here, the sublithostatic P f gradient in Scenarios 1-4 causes the stresses to increase with depth (to the left). Stress magnitude ranges widen progressively from Scenario 3 to Scenario 2 to Scenario 1, but the pattern is the same. As shown for Scenario 6 here, the lithostatic P f gradient in Scenarios 5 and 6 causes the stresses to be constant with depth. The stress magnitudes are larger in Scenario 5, but remain constant with depth. Below 23 km, the differential stress is tapered to zero in all scenarios (not shown).

Figure 4.
Initial shear traction (τ s ) and effective normal traction ( ′ ) on the megathrust in Scenarios 1-6. For each fault image, the shallowest part of the megathrust, near the seafloor, is to the left and the deepest part at 50 km depth is to the right. Note the depth-dependent ′ in Scenarios 1-4 with sublithostatic P f gradients applied vs. the nearly constant ′ in Scenarios 5 and 6 with lithostatic P f gradients. Both τ s and ′ vary with the non-planar fault geometry up to ≈5 MPa.
Mean slip and M w are similar in scenarios with the same P f magnitude, but different depth gradients, for example, in Scenarios 3 and 5 and in Scenarios 4 and 6. However, in Scenarios 5 and 6, in which the P f gradient is lithostatic and effective normal stress is constant with depth, maximum slip is shifted updip relative to the locations of maximum slip in Scenarios 3 and 4, in which the P f gradient is sublithostatic and constant effective normal stress increases with depth. Slip to the trench only occurs in Scenario 5, and slip is limited at the trench in Scenarios 3, 4, and 6. We discuss this further in Section 5.1 and Appendix E).
As with cumulative slip, peak slip rate PSR in these scenarios decreases as P f magnitude increases and the P f gradient influences its distribution along the megathrust. Mean PSR is 10 m/s in Scenario 3 with high P f and 5 m/s in Scenario 4 with very high P f . Mean PSR is 11 m/s in Scenario 5 with high P f and decreases to 6 m/s in Scenario 6 with very high P f . Comparing across P f gradients, we see that Scenarios 3 and 5 and Scenarios 4 and 6 have similar mean PSR values, but maximum PSR occurs below 35 km depth in Scenarios 3 and 4 and above 15 km in Scenarios 5 and 6. Thus, relative to depth-dependent effective normal stress under sublithostatic P f conditions, assuming a lithostatic P f gradient resulting in constant effective normal stress with depth shifts maximum PSR updip ( Figure 5). In addition, more of the megathrust experiences high PSR in Scenario 6 relative to Scenario 4, though maximum values are lower in Scenario 6.
We measure the mean dynamic stress drop (Δτ s ) as the average change in shear traction (τ s ) from the initial value to the dynamically reached value at the end of the earthquake. As for mean slip and PSR, P f has an inverse relationship with mean Δτ s . Mean Δτ s is 8 MPa in Scenario 3 and 7 MPa in Scenario 5, and 3 MPa in both Scenarios 4 and 6. The distribution of Δτ s varies with the P f depth gradient. In Scenarios 3 and 4, Δτ s is larger along the deeper fault, reaching values of 15 and 7 MPa, respectively, below 30 km depth ( Figure 5). In Scenarios 5 and 6, Δτ s is relatively constant along the central fault in the downdip direction. The highest values are farther updip near 20 km depth, at 12 and 5 MPa in these scenarios, respectively. In all scenarios, Δτ s is largest along the central portion of the fault along strike.
In contrast to the other earthquake characteristics, there is little variation in the distribution of Vr with P f depth gradient. However, an increase in P f magnitude overall causes a decrease in average rupture velocity, Vr, from 3,025 m/s in Scenario 3 to 2,370 m/s in Scenario 4 and from 3,206 m/s in Scenario 5 to 2,624 m/s in Scenario 6. Mean Vr is lower in Scenario 3 relative to Scenario 5, and lower in Scenario 4 relative to Scenario 6, suggesting that average Vr increases under conditions of constant vs. depth-dependent effective normal stress.
In all scenarios, average Vr is sub-Rayleigh relative to the lower velocity subduction channel surrounding the megathrust slip interface (V s = 3,500 m/s, Table 1). While Vr is below Rayleigh wave speed across most of the megathrust in all scenarios, exceptions of supershear rupture appear (a) propagating updip from the hypocenter at close to P wave speed triggered by energetic nucleation and (b) in the form of localized and relatively slow supershear fronts excited before the sub-Rayleigh rupture front at several isolated locations. In Scenario 5, where Vr is highest out of all scenarios, at these isolated locations Vr ≈ 70% of P wave speed. Vr that exceeds the S wave speed, but remains lower than the P wave speed, agrees with inferences and modeling for earthquake rupture in damaged fault zones (Bao et al., 2019;Harris & Day, 1997;Huang et al., 2016;Oral et al., 2020).

Post-Earthquake Stress Field
The dynamic rupture model utilized in these scenarios permits investigation of the post-earthquake absolute stress field. We compare principal stress orientations and relative magnitudes along a cross-section of the central part of the rupture in Scenarios 3-6 (see inset in Figure 6a). Figure 6a shows the orientations of the principal stresses (σ 3 < σ 2 < σ 1 , compression is negative) before the earthquake for all scenarios and Figure 6b shows the orientations after dynamic earthquake rupture in Scenario 4. The post-earthquake stress orientations for Scenarios 3-6 are shown in Figure F1. We summarize the post-earthquake stress orientations for all scenarios in stereonets focused on the hanging wall and footwall regions close to the fault in Figure 6c. We compare the mean

Table 3
Earthquake Characteristics Averaged Across the Megathrust Figure 5. For Scenarios 3-6: cumulative slip, peak slip rate (PSR), dynamic stress drop (Δτ s ), and rupture velocity (Vr) on the megathrust. For each fault image, the shallowest part of the fault is to the left and the deepest part (at 50 km depth) is to the right. A version with alternative colorbar limits that are set for comparison across scenarios is included as Figure C2.
In all scenarios, σ 2 and σ 3 have similar mean apparent rotations and rotate more than the minimum principal stress, σ 1 . The mean principal stress rotations in the hanging wall summarized in Table 5 vary with the magnitude of pore fluid pressure (P f ). As P f increases from Scenario 3 to Scenario 4 and from Scenario 5 to Scenario 6, mean rotations of each principal stress decrease in accordance with decreasing stress drop. Scenarios 4 and 6 have very similar apparent rotations for each principal stress, suggesting that the choice of P f depth gradient does not affect the amount of rotation when the P f magnitude is very high (97% of the lithostatic pressure, σ v ). Such similarity is not apparent when comparing Scenarios 3 and 5. Mean rotations in Scenario 5 are the largest of all scenarios, which we attribute to the high fault slip at the trench in this scenario.
To better understand the post-earthquake stress field, we also consider the effective principal stress magnitudes relative to one another. This is important to the stress rotation analysis, because magnitudes of two principal stresses that move closer to one another approach the condition for switching orientations, allowing for a larger amount of heterogeneity in the post-earthquake stress field. Figure 7 shows the maximum differential stress, , before and after the dynamic earthquake ruptures in Scenarios 3-6. Prior to each earthquake, the distributions of ′ 13 depend on the gradient in P f . Scenarios 3 and 4 have the same depth-dependent pattern of  and ′ 23 = ′ 2 − ′ 3 . As P f increases from Scenario 3 to Scenario 4 and from Scenario 5 to Scenario 6, pre-earthquake ′ 13 averages in the hanging wall decrease by ≈20 MPa. In each scenario, ′ 12 equals ′ 23 before the earthquake, as σ 2 is initially set to be halfway between σ 3 and σ 1 . Pre-earthquake, the magnitudes of these differential stresses differ from Scenario 3 to Scenario 4 and from Scenario 5 to Scenario 6 by ≈10 MPa.
In the plots of the post-earthquake ′ 13 distributions in Figure 7, contours indicate the amount and direction In all scenarios, there are larger changes in average ′ 23 than in average ′ 12 due to the larger coseismic decrease in the magnitude of ′ 3 relative to the decreases in ′ 1 and ′ 2 ( Table 5). The closeness of ′ 2 and ′ 3 before the earthquake therefore controls the amount of apparent post-seismic stress rotation here, and how likely these two principal stresses are to switch locations. In contrast, ′ 2 and ′ 1 have less apparent rotation, making them less likely to swap orientations.

Table 6
Differential Stress Before and After the Earthquake a Figure 7. Cosesimic change in maximum effective differential stress ( ′ 13 ) (a) before the earthquake in Scenarios 3 and 4, (b) after the earthquake in Scenario 3, (c) after the earthquake in Scenario 4, (d) before the earthquake in Scenarios 5 and 6, (e) after the earthquake in Scenario 5, and (f) after the earthquake in Scenario 6. Contours show change in ′ 13 from pre-to post-earthquake. Location is as shown in inset in Figure 6.

Discussion
We present six earthquake scenarios that vary in P f magnitude and depth gradient in order to explore the dynamic effects of different pre-earthquake P f levels and distributions in subduction zones. The model structure and input are consistent with conditions for the 2004 Sumatra-Andaman earthquake, using a base model following (Uphoff et al., 2017). We first discuss how the scenario earthquakes reflect observations of that event, as well as more general observations of earthquakes along megathrusts. Then, we discuss inferences from these scenarios relevant to fault mechanics. We analyze further the stress rotations from before to after these scenario earthquakes and compare them to observations following the 2004 Sumatra earthquake.

Earthquake Characteristics
Pre-earthquake conditions are not easily constrained by observations, here along the Sumatra-Andaman trench or elsewhere in the world. However, the observational matching of the base model by Uphoff et al. (2017) used here gives an ideal starting point to explore the effects of P f on earthquake dynamics. In addition, the 3-D physics-based forward modeling approach unifies the pre-earthquake conditions together with the earthquake dynamics to arrive at physically consistent earthquake characteristics, a capability of large-scale and geometrically complex computational models highlighted by Ulrich et al. (2022).
To first order, Scenarios 3 and 6 produce earthquakes with moment magnitudes similar to those inferred for the Sumatra earthquake of M w 9.1-9.3 (Shearer & Bürgmann, 2010), while the Scenario 4 earthquake is just below this range at M w 9.0 and the Scenario 5 earthquake is just above this range at M w 9.4 (Table 3). Maximum slip values from kinematic source inversions compiled by Shearer and Bürgmann (2010) range up to a maximum value of ≈35 m, suggesting that slip in the Scenario 5 earthquake, which averages 36 m, is too large. Seno (2017) estimates a mean stress drop of 3 MPa for this earthquake, which is matched by those for Scenarios 4 and 6. In contrast, Scenarios 3 and 5 have mean dynamic stress drops that are more than twice this value. The mean rupture velocities in Scenarios 4 and 6, respectively 2,370 and 2,624 m/s, are similar to the rupture velocity of 2,500 m/s inferred by Ammon et al. (2005) for the 2004 earthquake. In contrast, Scenarios 3 and 5 both have mean Vr exceeding 3,000 m/s. Seno (2017) estimates a subducted sediment thickness of 1.57 ± 0.12 km near Simeulue, in the southern region of the 2004 earthquake, which is high in comparison with other subduction zones. Correlation between subducted sediment thickness, stress drop and P f by Seno (2017) suggests that P f should be high and stress drop should be low in this earthquake, as in both Scenarios 4 and 6. This highlights the earthquakes in Scenarios 4 and 6 as more realistic.
Scenarios 4 and 6 both have very high P f at 97% of the lithostatic stress (σ v ), but differ in the way that P f is acting on the curved fault system. In Scenario 4, P f follows a sublithostatic depth gradient and the effective normal traction ( ′ ) increases with depth. In Scenario 6, following theoretical work by Rice (1992), P f follows the lithostatic gradient, maintaining a constant difference to σ v . As a result, ′ is close to constant with depth along most of the megathrust (varying only by up to 5 MPa due to variations in fault geometry). The good performances of Scenarios 4 and 6 relative to observations of the 2004 Sumatra earthquake suggest that megathrust earthquakes may occur under very high pre-earthquake P f magnitudes, resulting in low ′ magnitudes. Scenario 6 emerges as the scenario that best matches observations, as Scenario 4 has lower slip that results in a M w 9.0 event, smaller than the M w 9.1-9.3 2004 earthquake (Shearer & Bürgmann, 2010). This furth suggests that megathrust earthquakes may occur under conditions of a lithostatic P f depth gradient, resulting in relatively constant ′ along the megathrust.
These scenarios also are representative of variable conditions that may be present along a single megathrust at the same point in time, due to spatial variations in P f magnitude and/or gradient. Such variations in P f are one possible mechanism of conceptual seismic asperities, inducing heterogeneity in dynamic fault motion (Bürgmann, 2018;Lay et al., 2012). Sediments and high P f have been proposed as important mechanisms aiding stable sliding along geometric, frictional and rheological barriers, while (less effectively) thermal pressurization may provide a mechanism for stress-roughening slip events (Barbot, 2019;Gabriel et al., 2020;Perry et al., 2020;Wibberley & Shimamoto, 2005). Our presented scenarios serve as building blocks for future along-arc heterogeneous models that may be calibrated with site-specific friction and pore-fluid measurements to constrain dynamically plausible megathrust strength and P f gradients. For example, we find that very high P f leading to constant effective normal stress with depth produces a stress drop on the megathrust that is nearly constant with depth and pushes PSR updip on the megathrust. Also, earthquake magnitude and mean cumulative slip are larger for an equal or lower mean stress drop under these conditions. For a given subduction zone or megathrust event, such detailed conditions may be constrained by geodetic, geological, or tsunami observations (e.g., Ulrich et al., 2022).
High or very high P f that follows the lithostatic gradient favors higher slip at shallower depths, thus increasing the importance of near-trench strength and constitutive behavior in determining megathrust hazard. Widespread and high amplitude slip to the trench only occurs in Scenario 5, and slip is limited at the trench in Scenarios 3, 4, and 6. In all scenarios, near-trench behavior is influenced by the choice of on-fault cohesion, c, which is used as a proxy for near-trench behavior that we do not model explicitly here, such as velocity-strengthening during slip in shallow sediments (e.g., Kaneko et al., 2008) and the energy lost to rock yielding around the megathrust (off-fault plasticity, e.g., Gabriel et al., 2013). c is the same in all scenarios, but its relative contribution to the static fault strength increases as the P f increases and ′ magnitude decreases (Equation 2, Figure 4). Models that aim to capture natural co-seismic near-trench processes (e.g., Dunham et al., 2011;Lotto et al., 2019;Ma, 2012;Ma & Nie, 2019;Ulrich et al., 2022) can further discriminate governing factors of near-trench behavior (see also Appendix E). Specifically, Ulrich et al. (2022) focus on near-trench behavior during the 2004 Sumatra earthquake and its influence on the subsequent Indian Ocean tsunami.
Next, we look to general observations of stress drop from earthquakes on the subduction interface to further decipher between scenarios. Allmann and Shearer (2009) report depth-dependent stress drops when data is considered separately by region. Uchide et al. (2014) find an increasing stress drop from 30 to 60 km depth in a spectral decomposition analysis of smaller events occurring before the 2011 Tohoku earthquake. However, Bilek and Lay (2018) and Denolle and Shearer (2016) report very weak correlation between stress drop and depth. Abercrombie et al. (2021) re-evaluate previous studies based on the spectral decomposition method and show that when trade-offs between attenuation and depth-dependent sources are accounted for, the correlation between stress drop and depth from previous studies decreases and, in some cases, disappears altogether. We determine the dynamic stress drop on the megathrust in each scenario, which differs slightly from these observationally inferred values, but remains well within observational and methodological uncertainties. We find that dynamic stress drop varies more with depth in Scenarios 3 and 4 (up to 15 MPa), due to the depth-dependent effective normal traction resulting from the sublithostatic P f gradient ( Figure 5). In contrast, stress drop varies up to only 7 MPa in Scenarios 5 and 6, where effective normal traction is relatively constant along the megathrust resulting from the lithostatic P f gradient. Thus, a correlation between stress drop and depth is more consistent with high P f following a sublithostatic gradient, while a low dependence of stress drop on depth is more consistent with high P f following a lithostatic gradient. Should these end-member conditions be present in different locations along a single megathrust, deciphering a dependence of stress drop on depth observationally will be difficult. On the other hand, well-constrained observations of depth-dependent vs. depth-constant stress drops of small events may differentiate between locations of sublithostatic (e.g., Scenarios 1-4) vs. lithostatic (e.g., Scenarios 5 and 6) P f gradients along megathrusts.
Under a lithostatic P f gradient, the effective normal stress is constant and the effective normal tractions ( ′ ) are relatively constant, but variations of ≈5 MPa still arise due to variations in fault geometry. Bletery et al. (2016) attribute the location and extent of the 2004 Sumatra earthquake rupture to a region of relatively homogeneous megathrust shear strength. Homogeneity of ′ , and therefore of fault shear strength in these scenarios, is promoted by high P f that follows the lithostatic gradient with depth. Such homogeneous shear strength is more likely to be exceeded simultaneously over large areas, leading to the large earthquakes observed in subduction zones. However, it is interesting to note that conditions of relatively homogeneous ′ and shear strength may actually emphasize the influence of geometry on earthquake behavior, as geometry becomes the main control on shear strength variation along the megathrust. Both effects may be explored in future work focusing on variations in megathrust geometric complexity and cycles of fault slip (e.g., Perez-Silva et al., 2021) and by relaxing our assumption of a constant shear to effective normal traction ratio.

Inferences From These Scenarios Relevant to Fault Mechanics
Here, we consider the scenarios in light of inferences about fault mechanics, beginning with the initial shear traction (τ s ) on the fault, then discussing effective normal traction ( ′ ) magnitudes and variation with depth. τ s scales with ′ from scenario to scenario and the distribution of τ s / ′ is the same in all scenarios ( Figure B1). A static friction coefficient of 0.4 is applied in all scenarios.
From force-balance studies, Lamb (2006) finds that the crust above 7 out of 9 studied subduction zones sustains an average τ s of 7-15 MPa. This includes Sumatra, with an average τ s of 15.2 MPa (Lamb, 2006, Table 5), which is similar to the mean τ s prior to rupture on the megathrust in Scenarios 3 and 5. Brodsky et al. (2020, Figure 6) constrain τ s on the shallow part of the Tohoku megathrust prior to the 2011 Tohoku earthquake at ≈1.7 MPa using a friction coefficient derived from low-velocity friction experiments. Yao and Yang (2020) find the shear strength of the megathrust that ruptured in the 2012 Nicoya earthquake to be less than 7.5 MPa on average. In combination with observed low-stress drops of subduction megathrust events (Sibson & Rowland, 2003), low dynamic shear stresses during earthquake rupture (e.g., less than 1 MPa, Choy & Boatwright, 1995;Pérez-Campos & Beroza, 2001) also support low τ s on megathrusts prior to earthquakes, although this may include additional weakening from a variety of dynamic effects (Gao & Wang, 2014).
In this suite of six scenarios, more reasonable earthquakes emerge at higher pre-earthquake P f magnitudes and average initial τ s values in Scenarios 3-6 range from 5 to 11 MPa (Table 2). Thus P f higher than approximately 93% of the lithostatic gradient is consistent with inferences of low initial shear stress on the megathrust. As suggested by the analysis in Section 5.1, Scenarios 4 and 6 produce the most realistic earthquakes, supporting P f averaging at 97% of the lithostatic stress (σ v ) and consistent with mean τ s on the megathrust of 4-5 MPa. There are exceptions to inferences of low initial τ s , however. Lamb (2006) estimates values of 18.3 and 36.7 MPa on the Chile and Tonga megathrusts, respectively, while depth-dependence is inferred for the Tohoku and northern Hikurangi megathrusts with values ranging up to 80 MPa (Gao & Wang, 2014;. These values are more consistent with Scenarios 3 and 5. In all scenarios, the megathrust is moderately strong, with a static friction coefficient of 0.4. However, the low shear strengths (T fs , Equation 2) of the megathrust in the preferred scenarios can be used to classify the megathrust as weak. The megathrust also is dynamically weak, with friction dropping to 0.1 during sliding.
In these scenarios, high P f magnitude leads to low effective maximum differential stress and low τ s along the megathrust. However, this can occur independently of P f , for example, from absolute principal stresses that are close to one another in magnitude. We assume a least compressive principal stress, σ 1 , in our scenarios that is close to σ v . The other two principal stresses must be larger in magnitude in a thrust faulting regime, but are more difficult to constrain. σ 3 could vary from what we choose, which would then change τ s on the megathrust as well as the average τ s associated with a particular P f . More complicated stress conditions also are likely. For example, we choose to set σ 2 midway between σ 1 or σ 3 , but this is not necessarily the case in nature. In addition, principal stress magnitudes may vary in magnitude or orientation along the megathrust, both laterally and with depth. Past earthquakes may leave heterogeneous shear tractions on the megathrust and P f likely varies spatially in the vicinity of the megathrust (Heise et al., 2017). Close to the fault, there is field evidence of stress rotations within the damage zone that vary the principal stress orientations from those in the remote field (Faulkner et al., 2006) and this condition is supported by theory (Rice, 1992). It will be interesting to relate stress complexity with P f and additional along-arc heterogeneity in future work.

Off-Fault Results
It has been suggested that principal stress rotations are promoted by complete or near-complete stress drops that permit principal stresses to swap orientations (Brodsky et al., 2020(Brodsky et al., , 2017X. Wang & Morgan, 2019). However, by connecting 2-D stress rotations to the ratio of stress drop over pre-earthquake deviatoric stress magnitude, Hardebeck (2012Hardebeck ( , 2015 shows that partial stress release may generate moderate rotations. Scenarios 3 and 5 experience the largest rotations, but have larger initial differential stresses and larger post-earthquake differential stresses as well. The larger rotations in these scenarios appear to scale with fault slip and stress drop, both of which are larger than in Scenarios 4 and 6. X. Wang and Morgan (2019) attribute observed changes in stress orientations following the 2011 Tohoku earthquake to rapid weakening of a statically strong fault with μ s in the range of 0.3-0.6. K.  attribute rotations to a weak megathrust, with a low effective friction coefficient (0.032) and low shear stress in the forearc leading to low shear traction on the megathrust. These theories are compatible with one another, if the megathrust is considered to be statically strong, but dynamically weak, in terms of its dynamic friction coefficient, and if P f is high. This is supported by the scenarios presented here, with μ s = 0.4 and μ d = 0.1.
None of the scenarios results in a complete stress drop and yet we find that the postseismic stress field supports a variety of potential aftershock focal mechanisms. In all scenarios, σ 3 rotates toward parallel with megathrust strike and its plunge remains more or less unchanged, while the plunge of σ 2 increases and the plunge of σ 1 decreases. This post-seismic stress state supports a variety of aftershock mechanisms, including strike-slip faulting where σ 1 plunges more shallowly relative to σ 2 , and reverse faulting where σ 2 plunges more shallowly relative to σ 1 . Of 13 M w 6 or larger aftershocks with focal mechanisms solutions in the GCMT catalog (Dziewonski et al., 1981;Ekström et al., 2012) occurring along the central rupture within five years of the 2004 Sumatra mainshock (through 27 December 2009), 8 are reverse and 5 are strike-slip. We define the central rupture here as the region from 5° to 9° latitude, 91°-97.3° longitude, and 0-50 km depth, corresponding to the location of the fault-perpendicular slice in Figure 6. Out of 125 M w 5 or larger aftershocks occurring within 1 month of the mainshock in the same region, 63 have strike-slip focal mechanisms, while 29 have reverse, 31 have normal mechanisms and 2 cannot be categorized.
At Sumatra, Hardebeck (2012) finds rotations of the maximum compressive principal stress, which we call σ 3 , relative to the megathrust and in the 2-D plane perpendicular to the megathrust, to be up to ≈42° and increasing from South to North. Along the central rupture (zone B in Hardebeck, 2012), average σ 3 rotation is 26° ± 13°. Using the 2-D solution proposed by Hardebeck and Hauksson (2001), the ratio of the mean earthquake stress drop to the magnitude of the deviatoric stress, Δτ s /σ dev , can be estimated as a function of the pre-earthquake angle of σ 3 to the megathrust and its rotation. At Sumatra specifically, Hardebeck (2012) finds that this ratio varies from 0.6 along the southern part of the rupture to 0.8 along the central and northern part of the rupture. This implies that 60%-80% of the pre-earthquake deviatoric stress magnitude along the megathrust was relieved by the earthquake. The apparent rotations of σ 3 along the central rupture in the scenarios presented here (Table 5) are of similar magnitudes to those determined from data (Hardebeck, 2012), ranging from 36° to 55°, but are predominantly in the horizontal plane. We find similar average ratios of Δτ s to σ dev in these scenarios, of 0.6 in Scenarios 4, 5, and 6 and of 0.7 in Scenario 3, when taking σ dev from these 3-D models as the difference between σ 3 and the mean stress. We do not see correspondence between differences in Δτ s /σ dev and the amount of σ 3 rotation (Table 5), but note that this analysis is not directly comparable to the 2-D analysis by Hardebeck (2012), as σ 3 rotates out of the plane perpendicular to the megathrust in these scenarios.
Post-earthquake stress and aftershock focal mechanism heterogeneity would be further promoted in a model incorporating a heterogeneous initial stress field. In these scenarios, a laterally constant, depth-dependent regional stress tensor is applied, so P f and the resulting effective stress field are the same near to and far from the megathrust before the earthquake. Such similar on-and off-fault stresses are not likely in nature. Away from the megathrust, secondary faulting, the earthquake history, and material contrasts likely produce stress heterogeneities (van Zelst et al., 2020). Heterogeneity in the magnitude of the effective intermediate principal stress, ′ 2 , relative to the maximum and minimum effective principal stresses also would contribute to aftershock heterogeneity, by making it easier for different faulting regimes to be activated. For example, as we note in Section 4.2, the magnitude of ′ 2 relative to the other two effective principal stresses controls the ability for ′ 2 to switch places with ′ 1 or ′ 3 , thus affecting postseismic stress rotations. In addition, dynamic effects that decouple conditions on-and off-fault, such as thermal pressurization (Noda, 2008;Noda et al., 2009) during which P f increases rapidly due to reduced pore pressure diffusion in the fault zone during slip, may allow low effective normal tractions on the megathrust, even while a different stress state persists away from the fault. Considering more complex initial stress conditions off the fault and decoupling on-and off-fault stresses are clear next steps for this work.

Conclusions
We analyze the effects of pore fluid pressure (P f ) magnitude and gradient on pre-earthquake stress conditions and earthquake dynamics using 3-D high-performance computing enabled, physics-based dynamic rupture models that permit geometrically complex faults. The six scenarios presented, based on the 2004 M w 9.1 Sumatra-Andaman earthquake, have P f that varies from hydrostatic to lithostatic under sublithostatic vs. lithostatic gradients. These result, respectively, in either depth-dependent or constant effective normal stress near the megathrust and splay faults. As P f increases in these scenarios, moment magnitude, cumulative slip, peak slip rate, dynamic stress drop, and rupture velocity all decrease. A lithostatic P f gradient causes relatively constant effective normal tractions on the megathrust, moves peak slip and peak slip rate updip, and produces a more constant stress drop across the megathrust. This is consistent with theoretical analysis and observations inferring that the stress drops of smaller earthquakes in subduction zones are only weakly depth-dependent.
In comparison with a range of observations, we identify two preferred scenarios that both support the presence of very high coseismic P f on average over the ruptured area (here 97% of the lithostatic pressure). These have low mean shear and effective normal traction magnitudes of 4-5 and 22 MPa, respectively. The mean dynamic stress drop for these two scenario earthquakes is 3 MPa and the mean rupture velocity is 2,400-2,600 m/s, similar to observations of the 2004 Sumatra-Andaman earthquake. Although comparison with observations of the 2004 earthquake cannot conclusively differentiate between these two preferred scenarios, a lithostatic P f gradient, which causes constant normal stress near the megathrust, may be the theoretically more plausible condition under very high P f magnitudes. On weak megathrusts, in terms of low static shear strength and low dynamic friction during rupture, where P f follows the lithostatic gradient, near-trench strength and constitutive behavior are crucially important for megathrust hazard, as peak slip and peak slip rate occur at shallower depths.
Mean apparent rotations of the principal stresses in the hanging wall decrease as P f magnitude increases, but do not vary with P f gradient. Scenarios with the largest rotations have larger initial differential stress and larger post-earthquake differential stress as well. The larger rotations in these scenarios scale with fault slip and stress drop. Along the central rupture, maximum compressive stress rotations in the hanging wall average 36° ± 18° toward trench-parallel in the two preferred scenarios and the minimum principal stress rotates from near-vertical toward a shallower plunge. This post-earthquake stress field is consistent with the heterogeneous aftershocks observed following the Sumatra earthquake.
Variations in P f are one possible mechanism of conceptual seismic asperities, and our analysis may serve as guidance for future along-arc heterogeneous models. In addition, this work has implications for tsunami hazard, as the P f gradient is shown to influence the location of maximum slip and slip rate. Under conditions of a lithostatic P f gradient, relatively constant effective normal tractions downdip along the megathrust push maximum slip and slip rate toward the surface.

Appendix A: Model Mesh Resolution
Dynamic rupture simulations must resolve the cohesive zone width Λ, which spans the part of the fault across which shear stress decreases from its static to its dynamic value. In heterogeneous dynamic rupture simulations, Λ can vary considerably across the fault in dependence on initial stress, frictional properties, and propagation distance. Since Λ also changes dynamically across the fault, the number of elements per median Λ can also vary significantly across the fault for a given simulation. Selected findings from Wollherr et al. (2018) allow for a better understanding of how the numerical accuracy of SeisSol's ADER-DG scheme, in resolving on fault time-dependent parameters, is affected by mesh size and polynomial degree. By comparing the rupture arrival time, peak slip rate (PSR) time, final slip and the PSR averaged across 363 receivers with respect to a reference solution, Wollherr et al. (2018) show that errors are globally decreasing with mesh refinement and increasing polynomial degree. SeisSol resolves shear and normal stress and effective friction according to a friction law everywhere at the fault at (p + 2) 2 Gaussian quadrature points inside each fault element triangle, with p being the polynomial degree (and p = 3 in this study leading to fourth order accuracy in space and time, as measured in the L2 norm, of the ADER-DG scheme for seismic wave propagation).
In this study, we ensure that we resolve the median Λ, estimated at 1 km for Scenario 3, 1.6 km for Scenario 4, 0.9 km for Scenario 5, and 1.5 km for Scenario 6. In assessing sufficient resolution, we follow Day et al. (2005) and Wollherr et al. (2018). Day et al. (2005) define a dynamic rupture solution to be sufficiently close to the reference solution once the RMS errors reached the following thresholds: lower than 0.2% for rupture arrival time, lower than 7% for PSR and lower than 1% for final slip. We calculate Λ as the difference in distance between the rupture front arrival time and the first point in time at which shear stresses reach their dynamic value across the fault. The minimum Λ (approximated at the fifteenth percentile of all measured) varies across the scenarios as follows: 346 m in Scenario 3, 540 m in Scenario 4, 469 m in Scenario 5 and 627 m in Scenario 6. By analyzing Scenarios 3 and 6, with the longest and shortest Λ, we find that the errors for rupture arrival range from 0.09% to 0.20% and for final slip range from 0.68% to 1.2% across these four scenarios, which are sufficiently small with respect to the findings by Day et al. (2005). The expected errors for PSR are higher, ranging from 8.9% to 17%, above the 7% recommended by Day et al. (2005), however, Ramos et al. (2021) verify with higher resolution models that even with expected errors above 7% for PSR, megathrust slip is not affected in comparable SeisSol dynamic rupture models.

Appendix B: Prestress Ratio and On-Fault Frictional Cohesion
The relative prestress ratio, R, is the ratio of the fault stress drop (τ s − T fd ) to the breakdown strength drop (T fs − T fd ), where τ s is the initial shear traction, T fs is the static fault strength and T fd is the dynamic fault strength during sliding (Aochi & Madariaga, 2003). R varies along the megathrust with the non-planar fault geometry ( Figure B1), but is nearly the same across all scenarios since ∕ ′ is constant across all scenarios. The exception to this is with respect to the on-fault frictional cohesion, c. c is similar across all scenarios, but contributes differently to T fs in each scenario and these changes R slightly from scenario to scenario, particularly at shallow depths (see also Appendix E).
Cohesion, c, depends on local mineralogy and lithology. However, c is used here to limit slip in the absence of near-trench behavior, using the lowest value that restricts unrealistic slip and rupture dynamics (e.g., occurrence of supershear rupture) at the trench. We find this to be c = 0.4 MPa below 10 km depth and increasing linearly to 15 MPa at 0 km depth ( Figure B2). We tested two alternative c gradients from 0.4 MPa below 10 km to maxima of 1 and 10 MPa at z = 0, which lead to unrealistic near-surface behavior. As the base model used here does not capture the constitutive behavior of shallow sediments in the near-trench region, we do not draw conclusions about near-trench behavior or about realistic c values from these scenarios (see also Appendix E). Ulrich et al. (2022) take the work in this direction by incorporating slip-strengthening behavior near the seafloor, as well as off-fault plasticity, into models of the 2004 Sumatra-Andaman earthquake. Figure B1. (a) The ratio of the initial shear traction to effective normal traction ( ∕ ′ ) varies depending on the megathrust orientation relative to the local stress tensor, but the distribution on the megathrust is the same across all scenarios. (b) The prestress ratio, R, is shown here for Scenario 4, but is similar in all scenarios. Figure B2. Blue line is on-fault frictional cohesion, c, which is set to 0.4 MPa below 10 km depth and increases linearly to 15 MPa at 0 km depth. Due to topography, the intersection of the fault and the seafloor ranges between 3 and 5 km depth, so maximum c values on the megathrust and splay faults at the seafloor range from 8 to 11 MPa. The gray line shows this intersection between fault and seafloor on average.
When the fault is in tension and effective normal stress equals zero, the fault strength is equal to c. This is because tensile stresses are treated in SeisSol to prevent fault opening following a standard approach in the dynamic rupture community (Harris et al., 2018(Harris et al., , 2009). This procedure treats tension on the fault the same as if the effective normal stress equals zero.

Appendix C: Scenarios 1 and 2 Earthquakes
Slip, peak slip rate, dynamic stress drop, and rupture velocity are shown in Figure C1 for Scenarios 1 and 2, which have low and moderate P f , respectively. Slip proceeds to the trench in Scenario 5 and reaches maximum values there, which is clearly different from Scenarios 3, 4, and 6 ( Figure 5 and Figure C2). A similar difference between shallow slip in Scenario 4 and Scenario 6 also is visible in Figure 5. These differences are due not only to P f magnitude and depth gradient, but also to the contribution of the applied on-fault cohesion, c, to static fault strength, T fs (see also Appendix B). In all scenarios, c is constant below 10 km depth and linearly increases toward the surface above, contributing to T fs according to Equation 2. The influence of c on T fs increases as P f increases and the magnitude of ′ decreases. As a result, closeness to failure varies near the seafloor in all scenarios. Fault strength is overcome at the trench only in Scenario 5, while slip is restricted along the top of the fault in Scenarios 3, 4, and 6. This contrast is important because it highlights both that the influence of c on slip behavior at the trench increases as P f increases and c becomes a larger component of T fs , and that near-trench slip is encouraged by very high P f following a lithostatic gradient that causes conditions of constant ′ along the megathrust and pushes maximum slip and slip rate closer to the trench. In these scenarios, c is defined as the strength of the fault in the absence of τ n (Equation 2) and is used as a proxy for near-trench behavior that we do not model explicitly here, including the energy lost to damage around the megathrust (off-fault plasticity, e.g., Gabriel et al., 2013) and velocity-strengthening of the fault in shallow sediments (e.g., Kaneko et al., 2008). Further study of slip behavior at the trench requires that the appropriate physical processes near the seafloor are incorporated into the model (e.g., Dunham et al., 2011;Lotto et   slip strengthening and off-fault plasticity of lithified shallow sediments into coupled earthquake-tsunami models of the 2004 Sumatra earthquake and Indian Ocean tsunami to study near-trench slip, seafloor displacement, and tsunami genesis using a coupled tsunami model. Figure F1 shows the post-seismic stress field for all scenarios. While the rotation directions are similar in all scenarios, the amount of rotation is larger in Scenarios 3 and 5 than in Scenarios 4 and 6. Stereonets are included in the main text ( Figure 6).