Sequences of seismic and aseismic slip on bimaterial faults show dominant rupture asymmetry and potential for elevated seismic hazard

We perform numerical simulations of sequences of earthquake and aseismic slip on planar rate and state faults separating dissimilar material within the 2-D plane strain approximation. We resolve all stages of the earthquake cycle from aseismic slip to fast ruptures while incorporating full inertia effects during seismic event propagation. We show that bimaterial coupling results in favorable nucleation site and subsequent asymmetric rupture propagation. We demonstrate that increasing the material contrast enhances this asymmetry leading to higher slip rates and normal stress drops in the preferred rupture propagation direction. The normal stress drop, induced by the bimaterial effect, leads to strong dynamic weakening of the fault and may destabilize the creeping region on a heterogeneous rate and state fault, resulting in extended rupture propagation. Such rupture penetration into creeping patches may lead to more frequent opening of earthquake gates, causing increased seismic hazard. Furthermore, bimaterial coupling may lead to irregular seismicity pattern in terms of event length, peak slip rates,and hypocenter location, depending on the properties of the creeping patches bordering the seismogenically active part of the fault . Our results highlight robust characteristics of bimaterial interfaces that persist over long sequence of events and suggest the need for further exploration of the role of material contrast in earthquake physics and models of seismic hazard.


Introduction
Contrast in elastic material properties across fault surfaces is not uncommon. Examples include strikeslip faults across different rock formations with variable material properties in contact with one another [1], crustal fault zones that have accumulated more damage on one side of the principal slip surface than the other one [2], subduction zones joining continental and oceanic crustal blocks, and glaciers interfacing ice 5 and underlying bedrock [3,4]. Existence of such bimaterial interfaces has been confirmed by seismic imaging studies based on body waves [5,6,7,8] as well as head waves [9,10,11]. The contrast in seismic velocities across fault surfaces may range between less than 5%, as in some strike slip faults, to more than 30%, as in some subduction zones [12]. For example, Fuis et al. (2003) identified a wave speed contrast of 1.09 for the San Andreas fault [6]. For subduction zones, studies for velocity structure in Nankai, show a pressure wave 10 speed contrast of 25%, and a shear wave speed contrast of 45% [13].
Prior theoretical work indicates that mode II ruptures along a frictional bimaterial interface between linear elastic solids produce dynamic changes of normal stress on the fault that depend on the slip rate function, material properties, and direction of rupture propagation [3,14,15,16]. A sub-Rayleigh rupture propagating in the direction of slip on the compliant side of the fault causes a dynamic reduction of normal 15 stress across the fault, whereas a rupture propagating in the opposite direction experiences a dynamic increase of normal stress behind the rupture tip. While the shear stress changes associated with slip grows like the Hilbert transform of the fault-parallel slip gradient, the magnitude of the normal stress changes is proportional to the slip gradient directly. Specifically, normal stress changes increase with the degree of material contrast across the fault, up to the largest contrast for which the generalized Rayleigh wave speed 20 CGR exists, which is the case up to shear wave speed contrast of up to 30-40% for the range of material properties of common rock types.
Experimental work has also been carried out to investigate dynamic rupture propagation along bimaterial interfaces. Anooshehpoor & Brune (1999) performed sliding experiments with two different foam rubber blocks and confirmed that large ruptures propagate in the theoretically predicted preferred direction 25 alongside with the generation of Wrinkle-like slip pulse [17]. Xia et al. (2005) performed sliding experiments along a bimaterial interface for several loading configurations and obtained asymmetric bilateral ruptures [18]. Shlomei et al. (2016Shlomei et al. ( ,2020 studied the onset of frictional instability on bimaterial interfaces and the structure of the slip pulse with a focus on supershear rupture propagation [19,20]. McGuire (2002) and Henry & Das (2001) analyzed rupture properties of over 100 large global earth- 30 quakes and found that most are predominantly unilateral [21,22]. High-resolution locations of numerous small events on the San Andreas fault show a directional asymmetry that is compatible with a preferred propagation direction associated with the local velocity structure [23]. Dor (2006a, b) mapped rock damage at several scales, ranging from 1 to 100 m, in the structures of several faults in southern California including the San Andreas, Punchbowl and San Jacinto faults [24,25]. The results show considerably more damage 35 on the crustal blocks that have faster seismic velocities based on available velocity models [6,7,26,27,28], consistent with the theoretical and numerical predictions of asymmetric damage formation about bimaterial intefraces. However, in a study of 450 small earthquakes at Parkfield, Kane et al.(2013) showed a roughly equal number of events propagating to the SE and NW and only saw greater rupture directivity in the preferred direction (70% SE) if they limited their data sets to those earthquakes of larger magnitudes [29] . 40 The issue of rupture directivity and asymmetry is still a debatable one. Many numerical studies of dynamic ruptures have shown that the presence of bimaterial properties leads to a preferred rupture direction as well as asymmetric rupture features (sometimes propagating unilaterally), a preferred aftershock triggering and asymmetric off-fault damage [30,31,32,33,34,35,36]. Particularly, these numerous studies gave insight into how rupture direction and rupture mode on bimaterial interfaces 45 may be influenced by stress heterogeneity or co-seismic generation of off-fault plasticity. However, these conclusions rely on making apriori assumptions about the background stress and the nucleation procedure. It is thus unclear what part the bimaterial effect is robust over long sequence of seismic and aseismic slip that develop naturally with a smooth nucleation and self-consistent evolving prestress field.
To address this challenge, Erickson and Day (2016) numerically simulated a strike-slip fault governed by 50 rate-and-state friction where quasi-dynamic events nucleate spontaneously due to remote, tectonic loading [37]. They investigated the influence of material contrast over the course of many hundreds of years, and found that the presence of bimaterial properties, with contrasts ranging between 5-20%, influences the earthquake nucleation site, such that ruptures in the preferred direction are more favorable. For smaller values of the characteristic weakening distance in the rate and state friction law, partial ruptures that do 55 not span the whole fault length may emerge. With material contrast present, some of these small events propagate in the non-preferred direction, enabled by a favorable stress distribution left on the fault from previous ruptures. While the work of Erickson and Day (2016) filled an important gap in the literature by being the first work, to the best of our knowledge, to study the bimaterial effect on earthquake cycles in a rate and 60 state framework, it has important limitations. It neglected the effects of wave mediated stress transfer during seismic slip and approximated inertia by a radiation damping term limiting the investigation to the quasidynamic approximation [37]. While prior work on earthquake cycle simulation has highlighted significant differences between quasidynamic and dynamic models [38], including effects on peak slip rate, inter-event time, and event size distribution, inertia effects may be especially important for bimaterial interfaces as they directly influence the magnitude of the co-seismic normal stress changes which is central to the mechanics of rupture propagation.
Here we simulate sequence of earthquakes and aseismic slip on a strike slip fault governed by rate and state friction and separating two half spaces with different elastic properties using a recently developed methodology that resolves all stages of the earthquake cycle and incorporate full inertia effects during 70 seismic rupture. We consider contrasts in the shear wave speeds of the two materials that range between 10% and 30%. The remainder of the paper is organized as follows Section 2 introduces the governing equations. Section 3 outlines the model setup and briefly introduces the numerical method considered here. In Section 4 we evaluate the role of bimaterial interfaces in altering the earthquake cycle pattern, including both aseismic and coseismic contributions, and role of varying the material contrasts. Finally, in Section 5, 75 we conclude by discussing the implications of our results for source physics and seismic observations.

Governing Equations
We consider the two-dimensional domain Ω undergoing plane strain deformations. The domain is divided into two half spaces by a planar strike-slip fault interface S f . In the absence of body forces, the equation of motion for the 2-D plane strain problem are given by 80 ρu α,tt = σ αβ,β α, β = 1, 2 in Ω (1) with Dirichlet boundary conditions applied on S u and Neumann boundary conditions are applied on S T where u α is the displacement vector, σ αβ is the stress tensor, ρ is material density, n β is the outward normal vector (in 2-D) from the boundary. The stresses are given by the linear elastic constitutive law where ε ij is the strain tensor, and µ, and λ are the Lamé parameters, which may differ for each half space.

85
Assuming infinitesimal deformations, the strain tensor is given by Here, T ± o,α is the initial traction and △T ± α is the change in traction due to fault slip, on either side of the fault surface. Imposing continuity conditions at the fault surface we obtain the following jump conditions 90 and stress continuity conditions: where, δ is the slip, and ζ is the fault opening. Additionally, to ensure no interpenetration we enforce that ζ ≥ 0. In subsequent sections, we will use V to denote the slipping velocityδ.

95
Here, we adopt a rate-and-state frictional (RSF) formulation [39,40] used to describe friction in tectonic settings. The boundary condition on the fault surface is enforced by equating the fault shear stress to its strength: where the fault strength F is defined in terms of the effective normal stress σ n and the friction coefficient f . In the RSF, the friction coefficient depends on the slip rate V and state θ as: where L is the characteristic slip distance, f o is the reference friction coefficient defined at a slip rate V o . The state evolution is prescribed through the aging law [41], which is commonly used in earthquake cycle simulations [42,43,44,45] and defined as: This results in a steady-state solution of the state variable θ ss = L/V . The corresponding steady-state 105 friction coefficient is given by: Here, the parameter combination a − b > 0 describes a steady state rate-strengthening frictional response and a − b < 0 describes a steady state rate-weakening frictional response which can lead to unstable slip and stick slip sequences.

110
Rupture nucleation and process zone. RSF introduces a length scale for the nucleation size of earthquake that may be determined using an energy balance approach. Ampuero and Rubin (2008) established the following theoretical estimate for the nucleation size L nuc for a frictional crack under slow tectonic loading [46]: where, µ * = 1 1−ν µ for mode II rupture, µ is the shear modulus, and ν is Poisson's ratio. This nucleation size defines the critical wavelength that has to be resolved within the numerical scheme and is valid for a/b > 0.5. In addition to the nucleation size, Dieterich presented another characteristic length scale L b , which is associated with the process zone during the propagation of the rupture when V θ/L >> 1 and scales as b −1 [47]. The quasi-static estimate for process zone L b is given as: It is vital to properly resolve this length scale as it is more stringent than the nucleation zone's length. For dynamic simulations, continuously resolving the process zone becomes a more challenging ordeal as its size scales with the inverse of the Lorentz factor γ L (v r ) = 1 − v 2 r /c 2 s , where, v r is rupture speed, and c s is the shear wave speed. The Lorentzian contraction implies that the process zone will shrink with increasing 125 rupture propagation speed v r [48].
We further note, that for the case of a planar fault bisecting similar material, the effective normal stress remains constant in time σ n = σ o n , thus we can readily estimate the nucleation size and the process zone through expressions (13) and (14) respectively. For the bimaterial case, the induced deformations introduce perturbations in the normal stress, which alters the computation of these length scales. To this end, and for 130 the sake of representation, we normalize our results using L H This nucleation size is based on a homogeneous bulk with shear modulus µ H and initial normal stress σ o n .
Normal stress regularization. To account for rapid variations in normal stress that could occur due to the material mismatch, we utilize a RSF formulation featuring a delayed response of the shear stress according 135 to Prakash-Clifton law. This model fits observed frictional response better than the traditional formulation [49,16]. In this framework, the fault strength is given by the following [50]: The function ξ evolves exponentially with slip to the new value of σ aṡ where L P C is an evolution distance of choice. In our time stepping, we consider this as an additional evolution equation to compute ξ at any given time. In our analysis, we utilized a proportional scaling of the evolution distance relative to the characteristic length of RSF, such that L P C = 25L. This allows for a sufficiently smooth variation in the frictional strength without deviating substantially from the non-regularized version of RSF.

Model Setup
We consider sequence of earthquakes and aseismic slip on a strike slip fault governed by rate-and-state friction in 2-D plan strain condition. The fault separates two dissimilar half-spaces shown in (Figure 1a). On the fault, a potentially seismogenic patch borders regions steadily moving with a prescribed slip rate 150 V pl .
The parameters used in the simulations are listed in Table 1. We consider different material contrasts where c s and c p are the shear and pressure wave speeds, respectively, and the superscripts refer to the two half spaces. Additionally, Figure 1b shows the heterogeneous spatial distribution of friction parameters to create rheological transitions. Specifically, we consider the effect of the 155 frictional properties of the creeping regions on the patterns of seismicity.

Numerical Scheme
In our numerical model we utilize a coupled finite element and boundary element code FEBE to simulate sequence of earthquake and aseismic slip (SEAS) on a fault surface together with wave propagation in the adjacent medium. This approach was initially introduced to simulate spontaneous dynamic rupture 160 propagation for 2-D inplane problems by Ma et al. (2018) [51]. and later extended to SEAS by Abdelmeguid et al. (2019) in anti-plane setting [52].
Spatial discretization. The coupling procedure between the two methods ensure that continuity is enforced at the interface S s between the FE and BE domain through consistent communication of boundary conditions. This domain truncation approach allow for efficient, yet accurate modeling of the rupture propagation and 165 complex fault zone properties. Accordingly, the FE portion of the domain includes all the material and geometrical complexities associated with the fault zone which is properly resolved as demonstrated in Ma and Elbanna [53], and the homogeneous bulk is modelled using a BE approach known as spectral boundary integral method (SBIM) which is efficient and computationally inexpensive [54,55].
Temporal discretization. Instead of simulating the dynamic rupture propagation of a single earthquake that 170 involves artificially over-stressing the fault, we opt for a more self-consistent approach in which the entire sequence of earthquakes and aseismic slip is simulated. During the inter-seismic phase, aseismic slip is accumulated as the fault is stressed gradually through the background plate tectonic loading V pl applied on the edges of the interface. The accumulation of aseismic slip result in stress concentration that spontaneously nucleate an earthquake consistent with the frictional law and material response. Here, we utilize a quasidynamic approximation during periods of aseismic slip and switch to fully dynamic approach during the dynamic rupture period. During the quasi-dynamic periods we utilize an adaptive time marching scheme proposed by Lapusta et al. (2000) [42]. The time step during dynamic rupture periods is chosen to satisfy the CFL condition. This hybrid, quasi-dynamic-fully-dynamic, scheme is computationally efficient as it enables the choice of large time steps during periods of slow tectonic loading without compromising stability. 180 We switch between quasi-dynamic and fully dynamic solvers based on the value of the maximum slip rate. For the problem discussed below we switch from quasi-dynamic scheme to a dynamic scheme based on a threshold V QD = 1 mm/s , and from dynamic to quasi-dynamic based on a threshold V DQ = 0.5 mm/s. To evaluate the role of transition threshold, we estimate the ratio between radiation damping term given as s µ (1) ) and quasi-static shear stress τ qs . Neglecting the inertia effects 185 is justifiable as long as the magnitude of the radiation damping term is relatively small which is ensured by having the ratio R = ηV /τ qs much smaller than unity R << 1. The above thresholds ensure that ratio R < 10 −4 . Furthermore, we have performed numerical tests to confirm that the accuracy of the obtained results are independant of the threshold choice, as long as it is small enough as outlined above.

190
Slip on a planar bimaterial interface is fundamentally different from slip on a planar interface in a homogeneous bulk. In the bimaterial case, stress transfer between the two half-spaces couples local variations in the normal stress with interface slip. This coupling is independent of the frictional properties of the fault. Rather, it depends entirely on the elastic properties of the two half-spaces, namely, the material mismatch. Furthermore, the bimaterial coupling intensifies near the rupture tip. The sign of the coupling (i.e. whether it 195 introduces compressive or tensile normal stress perturbation) depends on the rupture propagation direction with respect to the sense of motion in the more compliant medium. Here, we adopt the notion that a preferred direction refers to the rupture propagation in the direction of the more compliant material. Along the preferred direction, normal stress decreases behind the rupture tip. In contrast, for the non-preferred direction, where the rupture propagation is in the direction of the stiffer medium particle motion, the normal 200 stress increases behind the rupture tip.

Earthquake cycle on bimaterial interface
Here, we focus on the role of material mismatch on the overall behaviour of the earthquake cycle, including its influence on rupture directivity, the complexity of earthquake sequence, and destabilization of stable creeping regions. We fix the frictional parameters of the rate weakening patch as we vary the shear 205 wave speed c s and pressure wave speed c p of one half-space relative to the other one. Fig. 2 shows the time history for slip rate comparing a homogeneous model (case (1)) to a bimaterial case with r c = 0.7 (case (2)). Figure 2a shows the time history (time steps) of the peak slip rate for case (1) versus case (2). The introduction of material mismatch contributes to higher peak slip rates per event. As the rupture propagates in the preferred direction, tensile change in the normal stress introduces a dynamic 210 weakening mechanism that favors the rupture propagation in that direction resulting in higher slip rates than the homogeneous case. Figure 2b-c illustrates the slip rate contours for both aseismic and co-seismic periods of the earthquake cycle. For the homogeneous case, shown in Figure 2b the problem is symmetric about the center point of the fault x = 0. However, for the bimaterial interface, the coupling between normal stress and fault slip result in asymmetric stress distribution during inter-seismic creep which dictates the 215 nucleation site of the rupture. Specifically, the nucleation site for the bimaterial case is shifted towards the right end of the fault. We observe for the homogeneous case ( Figure 2b) the pattern of earthquakes is periodic with identical seismic events and that the nucleation site remains the same throughout the whole time series.
In contrast, the introduction of material mismatch as shown Figure 2c perturb the earthquake cycle 220 periodicity and lead to several quantitative differences including: 1. Nucleation: For the bimaterial case, ruptures preferably nucleate on the right half of the fault and propagate towards the left in the preferable direction. There is some variability in the hypocenter locations due to the stress heterogeneity that evolve from one event to another.  bimaterial interface yields higher slip accumulation (highlighted in blue) when compared to the homogeneous case (highlighted in red) due to directivety of rupture propagation. In contrast, aseismic slip accumulation on a bimaterial interface is negligible at point x/L H nuc = −1.5 in which the slip rate drops substantially when compared to the homogeneous case at the same location as shown in Figure 3b. The substantial deceleration is mainly attributed to the higher shear stress drop carried by the rupture front propagating in 245 the preferred direction. Figure 3c shows that indeed following the dynamic rupture, the shear stress on the bimaterial interface is substantially lower leading to a near locking of the fault surface. Furthermore, the slip rates for a bimaterial interface are inherently larger as demonstrated in Figure 3b due to the enhanced dynamic weakening effect arising from the coupling between the changes in normal stress and fault slip in the preferred direction 250 In the non-preferred direction x/L H nuc = 1.5 we observe similarities between the bimaterial and the homogeneous cases in terms of slip accumulation during both coseismic and aseismic periods as shown in Figure 3d. Interestingly we observe small sudden accumulation (highlighted in blue and in the zoomed in figure) of slip during aseismic periods. Figure 3e-f elaborates further on this behavior as we observe aseismic accelerations V ∼ 1e −8 (highlighted in blue) and stress drops prior to the large event occurrences that are 255 absent from the homogeneous case. We note that the time between the previous event and this small scale instability is T c = 0.672 years, which is very similar to the recurrence time for the earthquake events on the homogeneous interface. This suggest that the occurrences of these instabilities is tied with an unsuccessful attempt for nucleation on the creeping front. Specifically, the local stress state may initiate the nucleation process, however the average stress state on the fault, resulting from prior events, may cause the arrest of 260 the emerging instabilities. This highlights the important role of emerging heterogeneous stress distribution in regulating slip instabilities. Figure 4a-b shows the slip accumulation during a portion of the the earthquake sequence representing the aseismic slip (blue) plotted every 0.1 years and the coseismic slip (red) plotted every 0.001 seconds. For the homogeneous case a periodic structure of events emerge as demonstrated in Figure 4a. Figure 4b for 265 the bimaterial case, highlights two crucial observations: (1) Slip accumulation in the preferred direction of propagation increases with rupture propagation distance. (2) Dynamic slipping of the stable creeping patch in the preferred direction of propagation suggesting a destabilization effect that is absent in the homogeneous case. The dynamic rupture penetrates within the VS region beyond the extent of the rate weakening patch. This penetration is facilitated by the tensile normal stress perturbation emerging from the coupling between 270 slip and normal stress changes in the preferred direction. While the penetration is initially accompanied by a transient slip deceleration localized at the boundary of the velocity strengthening region, the spacing of the slip lines within that region suggests that the dynamic rupture maintains approximately its speed as it penetrates through the VS region until the rupture is eventually arrested.

Dynamic weakening 275
The rupture penetration shown in Figure 2c and 4b into creeping patch is attributed to bimaterial coupling and subsequent dynamic weakening that occurs as the rupture propagates within the preferred direction. To explore this behavior further, Figure 5a shows the evolution of shear stress versus slip for eighth event in the cycle (highlighted in Figure 3b) at two locations x/L H nuc = 1.5, and x/L H nuc = −1.5. The clear signature of the dynamic weakening in the bimaterial case is highlighted, in which the shear 280 stress drops significantly during coseismic rupture due to the bimaterial coupling in the preferred direction. The dynamic weakening behavior further contributes to the lack of aseismic slip accumulation, and longer inter-seismic periods as shown in Figure 3a Similarly, the weakening effect is further emphasized in figure 5b showing evolution of shear stress versus the slip rate. This dynamic strength reduction in the preferred direction of propagation is the result of 285 rapid normal stress changes carried by the rupture tip due to the bimaterial coupling effect. This observed weakening of the fault strength is carried forward into the creeping patch as rupture continues to propagate in the preferred direction, and predominately contributes to the destabilization of velocity strengthening portion of the fault. Rupture arrest occurs when the dynamic weakening effects from the rupture propagation in the preferred direction can no longer sustain the propagation.

Role of material contrast
To evaluate the role of different material contrasts, we simulate two intermediate contrast ratios corresponding to case (1) r c = 0.9 and case (2) r c = 0.8 commonly observed in strike-slip faults. Figure 6a-b shows the slip rate contours evolution for the time period between 6.5 and 10.5 years. We note that the nucleation size decreases with an increase in material contrast. This is because for a bimaterial interface the 295 effective shear modulus of the medium is reduced due to the existence of a more compliant half space. This in turn reduces the nucleation size which is proportional to the effective rigidity. An estimate of the reduction in nucleation size is computed through the modulus M defined in  [56]. The modulus M depends on the material properties of the two halfspaces and is proportional to material contrast r c . The nucleation size then scales like L nuc ∝ M/µ * H L H nuc [57]. For r c = 0.9 the modulus ratio M/µ * H = 0.9, while 300 for r c = 0.8 the modulus ratio M/µ * H = 0.78 Furthermore, Figure 6c compares the peak slip rate for the two sampled events (I), and (II) (highlighted in the figure). The peak slip rate increases from V ∼ 25 (m/s) to V ∼ 34 (m/s) as the contrast in wave speed increases. This behavior may be explained by considering the following estimate for slip rate based on dimensional analysis V ∝ △τ v r /μ, where v r is the rupture speed,△τ is the stress drop, andμ is the 305 effective shear modulus. Both the rupture speed and the stress drop increase as the material contrast is increased.The effective shear modulus decreases as the material contrast is increased. Thus, the combination of these parameters lead to higher peak slip rates. Finally, the recurrence interval between earthquakes also increased from 0.72 to 0.75 years as the material contrast increase. Qualitatively, this maybe be explained as follows. Assuming a constant stress drop △τ and stressing rateτ , the recurrence interval is T c ∝ △τ /τ . 310 Figure 3 shows that the stress drop increases with the introduction of bimaterial interface. Similarly, the stressing rate drops due to the reduction in effective rigidity. Consequently, the inter-event time increases with increasing material contrast The drop in the normal stress in the preferred direction of propagation, shown in Figure 6d for the sampled events (I), and (II) (highlighted in Figure 6e), increase with material contrast as well as propagation distance 315 from the nucleation site. This increased reduction in normal stress contributes to the slightly extended rupture propagation distance seen in Figure 6a-b. We expect this effect to be amplified on longer faults where the rupture has the opportunity to propagate for larger distances and, hence, experience stronger reduction in the normal stress. 320 We explore the dependence of the normal stress perturbation on the slip rate. Figure 7 illustrates the change in normal stress drop due to changes in the peak slip rate , we observe that at a given material contrast the normal stress drop increases linearly with the peak slip rate △ξ ∝ C ′ V . The slope C ′ depends on the material contrast of the model (as shown in Figure 7). This observation is consistent with the theoretical results for a steadily propagating dislocation [3,58], in which the normal traction is also linearly 325 proportional to the slip rate σ n =μV /v r , where,μ is an algebraic function of the material properties and rupture velocity v r , (Weertman (1980)) which plays the role of an effective elastic modulus. This observation further suggest that the bimaterial effect may dominate the response of long enough faults since slip rate generally increases with propagation distance, at least until the rupture saturates the seismogenic depth.

Rupture penetration into the creeping region 330
In the presence of a bimaterial interface the interaction between normal stress perturbations and the frictional properties of the velocity strengthening patch alters the patch resistance to slip, and may consequently influence the earthquake sequence. Mainly, as mentioned earlier, the coupling between slip and normal stress perturbations could have a destabilization impact enabling rupture penetration within the velocity strengthening patch. This increases seismic risk, by enabling the rupture to propagate longer, 335 and also affect subsequent event dynamics, by influencing the post-event stress distribution. To study this interaction we focus on a specific material contrast r c = 0.9 and vary the frictional properties of the VS region. In addition to the case (1) (a V S − b) = 0.015 presented in the previous section, we consider case (2) where (a V S − b) = 0.01, and case (3) (a V S − b) = 0.005 the choice of parameters indicate a reduction in the fault strengthening behavior within the creeping region since △f ss = (a − b) ln(V /V o ). For case (2) 340 with (a V S − b) = 0.01 Fig.8a demonstrate the extent of the rupture propagation marked by the maximum slip rate contours along the bimaterial interface. The results suggest a larger rupture length due to slip penetration beyond the rate weakening patch and into the rate strengthening region, whose boundary is marked in the figure by dotted black lines, compared to case (1) shown in figure 5a.
Furthermore, in case (2) (a V S − b) = 0.01 (Fig.8a) the earthquake sequence pattern varies significantly 345 by altering the frictional properties of the velocity strengthening patch. Particularly, we observe that the earthquake sequence is not entirely uniform, rather we obtain complexity in the distribution of the nucleation sites (represented by the green stars). Interestingly, for this particular choice of parameters the rupture may nucleate within the left half of the fault indicating that weakening of the velocity strengthening zone could alter the stress state on the fault enabling occasional rupture nucleation in non-preferred sites.

350
The earthquake cycle, however, converge to a pattern of four preferred directional ruptures and one rupture propagating predominately in the non-preferred direction. For this particular a V S − b value, the counterpart homogeneous case events remain symmetric about the center of the fault. Fig.8b shows the rupture extent for case (3) (a V S − b) = 0.005 we observe that the rupture extent is increased further with deeper penetration in the preferred rupture direction. Unlike case (2), we do not 355 observe nucleation of events that predominately propagate in the non-preferred direction. The slip contours indicate that peak slip rate in the preferred direction stays quite large V > 15 m/s as the front penetrates deeper within the creeping region prior to the eventual rupture arrest.The rupture front propagating in the non-preferred direction, however, maintains a lower slip rate due to the increased fault resistance in that direction.

360
To explore the rupture penetration further, Figure 9a shows the extent of rupture penetration L pen for various cases of (a V S − b) and material contrasts. For a homogeneous case we observe that the rupture may also penetrate into the velocity strengthening patch as (a V S − b) decreases due to the reduction in patch resistance. The inclusion of a bimaterial contrast enhances this penetration further as observed for the r c = 0.9 and r c = 0.8 cases. This is primarily attributed to higher reduction in normal stress in the 365 direction of preferred propagation due to the increase in material contrast (as shown in Figure 7a) . For some combination of parameters, the rupture may penetrate through the entirety of the velocity strengthening patch. These observations imply that rupture penetration could be enhanced by either a reduction in the velocity strengthening frictional coefficients or an increase in the normal stress drop. The latter may be facilitated by an increase in material contrast or an increase in the propagation distance in the preferred 370 direction (as shown in Figure 7 and Figure 6d) This correlation may be deduced through estimating the resistance C of the VS patch to slipping, which was described by Kaneko et al. (2010) [59], and defined as: where, σ V S is the effective normal stress in the velocity strengthening patch, D V S is length of the VS patch, V dyn , and V i are the sesimic and interseismic slip velocities on the VS patch. The parameter λ * arises from 375 the short lived stress increases at the rupture tip, and approaches one for small critical slip distance L. Thus, the approximate resistance of the VS patch to slipping depends on the normal stress distribution, as well as, the frictional coefficients of the VS patch for a fixed D V S and similar seismic slip velocities. For a VS patch of a given length, the resistance to penetration can be attributed to either normal stress perturbations, or frictional coefficient changes or both. Utilizing the definition of the patch resistance we 380 can then normalize the minimum normal stress during an earthquake event at the boundary of the velocity strengthening region by the initial normal stress. Furthermore, we introduce a non-dimensional parameter β V S , defined as β V S = (a − b)/(a − b) max on the VS patch, we note that the choice of (a − b) max is arbitrary and made such that the base case considered in our analysis showing no penetration, for the homogeneous case,corresponds to ξ m β V S /σ o n = 1. This non-dimensional quantity ξ m β V S /σ o n couples the effects of both 385 normal stress perturbations and changes in the frictional parameters on rupture penetration. Figure 9b shows that the rupture penetration is indeed correlated with the normalized normal stress for all the simulated earthquakes with variable normal stress changes and frictional properties. The results collapse nearly on a master curve shown in Figure 8b. At lower values of ξ m β V S /σ o n we observe that rupture penetration could grow unbounded within the VS patch. This suggests that for ruptures propagating on longer faults, the 390 larger normal stress drops may cause destabilization of the entire length of the creeping patch.

Discussion
In this work, we consider a planar fault with a velocity weakening patch (VW) bordered by two velocity strengthening patches and separating two half spaces with dissimilar seismic velocities.. We simulate sequences of earthquakes and aseismic slip in which we consider the full inertial and wave-mediated stress 395 transfer during the dynamic rupture portion of the simulation. This enabled us to explore the effect of larger variations in the normal stress arising from the bimaterial coupling beyond what had been investigated earlier. Our results indicate that should the initial conditions on the fault be symmetric, the bimaterial coupling between slip and normal stress introduces heterogeneous stress conditions along the fault during aseismic slip. The heterogeneous stress distribution consequently results in a more favorable nucleation site 400 and asymmetric rupture propagation.
The aseismic loading of the bimaterial interface introduces normal stress reduction accompanying the creeping front propagating in the theoretically preferred direction favoring rupture nucleation on one side (x > 0) of the fault. Accordingly, the rupture nucleation promotes propagation predominately in a preferred direction along the bimaterial interface. This leads to asymmetric rupture propagation, with higher slip rates, 405 on one side of the fault, for the bimaterial versus homogeneous cases. This aligns with experimental results that show asymmetric rupture propagation along bimaterial interfaces [18]. Furthermore, we occasionally observe rupture nucleation in a non-favorable site, which propagates predominantly in a non-preferred direction, leading to smaller peak slip rates and lower slip accumulation. This implies that introducing material mismatch does not always guarantee preferred nucleation. Rather, other parameters, such as stress 410 heterogeneity and, interestingly, frictional characteristics of the velocity strengthening patches, may influence the nucleation. Similar observations were reported in earlier work [60,61,29,23]. A distinct feature of our model is that stresses and seismicity co-evolve naturally over long time scales without imposing arbitrary stress distribution or forced nucleation as routinely done in simulations of single dynamic ruptures.
Through studying different material contrast we observed that rupture asymmetry is enhanced as the 415 material contrast increase. This is mainly attributed to the increase in the normal stress drop which accompanies higher material contrast. Furthermore, changes to the nucleation site with increased material contrast favors longer rupture propagation in the preferred direction. Ruptures propagating in the preferred direction are accompanied by a tensile stress change, and thus achieve dynamic weakening that promotes higher slip rates. These observations are consistent with previous studies for dynamic rupture on a bimaterial 420 interface.
As the dynamic weakening introduced by the bimaterial coupling material increases with material contrast, we observe extended rupture propagation for similar frictional properties. This implies that the elastodynamic coupling introduced by the mismatch in material properties alters the frictional stability of the interface without introducing any changes to frictional parameters. Accordingly, we extended our in-425 vestigation to incorporate the impact of different rate strengthening properties of the creeping patches that border the rate weakening portion of the fault. Our analysis revealed several interesting characteristics. For example, we demonstrate that while the properties of the velocity weakening patch remain similar, and for the same material contrast, we may achieve a higher degree of asymmetry by altering (a V S − b) distribution. A particularly important observation in this study is the interaction between dynamic weakening 430 associated with normal stress perturbations in the preferred propagation direction and the stability of the fault creeping segment. Since the resistance of the VS patch to unstable sliding is correlated to the normal stress on the patch, the extent of rupture penetration is similarly proportional to normal stress at the onset of penetration. Such penetration is particularly hazardous for cases with alternating locked and creeping segments and may lead to more frequent opening of earthquake gates. Longer faults are more susceptible In all our simulations with the current framework utilizing a rate-and-state frictional law and the given fault dimension, we have only observed asymmetric crack-like rupture propagation. It remains to be investigated in future work whether other factors such as strong rate-weakening response or even larger normal stress drops could lead to the emergence of pulse-like ruptures. Furthermore, the nucleation mechanism for ruptures propagating on bimaterial interfaces plays an important role [62]. Prior studies that observed pulse-like transitions for single dynamic rupture event relied primarily on artificial nucleation of rupture [31,15], however, nucleation in the current framework was driven by the aseismic creeping on the fault surface which could play a role in the generation of pulses. Erickson et al. (2016) utilizing a rate-and-state frictional framework observed no pulse generation, however, the study was limited to quasi-dynamic analysis with small normal stress perturbations [37].
In the above study, we focused primarily on the role of frictional properties within the velocity strengthening patch while keeping other parameters fixed, it would be important to consider the interplay between frictional properties within the velocity weakening patch as well. Furthermore, we have limited our investigation to modeling the sequence of earthquakes and aseismic slip in linearly elastic heterogeneous domains 450 undergoing in-plane deformations. We recognize that the off-fault yielding may reduce the asymmetric behavior observed due to bimaterial interface as shown by [34,36]. This will be a focus of future work. and the National Science Foundation CAREER award No. 1753249 for modeling complex fault zone structures. This material is based upon work supported by the Department of Energy under Award Number DE-FE0031685.

Availability of data and materials
The code and input files used in the study are available upon request.