Strain and Stress Accumulation in Viscoelastic Splay Fault and Subducting Oceanic Crust

This study constructs a model to represent the strain–stress field surrounding a splay fault and demonstrates the accumulation processes of both strain and stress around the splay fault. We consider the case of multiple fault dislocations, construct the model as a function of the effective viscosity in the media, and investigate the influence of the effective viscosity on the strain and stress accumulation patterns. The results show that the strain and stress tend to accumulate in the splay fault and the subducting oceanic crust, and the rate of accumulation varies with the effective viscosity. The accumulation and relaxation of strain and stress are simultaneous, and the slower the effective viscosity, the slower the accumulation rate. We discuss the relationship between the splay faults, fluid, and intraslab earthquakes. Finally, the possibility that effective viscosity may contribute to the mode of occurrence of intraslab earthquakes at the Hikurangi subduction margin is discussed.

. (a) Regional tectonic map of New Zealand. (b) Drilled sites in International Ocean Discovery Program . (c) 05CM-04 seismic reflection image (Barker et al., 2018;Wallace et al., 2019). (d) Microseismic activity (Yarce et al., 2019) around the cross section of A′-A. Hypocenters located in a 15-km-width range on either side of the vertical plane defined by A′-A are plotted. The red dashed line indicates the plate boundary estimated by Williams et al. (2013). The entire computational domain (area enclosed by the green line), assumed splay fault (red line), and the primary target of the study (area enclosed by the dashed line) is layered. 10.1029/2023GL103496 3 of 9 as a function of temperature and investigate a model in which splay faults develop and dislocate on geological timescales. Conin et al. (2012) construct a model of an accretionary wedge containing a single splay fault as a function of friction coefficient, and show that the difference between active and inactive slip can be explained by the strength of friction between the plate boundary and the splay fault.
However, no previous studies have accurately considered for the influence of the rheological properties of bi-or tri-materials, such as subducting oceanic crust and foot wall and hanging wall of splay faults present in continental crust. As described in Fagereng et al. (2019) and Morgan et al. (2022), it is clear that the physical properties of the hanging wall and foot wall sides of the splay fault are different. To consider the temporal evolution of strain and stress and their relaxation process, especially on the differences in effective viscosity (bi-or tri-viscoelastic material effect), viscoelasticity should be assumed (Fukahata & Matsu'ura, 2006;Sun et al., 2014). This is an important assumption in the construction of the models that contribute to earthquake risk assessment. This is because it is possible that the stresses corresponding to the stress drop during an earthquake can be relaxed by viscoelastic response during several hundred years (Mitsui & Hirahara, 2001;Pollitz, 2009).
It is also necessary to construct at least a tri-material model to consider the interaction between plate boundaries and splay faults. Within the continental crust, including the splay fault, the physical properties are likely to be heterogeneous between hanging wall and foot wall due to deformation on geologic time scales (Conin et al., 2012;Morgan et al., 2022). Therefore, it is necessary to construct tri-material models with different physical properties.
Here, we propose a strain-stress accumulation model under a tri-viscoelastic material to clarify the influences on the effective viscosity of each medium. In particular, we investigate the effects of viscosity variations on the temporal evolution of strain and stress. Finite element modeling (FEM) is advantageous for numerical and experimental studies of fault interactions because it provides flexibility in dealing with complex fault configurations. We perform a parameter study to determine the magnitude of the effect of effective viscosity. Accordingly, we discuss whether the viscous effect should be considered in the evaluation of seismogenic potential. Finally, the distribution of intraslab earthquakes observed in the Hikurangi subduction margin are compared with the results obtained by the proposed model and we discuss the relationship between splay faults and intraslab earthquakes.

Model Space
We use the (Finite element) software PyLith (Aagaard et al., 2022) to construct a 2D viscoelastic model that includes the subducting oceanic crust, a splay fault and a splay fault junction (the meeting point of the splay fault and the plate boundary), where predefined faults are embedded in the finite elements of the two-dimensional (2D) plane ( Figure S1 in Supporting Information S1). We assumed a restricted rectangular region of 45 km × 10 km as shown in Figure 1d. In this study, we focus on the area enclosed by the dashed line within the entire computational domain and investigate the output in the area enclosed by the dashed line ( Figure 1d). Dirichlet boundary conditions are set at the right, left, and deepest boundaries of the computational domain to set the strain and stress to zero, while the top where the splay fault and the boundary of the computational domain are tangent is set as a free surface. Figure 2 shows the distinction of the three viscoelastic media (a-c) in the model and the associated fault segments (a-c). Fault (b) indicates the splay fault. Five cases are calculated, combining the cases where each of the faults of the media divided into three parts (media) is dislocated and the case where it is not dislocated. Case 1 is when the shallow part of the plate boundary is locked and the interseismic dislocation is partially compensated by the splay fault; cases 2 and 4 are when the shallow or deep part of the plate boundary is locked and the dislocation is not compensated by the splay fault; case 3 is when only the splay fault is dislocated due to external factors not near the junction; case 5 is when both the deep and shallow plate boundaries and the splay fault are simultaneously dislocated. We performed simulations assuming interseismic dislocation from the strainstress state, which is reset once after the occurrence of an interplate earthquake. These cases are based on the observation that dislocations at the plate interface occur partially or entirely during an interseismic period as SSE (Wallace et al., 2012), and that splay faults also dislocate at rates of comparable order to that of plate convergence (Tsuji et al., 2014). The x-component of the strain and stress is rotated to be with parallel to the plate boundary, and the y-component is normal to the plate boundary (left-hand coordinate system). We refine the triangular mesh near the intersection of the subducting oceanic crust and the splay fault ( Figure S1 in Supporting Information S1).
To simplify the problem and investigate the effects of the effective viscosity and viscoelasticity contrast, the Pand S-wave velocities of the structures are set to 2.5 and 0.8 km/s for media (a) and (b), and 5.3 and 3.0 km/s for 4 of 9 medium (c) (the subducting oceanic crust), respectively, and are assumed to be uniform across the media. At the beginning of the calculation, the initial strain and stress are set to zero, and the density is set to 2.5 g/cm 3 for media (a) and (b), 3.0 g/cm 3 for medium (c) (the subducting oceanic crust). The shear and bulk moduli are calculated from the velocity structure, density, and Poisson's ratio. Ruh et al. (2016) and Leah and Fagereng (2022) predict the effective viscosity in the shallow part of the subduction margin to be 1.5 × 10 19 -1.5 × 10 23 Pa⋅ s. Accordingly, the present calculations are performed by varying the effective viscosity in this range of 1.5 × 10 19 -1.5 × 10 21 Pa⋅ s. The effective viscosities of media (a) and (b) are kept constant at 1.5 × 10 21 Pa⋅ s and 1.5 × 10 19 Pa⋅ s, respectively, and the viscosity of medium (c) is varied. The friction coefficient is not set at the faults, and the faults are displaced by forced displacement, corresponding to solving a quasistatic problem. We consider the time interval from 0 to 1,000 years, and the time-step interval is set to 1 year. Fault (b) is displaced at 16 mm/yr and faults (a) and (c) are displaced at 40 mm/yr (Wallace et al., 2012). Physical properties are assumed to be constant at each time step.

Maxwell Model
The Maxwell model is used to determine the viscoelasticity of each medium (Aagaard et al., 2022). The model splits the stress into volumetric and deviatoric components, where the volumetric component is identical to that of an isotropic elastic material, and the deviatoric component is expressed as where τ tot is the total shear modulus of the model and τ 0 is the fraction of the shear modulus accommodated by the elastic spring parallel to the Maxwell models; τ m denotes the fraction of the shear modulus accommodated by each Maxwell model; dev and dev represent the deviatoric stress and strain, respectively. The viscous strain is expressed as follows: where t n is the time between t = 0 and t. ∆ℎ is the viscous strain between t = t n and t. t m is the Maxwell time.

Results
The temporal evolution of the shear strain (ϵ xy ) calculated for cases 1-5 is shown in Figure 2, where a characteristic butterfly-shaped pattern is identified in cases 1, 2, and 4. This simulation confirms the temporal evolution of the shear strain at the edge. Compared to medium (a), the shear strain tends to accumulate in media (b) and (c). As observed, the shear strain tends to accumulate in medium (c) in cases 1, 2, 4, and 5. Comparing cases 1 and 2, the shear strain accumulation is less in case 1 in the vicinity of fault (a) (∼600 years). This directly indicates that the dislocation of fault (b), in addition to fault (a), eliminates a certain amount of the accumulated shear strain near the edge. In case 5, a large shear strain is not to be accumulated at the junction compared to cases 1, 2, and 4 (see red and yellow patterns). The results of the simulations assuming elasticity in all media are shown in Figure  S2 in Supporting Information S1. The results for the continental crust with splay fault are consistent with those presented in Conin et al. (2012). On the other hand, when all media are viscoelastic, shear strain tends to accumulate more easily than when all media are elastic, especially in the subducting oceanic crust.
The temporal evolution of the differential stress in the profiles near the fault (A-D and A′-D′) is shown in Figure 3, which indicates a faster accumulation rate for cases closer to the junction. However, the mode of accumulation varies with the profile. For example, in the case of A-D, the differential stress near the junction accumulates at a linear rate as a function of time compared to another line. Conversely, in the case of A′-D′, the rate of accumulation gradually decreases, which is prominently different from that of another profile. In the case of low effective viscosity, the differential stress accumulates at a slow rate. The same trend is observed in cases 2-4 ( Figures S3-S6 in Supporting Information S1). This is because relaxation proceeds with the accumulation of strain and stress. The differential stress accumulated over 1,000 years differs by ∼50 MPa between the cases with an effective viscosity of 1.5 × 10 19 Pa⋅ s and 1.5 × 10 21 Pa⋅ s.

Dislocation Patterns
The process of strain and stress accumulation varies with the dislocation pattern. Cases 1-5 simulate the possible dislocation patterns that occur near the splay fault and plate interfaces. Generally, the strain and stress tend to accumulate at the edge of the dislocation area, which has been confirmed in both dynamic fracture simulations (DeDontney et al., 2012) and quasistatic simulations (Conin et al., 2012), similar to the present results. Some results are similar when comparing the assumption of a viscoelastic medium with the assumption of an elastic medium ( Figure S2 in Supporting Information S1). In detail, the shear strain tends to accumulate around the foot wall side of the splay fault and the subducting oceanic crust. The subducting oceanic crust tends to accumulate more strain when a viscoelastic medium is assumed.
The strain accumulated in the splay fault (fault (b)) is partly because of the differences in the dislocation rates of fault (a), fault (c) and the splay fault (fault (b)), with compression and tension occurring near the junction (Figure 2). Since the dislocation is treated as a forced displacement in our model, the balance of bulk forces such as the critical taper model could not be considered. However, it is obvious that the spatial pattern of the accumulated shear strain is partly due to the variations in the dislocation rates between the plate interface and splay fault. Considering our results and those of Conin et al. (2012), the combination of the dislocating fault planes, dislocation rate, force balance, fault plane conditions and effective viscosity can form the spatial shear strain patterns.

Effect of Effective Viscosity on Temporal Evolution of Strain and Stress
The influence of effective viscosity increases with time, with differences in the magnitude of accumulated strain and stress depending on the effective viscosity (Figure 3). For the time scale in this study, when the effective viscosity is 1.5 × 10 21 Pa⋅ s, the results are almost the same as when the medium is assumed to be elastic. On a longer time scale, even when the effective viscosity is greater than 1.5 × 10 21 Pa⋅ s, the difference in the accumulation rate of the differential stress is significant depending on whether the medium is assumed to be elastic or viscoelastic. While this argument is generally true in all cases, it is difficult to discuss the case where the differential stress becomes a sigmoidal curve as a function of time, as in case 2. This behavior depends on the fault configuration and the location of the output relative to that of the fault edge. The accumulation rate of the differential stress could be faster, depending on the adjacent media type or physical properties (Hirano & Yamashita, 2011). This suggests that when the medium is viscoelastic and the effective viscosity is heterogeneous, the seismogenic potential may increase locally. A more detailed study is needed to verify this possibility.
In the case of steady-state strain, such strain accumulation and relaxation are compatible. In other words, relaxation occurs simultaneously with the accumulation of strain and stress due to the viscoelastic response to stepwise forces such as dislocation (Fukahata & Matsu'ura, 2016; see Supporting Information S1). Because viscoelasticity depends on the time constant of the external force with its spatiotemporal frequency dependence (Sumita & Manga, 2008;Sun et al., 2014), our simulation differs from the relaxation for pulsed dislocations such as after slip. In addition, the effective viscosity of the medium varies depending on the geology. Therefore, there are viscous heterogeneous and spatiotemporal heterogeneous relaxation phenomena on very small spatial scales in a natural environment (Barnes et al., 2020;Leah & Fagereng, 2022). The relaxation phenomena with different time constants (Maxwell time) are expected to interact with each other. If we will further discuss the interaction of the relaxation phenomena, we should consider the ratio of effective viscosity between each structure. In other words, it will be necessary to set up a bi-or tri-material problem in which multiple patterns are tested for the effective viscosity ratio between each structure, rather than just changing the effective viscosity of one structure.
The variation of the accumulated differential stress depends on the effective viscosity of the medium. Moreover, the accumulated differential stress is in the same order of magnitude or greater than the amount of stress drop of earthquakes (Figure 3). Our results show quantitatively that the effective viscosity in the medium significantly affects the strain and stress accumulation process. However, the actual evolution of the relaxation in a natural environment would not be easy to determine. More potentially, the response after the completion of relaxation is crucial. In the future, the evaluation of a multifaceted approach, such as geodetic observations using strain gauges, observation of core microstructure, and velocity structure analysis, will provide a better understanding of the actual relaxation occurring at the subduction margin.

Implications for the Hikurangi Subduction Margin
When all media are viscoelastic, shear strain tends to accumulate more easily than when all media are elastic, especially in the subducting oceanic crust (Figure 2). This suggests that the viscoelastic medium may contribute to the occurrence of intraslab earthquakes. Using the Hikurangi subduction margin as an example, we discuss the relationship between the strain-stress field around the splay fault and the intraslab earthquakes. In this section, we focus on the output in medium (c), which corresponds to a subducting oceanic crust. For example, as shown in Figure 1, numerous microearthquakes have occurred in the subducting oceanic crust at the Hikurangi subduction margin.
In principle, intraslab earthquakes occur because of the ease of stress drop due to constituent materials and the fluid migration caused by episodic SSE in a subduction margin (Kita & Katsumata, 2015;Warren-Smith et al., 2019). In addition to these factors, our result suggests that the shear strain concentrated at the ends of faults may contribute to these occurrences ( Figure 4). The presence of fluid migration as well as the geological diversity of the Hikurangi subduction margin complicates the stress accumulation and relaxation (Barnes et al., 2020). Our results suggest that the stress perturbations at the fault edges and around the splay faults, caused by displacements on the order of hundreds of years, extend into the subducting oceanic crust and vary considerably with the magnitude of its effective viscosity. It should be noted that our model is subjected to many assumptions and simplifications that may not apply in a natural environment. However, our model suggests that the ease or difficulty of occurrence of intraslab earthquakes may be partially controlled by the effective viscosity of the constituent materials, which is indirectly supported by Kita and Katsumata (2015).
Our results for the spatial distribution of the stress accumulation in the subducting oceanic crust at the end of the splay fault are competitive with the stress field found by Warren-Smith et al. (2019) using the local earthquake catalog, which is completely independent of the stress accumulation expected in our model (Figure 4). In most cases, the differential stress accumulation tends to plateau at medium (c). In other words, in a natural environment, stress accumulation on a longer time scale than that considered in this study may tend to plateau in the subducting oceanic crust (Figure 3; Figures S3-S6 in Supporting Information S1). In contrast, if the stress perturbation, or stress change around the fluid filled fracture is caused by the SSE, the spatial pattern of stress would be complicated because the stress caused by the steady-state dislocation would be partially relaxed. In addition, a large number of splay faults are actually present, which could influence different slip modes in the subduction margin (Fagereng et al., 2019). The interaction of these influences is thought to be responsible for the widespread distribution of intraslab earthquakes. Quantitative studies using Coulomb failure stress change for multiple events is a future work. It is expected that the universality of the results of this study will be ensured by more global knowledge, such as computational results using other viscoelastic models.

Conclusions and Perspectives
Quantifying of the accumulation and relaxation of strain and stress around the splay fault during earthquake cycles is critical to understanding the fundamental physics of the subduction system. This study constructs a viscoelastic finite element model to enhance the representation of the strain-stress field around the splay fault, based on which the strain and stress accumulation processes around the splay fault are discussed. Our result suggests that strain and stress tend to accumulate along the foot wall side of the splay fault and the subducting oceanic crust. In a viscoelastic medium, the accumulation and relaxation of strain and stress occur simultaneously, and the accumulation rate varies with the effective viscosity. This suggests that the strain and stress concentrated at the ends of splay faults could contribute to intraslab earthquakes such as those observed at the Hikurangi subduction margin, and that the mode of occurrence depends on the effective viscosity. Although this study focused on interseismic strain and stress, the seismic potential can be further evaluated by simulating different cases and accumulating knowledge. Our model suggests that a more representative model can be constructed by systematically and quantitatively investigating the individual effects of each element.