The role of heterogeneous stress in earthquake cycle models of the Hikurangi–Kermadec subduction zone

SUMMARY

Seismic and tsunami hazard modelling and preparedness are challenged by uncertainties in the earthquake source process. Important parameters such as the recurrence interval of earthquakes of a given magnitude at a particular location, the probability of multifault rupture, earthquake clustering, rupture directivity and slip distribution are often poorly constrained. Physics-based earthquake simulators, such as RSQSim, offer a means of probing uncertainties in these parameters by generating long-term catalogues of earthquake ruptures on a system of known faults. The fault initial stress state in these simulations is typically prescribed as a single uniform value, which can promote characteristic earthquake behaviours and reduce variability in modelled events. Here, we test the role of spatial heterogeneity in the distribution of the initial stresses and frictional properties on earthquake cycle simulations. We focus on the Hikurangi–Kermadec subduction zone, which may produce M_w > 9.0 earthquakes and likely poses a major hazard and risk to Aotearoa New Zealand. We explore RSQSim simulations of Hikurangi-Kermadec subduction earthquake cycles in which we vary the rate and state coefficients (a and b). The results are compared with the magnitude-frequency distribution (MFD) of the instrumental earthquake catalogue and with empirical slip scaling laws from global earthquakes. Our results suggest stress heterogeneity produces more realistic and less characteristic synthetic catalogues, making them particularly well suited for hazard and risk assessment. We further find that the initial stress effects are dominated by the initial effective normal stresses, since the normal stresses evolve more slowly than the shear stresses. A heterogeneous stress model with a constant pore-fluid pressure ratio and a constant state coefficient (b) of 0.003 produces the best fit to MFDs and empirical scaling laws, while the model with variable frictional properties produces the best fit to earthquake depth distribution and empirical scaling laws. This model is our preferred initial stress state and frictional property settings for earthquake modelling of the Hikurangi–Kermadec subduction interface. Introducing heterogeneity of other parameters within RSQSim (e.g. friction coefficient, reference slip rate, characteristic distance, initial state variable, etc.) could further improve the applicability of the synthetic earthquake catalogues to seismic hazard problems and form the focus of future research.

1 INTRODUCTION

Seismic hazard assessment is an important basis for building codes, risk assessment, insurance rate calculation, site selection for infrastructure and the assessment of earthquake-related hazards such as tsunami and shaking-induced landslides (Chan et al. 2020; Kumar et al. 2022; Gerstenberger et al. 2022a). It can also provide information for real-time decision-making and incident response immediately after an earthquake. However, the characterization of where and how earthquakes may occur—including information such as recurrence intervals of earthquakes of a given magnitude at a given location, or the probability of multifault or multiple-segment rupture—can be highly uncertain. It is especially difficult to characterize the hazard from large earthquakes (M > 7.0) because the time period for which modern seismological observations are available is much shorter (∼50–120 yr) than the recurrence interval of most active faults (∼100–10 000 yr). Although palaeoseismic data provide coverage of large earthquakes further back in time, their use in hazard modelling is somewhat complicated by uncertainties in their ages, the dimensions of the fault rupture area and how to extrapolate from individual palaeoseismic sites to the rest of a fault surface (Dieterich & Richards-Dinger 2010). Furthermore, the complexity of the earthquake source process is not fully taken into consideration or, if it is, carries significant uncertainty in seismic hazard assessment (Stirling et al. 2012; Field et al. 2014; Schwartz 2018; Page 2020; Gerstenberger et al. 2022a). The use of uniform source parameters (e.g. uniform or smoothly varying coseismic slip), simple scaling relations (e.g. average slip versu. magnitude) and empirical ground motion models nominally provide the results of an average condition, which might not be suitable for hazards that are sensitive to variability in the earthquake source, such as local-scale ground motions and tsunami inundation.

Physics-based earthquake simulators offer one way to assess uncertainties in hazard-related input parameters (Field 2019). They generate long-term catalogues of synthetic earthquakes based on simple models of fault physics. Previous studies have shown that in California, ground motions modelled using synthetic earthquake catalogues are comparable to those from traditional earthquake-rupture forecasts, that is, earthquake simulators can generate similar earthquake catalogues to the instrumental catalogues (Richards-Dinger & Dieterich 2012; Tullis et al. 2012; Shaw et al. 2018). However, in these studies, uncertainties in several parameters, such as fault geometry and initial stresses, mean that model inputs are often simplified, with relatively smooth fault geometries and spatially uniform fault stresses. Collectively, these simplifications tend to produce ‘characteristic’ earthquakes, commonly recurring events with similar rupture evolution and final slip distribution on a given fault over multiple seismic cycles (Shaw et al. 2022). This characteristic behaviour might mean that these models do not consider a sufficiently wide variety of earthquake sources, affecting hazard assessment results. Delogkos et al. (2023) have shown the effect of the complexity of fault geometry and slip rates on simulated earthquake catalogues using an earthquake simulator. The investigation of the sensitivity of synthetic earthquake catalogues to the variation of initial stress has not been studied. Some researchers have shown the sensitivity of conceptual earthquake cycle simulations to frictional parameters such as rate-and-state coefficients and fault bend (Kaneko et al. 2010; Yang et al. 2012; Noda & Lapusta, 2013; Sathiakumar et al. 2020). In this study, we deviate from these previous studies by focusing on a more complex model of the fault system and varying its initial stress conditions in an effort to (1) understand the sensitivity of large-scale earthquake cycle simulation for realistic fault systems and (2) obtain results that can directly be used in future applied hazard research such as that presented in Milner et al., (2021) and Shaw et al. (2018).

In this study, we focus on the effect of the complexity of stresses on faults by introducing heterogeneity into the distribution of initial normal and shear stresses of the physics-based earthquake simulators. To assess the impact of heterogeneous initial stresses, we compare earthquake catalogues for homogeneous and heterogeneous initial stresses on the Hikurangi–Kermadec subduction interface, which extends north–northeast from the eastern North Island of New Zealand (Fig. 1). This subduction zone has the potential to generate large (M_w > 9.0) subduction earthquakes and is therefore very important to consider in seismic and tsunami hazard assessments for New Zealand and the Pacific (UNESCO/IOC 2020). Despite this, seismic and tsunami hazard assessments are limited by a lack of constraints on the characteristics of large earthquakes on the Kermadec segment, primarily because of the lack of nearby land to record past great earthquakes or deploy monitoring instruments. Thus, we use models of the Hikurangi–Kermadec subduction zone to test the effect of stress heterogeneity on earthquake-cycle behaviour with an earthquake simulator.

2 THE RSQSIM SIMULATOR

We use RSQSim, a physics-based earthquake simulator based on rate-and-state friction laws (Richards-Dinger & Dieterich 2012), to generate long-term earthquake catalogues for the Hikurangi–Kermadec subduction zone. Dieterich (1979) introduced the constitutive law of rate-and-state friction from experimental observations of fault-sliding behaviour,

where ${\tau }^{\mathit{\text{fric}}}$ is the frictional shear stress-resisting motion, $\sigma$ is the effective normal stress, ${\mu }_0$ is the steady-state friction coefficient at the reference slip rate and constant normal stress, V and $\theta$ are the slip rate and state variable, $V_0$ is the reference slip rate, $D_c$ is the characteristic distance and a and b are dimensionless constants for rate and state dependence of friction, respectively. At steady state, $\theta =\frac{D_c}{V}$⁠, the frictional shear stress will be:

If (b − a) > 0, a fault exhibits velocity weakening, implying that friction decreases as sliding velocity increases, while if (b − a) < 0 it exhibits velocity strengthening and stable sliding. Velocity weakening is the behaviour associated with most stick-slip earthquakes, while velocity strengthening is the behaviour associated with creeping on the fault surface. RSQSim was developed by Dieterich & Richards-Dinger (2010) based on the computational scheme of earthquake faulting from Dieterich (1995). In RSQSim, the seismic cycle based on rate-and-state friction on each fault element is simplified to four states, healing (State 0), nucleating (State 1), seismic rupture (State 2) and stable sliding (State 3). Note that stable sliding only applies to fault elements with velocity strengthening behaviour (b − a < 0). Quasi-static elastic interactions of elements and constant slip rate during dynamic sliding are also considered in the program.

3 SIMULATING SEISMICITY ALONG THE HIKURANGI–KERMADEC SUBDUCTION INTERFACE

The Hikurangi–Kermadec subduction interface dips northwest and extends ∼2100 km north–northeast from near Kaikōura, Aotearoa New Zealand, to around 500 km north of Raoul Island (Fig. 1). We adopt the geometry of Hikurangi–Kermadec subduction interface from the 2022 New Zealand National Seismic Hazard Model (NZ NSHM, Van Dissen et al. 2022). Models of interface geometry are provided for the Hikurangi portion of the subduction system by Williams et al. (2013) and for the Kermadec portion by Hayes et al. (2018). In the 2022 NZ NSHM, these were blended into a single Hikurangi–Kermadec geometry. This geometry was meshed into triangles with a side length of ∼10 km in this study.

We impose stressing rates using backslip loading and the slip-deficit model developed for the combined Hikurangi–Kermadec interface for NZ NSHM. This slip-deficit model (shown in Fig. S1, Supporting Information) is primarily constrained by geodetic data in the Hikurangi portion of the study area. Further north, slip-deficit rates in the model are poorly constrained, but are based on modelled relative Australia-Pacific plate motions and an approximate coupling distribution inferred from seismicity rates and geodetic observations of deformation at a single site at Raoul Island. Like the 2022 NZ NSHM, we treat this slip-rate deficit model as a representation of the total slip rate that is released coseismically in earthquakes (Gerstenberger et al. 2022a;b). We therefore use slip rate deficit as the direct input for calculating stressing rates, and do not explicitly include any contribution to subduction interface slip from aseismic creep or slow-slip events. Note that since the degree of locking close to the trench in the Hikurangi subduction zone is unknown, Van Dissen et al. (2022) created two alternative ‘trench-locked’ and ‘trench-creeping’ slip-deficit models. We use the trench-creeping model in our study, but do not expect this model selection to impact our results significantly because the Hikurangi margin represents a small, slow-slipping portion of the entire modelled subduction interface. Instead, the high stressing rates north of Aotearoa New Zealand (which are identical between the two alternative portions) account for most of the M_w > 7.0 events in our models.

4 MODELS

4.1 Computation of initial stresses

As described in the following sections, we consider a number of different types of initial stress variations for our models. We vary the initial values for both the effective normal and the shear stresses. When assigning initial stresses, we assume that the fault patches are in a reverse faulting environment and that the maximum principal stress (⁠${\sigma }_1$⁠) is aligned with the local convergence direction. Convergence directions are obtained from the results of Wallace et al. (2012) and Power et al. (2012) for the Hikurangi and Kermadec portions of our models, respectively. The initial stresses required for RSQSim are the effective normal stress (⁠${\sigma }_n^{\prime }$⁠) and the shear stress (⁠$\tau$⁠) on the subduction interface. These stresses are related through a simple static friction law (Sibson 1985):

where $\mu$is the long-term friction coefficient (assumed here to be 0.6) and P is the pore-fluid pressure. Assuming the minimum principal stress (⁠${\sigma }_3$⁠) is equal to the gravitational stress, the effective minimum principal stress (⁠${\sigma }_3^{{}^{\prime }}$⁠) would be:

where $\rho$ is the density of crust (3000 $\mathrm{kg}{\mathrm{m}}^{-3}$⁠), g is gravitational acceleration, z is the depth in the crust and $\lambda$is the pore-fluid pressure ratio.

Following Sibson (1985), we can formulate eq. (3) in terms of the effective principal stresses and obtain the ratio of the maximum and minimum effective principal stresses for fault reactivation,

where ${\sigma }_1^{\prime }$ is the effective maximum principal stress and $\theta$ is the dip angle of the fault. In terms of the effective principal stresses, the initial effective normal (⁠${\sigma }_n^0$⁠) and shear (⁠${\tau }^0$⁠) stresses on fault are:

We also consider the possible depth dependence of fault rheology that may arise due to increasing pressure and temperature (e.g. Blanpied et al. 1991; Scholz 2019), resulting in a brittle–ductile transition. The ductile stress variation is calculated by considering the thermal profile of the Hikurangi subduction margin from Fagereng & Ellis (2009). Since we have no information about the thermal condition of the Kermadec segment, we extend the thermal profile of the Hikurangi margin along strike to the northern end of the Kermadec segment. We assume the viscous shear stress along the fault is controlled by a steady-state dislocation creep flow law in simple shear (Fagereng & Ellis 2009):

where $\dot{\gamma }$ is the shear strain rate, Q is the activation energy, R is the ideal gas constant, T is the temperature in degrees Kelvin, A is a material parameter, $f_{{\mathrm{H}}_2\mathrm{O}}$ is water fugacity, and m and n are power-law exponents for $f_{{\mathrm{H}}_2\mathrm{O}}$ and stress, respectively. We adopt the values of the parameters in eq. (8) from Fagereng & Ellis (2009) and Hirth et al. (2001): $\dot{\gamma }={10}^{-12}{s}^{-1}$⁠, Q = $135\mathrm{kJ}\mathrm{mol}{}^{-1}$⁠, R = 8.31, A = ${10}^{-11}\mathrm{MPa}{}^{-n}{\mathrm{s}}^{-1}$⁠, $f_{{\mathrm{H}}_2\mathrm{O}}$ = 37 MPa, m = 1 and n = 4.

4.2 Models with constant values for rate-and-state coefficients and pore-fluid pressure ratio

We first consider models where we use constant values for the rate-and-state coefficient difference (b − a) and the pore-fluid pressure ratio (λ). Our simplest models (hereafter referred to as ‘homogeneous’) assume that the initial stresses are constant, while a slightly more complex model (hereafter referred to as ‘heterogeneous’) assumes depth-varying gravitational stress with initial stresses as defined by eqs (6) and (7) and a constant pore-fluid pressure ratio. The homogeneous stress model has initial values of uniform effective normal and shear stress on the entire subduction interface of 248 and 149 MPa, respectively (the averages from the heterogeneous stress model, Fig. 2a). In the heterogeneous stress model, the initial stresses vary according to the depth and geometry of the fault patches (Fig. 2b).

For the heterogeneous stress model, we assume that the pore fluid is overpressured and has a constant ratio (⁠$\lambda$⁠) of 0.7 for the subduction interface according to the suggested range of $\lambda$ in Ellis et al. (2015) and Arnulf et al. (2021). We also propose stress models with hydrostatic conditions to understand the effect of the change of pore-fluid pressure in our stress models. The effective initial stress values in hydrostatic models are roughly twice those of the overpressured models at the same depth (see blue and light blue lines in Fig. 2d). Finally, we consider a model where the stresses below the brittle–ductile transition are defined by eq. (8); we refer to this as the ‘brittle–ductile’ model (Fig. 2c). Unlike the heterogeneous model, the shear stress values of the brittle–ductile model decrease to an extremely small value below the brittle–ductile transition (Fig. 2d). The effective normal stresses at depth fall between those calculated for the heterogeneous and homogeneous models. Note that the stresses towards the trench approach zero because of the shallow depth. The stresses are not tapered.

As per eq. (2), the coseismic slip and stress drop in an earthquake are affected significantly by the difference of the rate and state coefficients (b − a). Thus, we test four (b − a) values by fixing a to 0.001 and exploring b = 0.004, 0.003, 0.002 and 0.0015 for all the stress models. In order to understand the effect of the variation of the initial stresses on the earthquake cycles exclusively, we assume that (b − a) is uniform for the entire subduction interface in this set of models.

We adopted the default values from RSQSim for the other parameters in eqs (1) and (2): ${\mu }_0$= 0.6, $D_C$= ${10}^{-5}$ m, $V_0$ = ${10}^{-6}$ m s⁻¹ and initial $\theta$ = $2\times {10}^8$ s. To focus on testing the effect of stress heterogeneity, these parameters were also set as constant along the subduction interface. Only the initial stresses were varied according to the depth and geometry of each patch. Each model run is stopped when either the number of earthquakes reaches 2 000 000 (6 000 000 for the models with b of 0.0015 to reach the desired catalogue length), or the length of the synthetic catalogue reaches 110 000 yr. We assume, conservatively, that the first 50 000 yr of each run are a ‘spin-up’ period and discard them from our analysis. We recognize that the homogeneous and heterogeneous initial stress models will all have heterogeneous stresses following the spin-up period of the simulations. However, for the remainder of this paper, we refer to those models that had initial homogeneous stress as ‘homogeneous’ and those that had initial heterogeneous stress as ‘heterogeneous’.

4.3 Models with variable rate-and-state coefficient and pore-fluid pressure ratio

In the previous section, we assume the rate and state coefficients (a and b) as constant for the entire subduction interface. However, the properties of the material may vary with depth on the subduction interface. The values of a and b in the rate-and-state friction law are determined by the material properties. As we mentioned in Section 2, a fault can exhibit velocity weakening or strengthening when (b − a) is positive and negative, respectively. When (b − a) is positive, the materials are more brittle and unstable forming a seismogenic zone on the fault surface. On the other hand, materials with negative (b − a) are more ductile and stable creating a creeping zone on the fault surface. Scholz (1998) proposed a generic (b − a) model for subduction interface indicating that the materials at the shallow and deep parts of the subduction interface tend to behave as velocity strengthening with (b − a) < 0 because of the unconsolidated sedimentary wedges and the high temperature and pressure at depth, respectively. The middle range of the depth of the subduction interface displays velocity weakening with (b − a) > 0. This variation of (b − a) values with depth can be observed in the distribution of the earthquake depths, that is, more seismicity at the seismogenic zone in the middle of the fault and less seismicity at the top and bottom of the fault. The model proposed in Scholz (1998) only provides the synoptic shape of (b − a) along the depth. To obtain the range of the (b − a) values, we refer to the (b − a) model as a function of temperature from the frictional sliding experiments for illite-rich gouge from den Hartog & Spiers (2013). We incorporate the thermal profile of the Hikurangi subduction interface from Fagereng & Ellis (2009) into the (b − a) model to obtain the relationship between (b − a) and depth. The shallower transition temperature of velocity strengthening and weakening is around 250 °C in the (b − a) model from den Hartog & Spiers (2013), which corresponds to a depth of around 20 km according to the thermal profile of Fagereng & Ellis (2009). This transition depth does not match the depth distribution of the seismicity in the study area (Fig. 10). We then modified the top transition temperature of velocity strengthening and weakening to around 120 °C according to the depth distribution of the seismicity in the study area. Fig. S2 (Supporting Information) shows the modified (b − a) model we use in this study. Due to the poor performance of RSQSim when the (b − a) value is close to zero, we modified the (b − a) model to skip the (b − a) values smaller than 0.001 at depths of around 10 and 40 km. In the construction of rate-and-state models in Arnulf et al. (2021), (b − a) < 0.001 is considered as nearly velocity neutral and is suitable for simulating slow-slip events. Skipping the (b − a) values between −0.001 and 0.001 will possibly eliminate the generation of slow-slip events in the synthetic catalogues. Although slow-slip events on the Hikurangi margin are well-known, the focus of this study is to examine the effects of the variation of initial stress state and frictional properties on the resulting catalogues of RSQSim. Including velocity neutral (b − a) in the synthetic earthquake simulation also requires a lot more computational resource making the modelling more difficult to converge. This elimination of (b − a) values between −0.001 and 0.001 causes two ‘steps’ in the otherwise smooth profile in Fig. S2 (Supporting Information).

Our heterogeneous stress model presents the increase of stress with depth using only one slope because of the constant pore-fluid pressure ratio. Saffer & Tobin (2011) suggested that the effective normal stress increases gradually at shallow depth (within 10 km) based on borehole observations and calculated pore-fluid pressure from P-wave interval velocity for the Nankai Trough. Saffer & Tobin (2011) further proposed a conceptual pore-fluid pressure ratio distribution along the depth for subduction zones including two peaks of pore-fluid pressure ratio—one shallow, one deep—based on geophysical observations of subduction zones from previous studies. We generate a pore-fluid pressure ratio model with peaks at 12 and 40 km according to the modelled pore-fluid pressure ratio profile of the Hikurangi subduction zone from Ellis et al. (2015) and the model from Saffer & Tobin (2011), respectively (Fig. S3a, Supporting Information). Fig. S3b (Supporting Information) shows the stress profile incorporated with the variation of pore-fluid pressure ratio. By adding the peaks of pore-fluid pressure ratio, the effective normal stress increases gradually for depths shallower than 15 km in our model.

We first incorporate the variation of (b − a) and pore-fluid pressure ratio along depth into our stress models separately for the realization of the effect of each parameter. A model with the variation of both (b − a) and pore-fluid pressure ratio was also performed for the comparison of the levels of complexity of parameters. The results of the models mentioned above and the comparison with the models with constant frictional properties are discussed in Section 6.4.

5 RESULTS FOR MODELS WITH CONSTANT FRICTIONAL PROPERTIES

5.1 Magnitude–time and magnitude–frequency distributions

Magnitude-time and magnitude–frequency distributions (MFDs) were developed for the homogeneous, heterogeneous and brittle–ductile stress models described in Section 4.2 (Figs 3 and 4). The magnitudes of the smaller earthquakes (5.4 ≤ M < 6.1) are controlled by the (b − a) values and the sizes of the triangle patches, which are also mostly uniform along the interface (30–55 km²). Therefore, we could only generate events with certain magnitudes corresponding to the rupture of one or two triangles in lower magnitudes. Though there are still events clustered at magnitudes around 5.7–5.8 in the heterogeneous and brittle–ductile models with b of 0.003 and 0.004, the variation in stresses of the heterogeneous and brittle–ductile stress models allows us to produce earthquakes in the gaps between certain magnitudes. The events in the heterogeneous stress models are more evenly distributed by magnitude, especially from magnitude 7 to 8.5 (Fig. 3, centre and right). This feature can also be seen in the MFDs (Fig. 4), where differences between the instrumental and synthetic MFDs of the heterogeneous models are much smaller than the homogeneous models. Both minimum and maximum magnitudes of the synthetic catalogues become smaller when the state coefficient b is decreased. The homogeneous model with b = 0.0015, heterogeneous model with b = 0.0015, and 0.002 and brittle–ductile stress model with b = 0.0015 and 0.002 do not produce any events larger than magnitude 9 (Fig. 3).

To compare with the model MFDs, we compute the MFD of the instrumental earthquake catalogue along the Hikurangi–Kermadec subduction zone, using the earthquakes identified as being on the Hikurangi–Kermadec subduction interface (Rollins et al. 2022). The interface events are defined based on the distance between the events and the subduction interface and the similarity of the focal mechanism to the interface geometry. The period of the filtered catalogue is from 1917 to 2018 with magnitudes filtered from 5 to 8.17, where a single event with magnitude over 8 happened in 1917. To make this comparison quantitative, we estimate the Gutenberg–Richter b-value (Gutenberg & Richter 1944) (b_GR in Fig. 4) of the instrumental catalogue and the negative exponential parts of the models. Considering the instrumental catalogue only contains events for around 100 yr, it is not sufficiently long in duration to constrain the upper magnitude limit and we only compare the slopes of the MFDs. Since the magnitude range varies in different synthetic catalogues, a grid search for determining the magnitude of completeness was applied when the b-values were estimated (Figs S4–S7, Supporting Information). We estimate the b-value of the instrumental catalogue as 1.12 (Fig. 4). The homogeneous stress models have estimated b-values ranging from 1.368 to 2.013 and their strongly characteristic-like MFDs are a poor match to the instrumental catalogue, that is, large deviation from the instrumental catalogue in M7–8.5 (Fig. 4a). The heterogeneous stress models have b-values (1.228–1.672) and MFDs much closer to the instrumental catalogue, that is, the deviation from the instrumental catalogue in M7–8.5 does not exist in the MFDs of the heterogeneous models (Fig. 4b, solid lines). The models with higher state coefficients (b = 0.004 and 0.003) fit the instrumental catalogue better than the models with lower state coefficients (b = 0.002 and 0.0015). The hydrostatic models have similar b-values to the corresponding heterogeneous stress models but produce higher maximum magnitudes (M_Max) than the heterogeneous stress models. The hydrostatic model with state coefficient b of 0.004 has an MFD which is the closest to the seismicity rate model from NZ NSHM in larger magnitudes (M > 8, dashed grey line in Fig. 4). Among the brittle–ductile transition models (Fig. 4c, dashed lines), the b-values (1.351–1.901) are higher than the heterogeneous models and the instrumental catalogue due to the high seismicity rate in lower magnitudes. Overall, adding heterogeneity to the initial stress produces synthetic catalogues that more closely match the instrumental catalogue, particularly if a higher state coefficient is used. We also see that by increasing (b − a), we are able to reduce the deviation from the observed catalogue at the smallest magnitudes.

We also examine whether the seismic moment rates of the synthetic and the instrumental catalogues are identical. The seismic moment rate of the instrumental catalogue is $8.46\times {10}^{19}$ N·m yr⁻¹. All of the synthetic catalogues have a nearly identical seismic moment rate, $2.3\times {10}^{20}$ N·m yr⁻¹ if all magnitudes are included in the calculation. When only the modelled events in the magnitude range of the instrumental catalogue (M5–8.17) are included, the seismic moment rates range from $4.34\times {10}^{19}$ to $9.81\times {10}^{19}$ N·m yr⁻¹, which are more comparable to the moment rate of the instrumental catalogue.

5.2 Scaling of coseismic slip, rupture area and stress drop

Average coseismic slip and rupture area are important earthquake source parameters for ground motion and tsunami simulation. Here, we check our synthetic model outputs for these values by comparing their scaling with recent empirical scaling laws for subduction interface earthquakes from Murotani et al. (2013), Allen & Hayes (2017), Thingbaijam et al. (2017) and Skarlatoudis et al. (2016). Stirling et al. (2024) combined the scaling relations from the studies we mentioned above and proposed the scaling relations of maximum and minimum values for application in 2022 NZ NSHM. We also include the scaling relations from Stirling et al. (2024) in our comparisons. The events in the synthetic catalogues have larger coseismic slip than the mean empirical scaling laws would predict (Fig. 5). Lower values of the state coefficient reduce this discrepancy and bring some of the models within the uncertainties of the empirical scaling laws. Compared to the homogeneous stress models, the slopes of the scaling of synthetic catalogues of the heterogeneous and brittle–ductile stress models are more like the empirical scaling laws. For the rupture area, the events in the synthetic catalogues from models with higher state coefficients (b = 0.004 and 0.003) have smaller rupture area than the scaling equations, while the models with state coefficient b of 0.002 and 0.0015 have rupture areas that match the empirical scaling equations (Fig. 6).

As we found for the scaling of coseismic slip, the stress drops also decrease with decreasing state coefficients b (Fig. 7). The modelled stress drops lie mostly in the range of global subduction interface events (∼0.4 to ∼40 MPa) derived by Ye et al. (2016a). We also find that the modelled average stress drop of our events are relatively stable over all magnitudes, which is similar to the stress drop of global subduction interface events in Ye et al. (2016a). At lower magnitudes, there is a larger variation in stress drops. The variation reduces with increasing magnitude; this could be related to small-scale variations in the heterogeneous stress field, which may be similar to the natural world. At larger magnitudes, these stress variations are smoother, yielding a smaller range of stress drop values (Fig. 7). The lack of systematic depth dependence of stress drops in Ye et al. (2016b) differs from the stress drops of the synthetic catalogues in this study, which have a clear depth dependence when the heterogeneous stress distribution is incorporated (Fig. S8, Supporting Information).

The scaling of the hydrostatic models is shown in Fig. S9 (Supporting Information). Due to the high initial effective normal stress in the hydrostatic models, the events have higher coseismic slip and stress drop and smaller rupture area than the overpressure models. The differences between the modelled scaling and the empirical scaling of the hydrostatic models are also larger than the overpressure models. The high-slip values generated in the hydrostatic models, especially the model with state coefficient (b) of 0.004, could cause overestimation of ground motion and tsunami height and make the hazard assessment less reasonable. Therefore, the hydrostatic models are not considered as preferred models in this study although the larger M_Maxs from these models could be helpful to the hazard assessment of extreme seismic events, particularly for the field of tsunami hazard assessment and disaster risk reduction.

5.3 Slip distributions of synthetic earthquakes

We have also tested our models by examining slip distributions of individual ruptures in the synthetic catalogues. Fig. 8 shows the slip distributions of magnitude 8.71–8.75 earthquakes in the synthetic catalogues from the models using rate and state coefficients of a = 0.001 and b = 0.004. Interface ruptures for the homogeneous stress model are similar with an approximately rectangular shape and a single elongate slip maximum (Fig. 8, left), while rupture models for the heterogeneous stress model have more variability in slip distribution and, in some cases, have multiple high-slip areas (Fig. 8, middle). The large deeper slip in the heterogeneous models arises due to the increase in stress with depth, as these models do not account for depth-dependent rheology. The brittle–ductile heterogeneous models correct for this effect and also retain the variability in slip distribution between individual ruptures (Fig. 8, right).

6 DISCUSSION

6.1 The effect of heterogeneous initial stresses

Our results show that prescribing spatially uniform initial stresses on faults in RSQSim tends to result in synthetic earthquake catalogues which have repeating (characteristic) ruptures and MFDs that are difficult to reconcile with observations from the instrumental seismicity record. By contrast, using spatially heterogeneous initial stresses generates earthquake catalogues that have more varied ruptures (both shapes and slip distributions), and MFDs that are in better agreement with the instrumental catalogue.

The reason why the spatial variation of initial stresses plays such an important role in the synthetic earthquake catalogue simulation is the difference of the evolving speeds of effective normal and shear stresses. Although normal and shear stresses both evolve when earthquakes are generated in RSQSim, the earthquake ruptures contribute more to changes in the shear stresses than the normal stresses. Also, since the initial shear stresses are just a fraction of the initial effective normal stresses, the relative perturbations are larger for the shear stresses. Therefore, the spatial distribution of initial effective normal stress dominates the shape of the resulting synthetic catalogue. To evaluate this point, we propose two test models by adding heterogeneity to either normal or shear stresses at one time, that is, heterogeneous normal stress and heterogeneous shear stress models, and examine their stress evolution in time. The spatial variation of initial normal stress forces the homogeneous shear stress to evolve into a heterogeneous distribution (with smaller values) even though the shear stress is initially homogeneous (Fig. S10, Supporting Information). On the other hand, the homogeneous initial normal stress makes the shear stress develop into a spatially less variable distribution at the end of the catalogue although the heterogeneity was added into the initial shear stress (Fig. S11, Supporting Information). The heterogeneous normal stress model also produces a more Gutenberg–Richter-like MFD, which is also closer to the MFDs of the heterogeneous models (Fig. 4), than the heterogeneous shear stress model (Fig. S12, Supporting Information).

The test of different state coefficients (b) with fixed rate coefficient of a = 0.001 indicates that smaller state coefficients result in smaller coseismic slip and larger rupture area. Using small values of the state coefficient helps bring the average coseismic slip and rupture area of ruptures within the uncertainty bounds from empirical scaling laws (Figs 5 and 6). Some models can even match the mean values of the empirical scaling laws; however, the plurality of smaller earthquakes and paucity of larger earthquakes in these same models contrast with the instrumental catalogue (Fig. 4). It is therefore necessary to strike a balance between fit to slip and area scaling relations on one side and the fit to observed seismicity and expected MFDs on the other. We note that the uncertainties presented in Figs 5 and 6 only include the standard deviation of the y-axis, and the uncertainties of the empirical scaling laws could be larger if the standard deviations of the regression coefficients are also included. At the same time, catalogue incompleteness may bring down the rates of moderate earthquakes in the instrumental Hikurangi–Kermadec catalogue, meaning that a comparatively higher rate of moderate earthquakes in a synthetic catalogue is not necessarily problematic. Considering these uncertainties and the similarities of the MFDs between the instrumental and synthetic catalogues, we suggest that the model with a state coefficient of 0.003 yields a reasonable compromise between fit to the instrumental MFD and fit to the slip and rupture area scaling relations.

6.2 Controls on synthetic earthquake rupture in subduction zones

The earthquake rupture of the subduction zone plays an important role for the estimation of seismic and tsunami hazard. In previous sections, we found that the initial stress distribution could be one of the factors that control the rupture on the subduction interface. Fig. 9 shows the ruptures along latitude of the large earthquakes (M ≥ 8.5) from the homogeneous, heterogeneous and brittle–ductile stress models with constant b of 0.004. The earthquake ruptures of the other models with constant frictional properties are presented in Figs S13–S15 (Supporting Information). We found more earthquake ruptures at the north Kermadec segment. Since the stress evolution in RSQSim is dominated by the stressing rate calculated from the slip-deficit rate, earthquakes tend to occur more at the patches with higher slip rates. Although all models have concentrated ruptures at the north Kermadec segment, the distribution of earthquake ruptures of the heterogeneous model is more smoothed along latitude. We also found that there are no earthquakes that rupture the entire subduction interface in any of the models. Only some of the M9⁺ earthquakes from the homogeneous stress model ruptured the whole Kermadec segment. Two rupture barriers at the intersection of the Hikurangi and Kermadec segments (−37.5°) and near Raoul Island (−32.5°) can be clearly seen in Fig. 9, especially in the homogeneous stress model (Fig. 9a). We found these barriers strongly correspond to the abrupt increases of the slip-deficit rates at these two locations. Although the slip-deficit rates have been smoothed when the Hikurangi and Kermadec segments were combined in the 2022 NZ NSHM (Van Dissen et al. 2022), the differences of the slip-deficit rate still act to stop the earthquakes rupturing the entire Hikurangi–Kermadec subduction zone. Despite the absence of the records of historical or palaeo-earthquakes that ruptured the entire subduction zone, it will be necessary to include the largest earthquakes to estimate the seismic and tsunami hazard more completely.

6.3 M_Max of earthquakes on the Hikurangi–Kermadec subduction interface

The estimated M_Max of the synthetic catalogues from the suggested model [heterogeneous stresses with constant state coefficient (b) of 0.003] is 9.17, which is much lower than the M_Max (9.5) that NZ NSHM estimated in their seismicity rate model (dashed grey lines in Fig. 4). Due to the different methods and parameters used between this study and the NZ NSHM, it is hard to pin down the main factor that caused the differences of the M_Maxs. However, as we discussed in the previous section, the preferred models have no earthquake that ruptures the entire subduction interface because of the barriers at the certain points of the subduction interface. The M_Max from NZ NSHM is an average of the inversion-based models with the constraint of the maximum rupture area of the subduction interface and an empirical scaling relation (Gerstenberger et al. 2022b). Considering the total area of the subduction interface of the geometry model used in this study is 378 406 km², the potential M_Maxs estimated by the empirical scaling relation from Allen & Hayes (2017) and Thingbaijam et al. (2017) are 9.59 and 9.35, respectively. If we examine the earthquakes generated in the spin-up period—before the model has introduced spatial heterogeneities in the stress state—there are a few earthquakes where the entire interface ruptures from the homogeneous models, with a magnitude range from 9.12 to 9.63. The models proposed in this study are therefore able to produce earthquakes with these potential M_Maxs if there is full rupture on the subduction interface. However, such full interface ruptures do not appear in the later parts of our catalogues, after the spin-up period. Future studies should seek to determine whether they do occur, and if so, how often.

6.4 Effect of the variation of (b − a) and pore-fluid pressure ratio on the earthquake depth distribution

The variation of (b − a) and pore-fluid pressure ratio controls the depth distribution of the earthquakes. We compare the density of the depth distribution of the earthquakes with magnitude ranges of 5–6, 6–7, 7–8, 8⁺ from the instrumental catalogue (black bars) and the synthetic catalogues (blue steps) in Fig. 10. Note that the instrumental catalogue only has one M8⁺ event, therefore the depth distribution in this magnitude range is for demonstration of the depth distribution of the large events in the synthetic catalogues.

Fig. 10(a) shows the depth distribution of the model with the heterogeneous stress but constant (b − a) and pore-fluid pressure ratio. High seismicity density can be observed at the top (0–15 km) and bottom (45–50 km) of the subduction interface in the synthetic catalogue. This behaviour makes sense given that in RSQSim, stresses are imposed by backslip loading; this combination of loading and small imposed effective normal stresses at shallow depths promotes earthquake nucleation in the model. By adding the variation to the (b − a) and pore-fluid pressure ratio, the peak of seismicity at the top of subduction zone moves down to 5–15 km, at the same time, the bottom peak of the seismicity is not seen in the model with variable frictional properties (Fig. 10b). The effect of the separate variation of only (b − a) and only pore-fluid pressure ratio on the depth distribution is shown in Fig. S16 (Supporting Information). The variations of (b − a) and pore-fluid pressure ratio play similar roles on the control of seismicity depth. Although the depth distribution of the model with variable frictional properties does not perfectly fit the observations, considering the depth uncertainty in the instrumental catalogue ranges mostly from 0 to 20 km (Fig. S17, Supporting Information), our model has provided a decent result by moving down the seismicity in depth.

The MFDs of the models with variation of frictional properties are shown in Fig. S18 (Supporting Information). The seismicity rates in lower magnitudes (M5.5–7) of the models with variation of frictional properties (green, blue and red lines) are higher than the instrumental catalogue. This difference is due to the low (b − a) values near the transition depth and the low effective normal stress at the top and bottom of the subduction interface. At the higher magnitudes (M7⁺), the MFDs of the models with the variation of frictional properties have higher deviation from the instrumental catalogue than the model with uniform frictional properties. Fig. S19 (Supporting Information) shows the scaling of average slip, rupture area, stress drop and depth dependence of stress drop for the model with variation of (b − a) and pore-fluid pressure ratio. Note that RSQSim by default does not include patches in State 3 when computing rupture area and moment; however, we include patches in State 3 that have slip during the duration of the event. By doing this, we allow ruptures to propagate into velocity strengthening regions, which would not happen otherwise. By adding the variation of frictional properties, we now have more average slips and rupture areas from the synthetic catalogue lying in the uncertainty range of the empirical scaling laws (Figs S19a and b, Supporting Information) comparing to the results of the models with uniform frictional properties (Figs 5 and 6). For the stress drop, the gaps of the synthetic earthquake (middle column in Fig. 7) and the lower limit of the global subduction interface events are filled (Fig. S19c, Supporting Information). Since (b − a) controls the value of stress drop in the rate-and-state friction law, the depth dependence of stress drop is identical to the (b − a) variation along the depth in the middle of subduction interface (the velocity strengthening portion).

Overall, incorporating both variations of (b − a) and pore-fluid pressure ratio into the stress models improves the results in the perspective of earthquake depth distribution and scaling of coseismic slip, rupture area and stress drop. Although there is still high seismicity at shallow depth, depending on the intended use case, the nucleation depth behaviour may not matter. For example, tsunami hazard analyses rely on slip distributions rather than hypocentres and seismic hazard analyses typically use earthquake magnitude and a 3-D representation of the rupture, especially for large earthquakes. For such analyses, the improvements to the synthetic catalogues discussed in Section 6.1 are likely more important than realistic hypocentres. Since we only consider simple 1-D models for the variation of the frictional properties along depth in this study, there are limitations on fitting all the components perfectly, such as the deviation of the MFDs of the synthetic catalogues and the instrumental catalogue. Future studies can further consider the variation of frictional properties along strike for a better match of the locations of earthquakes.

6.5 Raoul island doublets

In 2021, a pair of relatively large earthquakes (M_w 7.4 and 8.1) occurred within 2 hr near Raoul Island. Similar pairs of earthquakes were also recorded in 1917 and 1967 making the area around Raoul Island a repetitively ruptured asperity (Lythgoe et al. 2023). We scan the events in this area in our synthetic catalogues to examine whether the earthquake simulator can generate doublets like the 2021 events. To perform a more thorough search around Raoul Island, any pair of events with interevent time within 3 hr and locations in this area (181°–185° and −32° to −27°) is considered as ‘doublet’. Figs S20–S23 (Supporting Information) show the selected pair of earthquakes from the homogeneous, heterogeneous brittle–ductile stress models and the model with the variation of frictional properties. The pair of events from the homogeneous model has an M7 event as the first earthquake and an M9 event as the second earthquake (Fig. S20, Supporting Information). Since doublets usually consist of two events with similar magnitudes, the pair from the homogeneous model is more likely to be considered as a pair of foreshock and main shock. The pairs from the models with the variation of stresses and frictional properties have magnitudes that are more similar to each other and the 2021 doublet. This indicates that the earthquake simulator can generate doublets if the variations in stress and frictional properties are applied. However, the asperities of the coseismic slip are not fully locked in the source area of the 2021 doublet at Raoul Island (see Fig. 1 and Figs S23–S23, Supporting Information) because our models do not contain the variation of stress and frictional properties along strike. A more thorough study of the characteristics of the hazards of doublets along the Hikurangi–Kermadec subduction interface could be performed in the future when the variation of the frictional properties along strike is applied to the earthquake simulator.

7 CONCLUSIONS

First-generation multicycle earthquake simulators typically use uniform input parameters, with initial stresses in these simulations typically prescribed as a single uniform value. Here, we test the role of spatial heterogeneity in the initial stresses for RSQSim earthquake cycle simulations on the Hikurangi–Kermadec subduction zone, Aotearoa New Zealand. Heterogeneous stress was calculated from the three principal stresses on each fault patch (10 km side length) and varying rate-and-state coefficients; we test four (b − a) values by fixing a to 0.001 and varying b values to 0.004, 0.003, 0.002 and 0.0015. The results for heterogeneous and homogeneous initial stress models are compared with the MFDs for the instrumental earthquake catalogue and empirical scaling laws from global earthquakes. Our results suggest stress heterogeneity produces more realistic and less characteristic synthetic catalogues than the homogeneous model. Heterogeneous stresses with pore-fluid overpressure (uniform pore-fluid pressure ratio) and constant state coefficients (b) of 0.003, produce the best fit of the heterogeneous stress models to MFDs and empirical scaling laws. We further add heterogeneity to the frictional properties, the rate-and-state coefficients (b − a) and pore-fluid pressure ratio. The variation of the frictional properties along depth improves the fitting to the depth distribution of earthquake nucleating points and the empirical scaling laws. Heterogeneity of other parameters (e.g. friction coefficient, reference slip rate, characteristic distance, initial state variable, etc) within RSQSim could further improve the applicability of the synthetic earthquake catalogues to seismic hazard problems and form the focus of future research.

We also examine the factors that affect the earthquake ruptures in our results. None of the models generate earthquake ruptures along the entire subduction interface due to the abrupt increases of slip-deficit rate at the intersection of the Hikurangi and Kermadec segments and near Raoul Island. This might be one of the reasons that the M_Max of the synthetic catalogues from the suggested models (9.17) are smaller than the potential M_Maxs estimated from the empirical rupture area scaling equations (9.59 and 9.35) and the M_Max estimated by NZ NSHM (9.5).

ACKNOWLEDGEMENTS

We would like to express our sincere gratitude to Resilience to Nature's Challenges for the support of this research project. Their funding made it possible to conduct this study and achieve our objectives. We also extend our sincere appreciation to the editor and the anonymous reviewers for their insightful comments and constructive feedback on our manuscript. Their thoughtful critiques significantly enhanced the quality of our work and helped clarify our findings. We are grateful for their time and effort in the review process. We would also like to thank our colleagues in GNS Science and University of Canterbury for their invaluable support and collaboration.

DATA AVAILABILITY

The data underlying this article are available in GNS Science data set catalogue, at https://doi.org/10.21420/k8ab-2671.

REFERENCES