Impact of Variable Fault Geometries and Slip Rates on Earthquake Catalogs From Physics‐Based Simulations of a Normal Fault

Physics‐based earthquake simulators have been developed to overcome the relatively short duration and incompleteness of historical earthquake and paleoseismic records, respectively. These simulators have the potential to be a useful addition to seismic hazard assessment as they produce millions of synthetic earthquakes over thousands to millions of years using predefined fault geometries and slip rates. Due to the sparsity of fault data and the computational expense of the modeling, it is common to simplify earthquake simulator input parameters. Fault surfaces are typically characterized by simplified mapped traces and constant dip with depth, with slip rates that are often assumed to be uniform over the fault plane. This study uses the well‐defined 3D geometry of an active normal fault in offshore Aotearoa New Zealand, the Cape Egmont Fault, to demonstrate the impact of using nonuniform slip rates and nonplanar fault geometries on the resulting synthetic earthquakes from the earthquake simulator RSQSim. Adopting variable slip rates that decrease to zero along the fault tipline (rather than uniform slip rates) reduces unrealistically high nucleation rates of seismicity along fault edges. Introduction of complex 3D fault geometries, including fault segmentation and bends on kilometer scales, and variable slip rates produces less characteristic earthquake populations, increasing the number of M6‐7 moderate‐large magnitude events. Incorporation of variable fault geometries and slip rates in physics‐based simulators may modify the seismic hazards estimated from synthetic earthquake catalogs.

• Variable, rather than constant, slip rates reduce high seismicity rates along the fault edges and can improve the depth distributions of events • Introduction of complex 3D fault geometries produces less characteristic earthquake populations • Variable fault geometries and slip rates in physics-based simulators can improve the synthetic earthquake catalogs Supporting Information: Supporting Information may be found in the online version of this article.
simulators allow us to consider a wide range of possible slip distributions for future large earthquakes, including multi-fault earthquakes, and therefore provide additional information to constrain associated uncertainties (Field, 2019).Finally, earthquake simulators represent an opportunity to test or augment more conventional seismic hazard models, which use the same inputs (fault geometry and slip rate), but rely on empirical relationships rather than approximations of the physics that govern fault slip (Shaw et al., 2018).However, to maximize the utility of simulators for these purposes, the simulators and resulting synthetic earthquakes need to be as realistic as possible.
Due to the difficulties associated with detailed characterization of the three-dimensional fault properties, such as, geometry, slip rates, and frictional behavior, in combination with the computational expense of earthquake simulators, it is common to simplify input parameters by, for example, using planar faults with uniform slip rates (e.g., Tullis et al., 2012).Previous studies of faults, however, indicate that their geometries are rarely planar.Instead, faults have been shown to have complex, often segmented, surfaces on a wide range of scales in both map-view and cross-section (e.g., Childs et al., 1995;Manighetti et al., 2015Manighetti et al., , 2021;;Peacock & Sanderson, 1991;Roche et al., 2021).Furthermore, detailed three-dimensional mapping indicates that many faults comprise complex 3D anastomosing structure of sub-parallel and bifurcating fault surfaces, that locally vary in both strike and dip, and display varying degrees of connectivity (e.g., Camanni et al., 2019;Delogkos et al., 2020;Deng et al., 2020;Marchal et al., 2003;Roche et al., 2021;Walsh et al., 2003).These changes in fault geometry have the potential to impact the size and magnitude-frequency distribution (MFD) of both natural and synthetic earthquakes (Howarth et al., 2021;Schwartz & Coppersmith, 1984), and their application to simulators warrants further investigation.
Spatially uniform or near uniform slip rates imposed on fault surfaces for the purposes of simulating earthquakes are also a simplification: one which has the consequence of generating infinite slip gradients and associated modeled stress singularities near fault tips.Slip rates on the modeled faults that serve as inputs for earthquake simulators often taper to zero at fault edges to avoid the creation of such stress singularities (e.g., Howarth et al., 2021;Shaw, 2019), but this tapering is rarely well constrained by observations.In nature, observed incremental slip during earthquakes, cumulative fault slip on geological timescales, and fault slip rates all vary along faults, typically reaching a maximum toward the fault center with a gradual decrease toward the fault tips with distributions spanning from strongly asymmetric to symmetric (e.g., Manighetti et al., 2005Manighetti et al., , 2001;;Marchal et al., 2003;McPhillips & Scharer, 2023;Nicol et al., 2005Nicol et al., , 2020;;Roberts & Michetti, 2004;Torabi et al., 2019;Walsh & Watterson, 1988;Wesnousky, 2008).Departure from this first-order pattern of slip and slip rates can occur along fault traces at geometrical irregularities such as fault bends or segment boundaries (e.g., Iezzi et al., 2018Iezzi et al., , 2020;;Walker et al., 2009).The impact of these irregularities may be transient, while slip rates on individual faults can also vary temporally on a range of timescales (Mouslopoulou et al., 2012(Mouslopoulou et al., , 2009;;Nicol et al., 2009).These slip-rate variations also have the potential to impact synthetic earthquake catalogs and their inclusion in earthquake simulators also requires investigation.
In this study, we examine the effect of fault geometry and spatial variability of fault slip rates on the synthetic earthquake catalogs produced by RSQSim, a physics-based earthquake simulator (Richards-Dinger & Dieterich, 2012).Published studies demonstrate that the MFDs of simulated earthquakes using RSQSim can be similar to historical earthquakes in many regions of active faulting (e.g., Herrero-Barbero et al., 2021;Shaw et al., 2022).These similarities support the use of RSQSim for either estimating seismic hazard or for constraining parameters input to seismic hazard models.Shaw et al. (2018) found that RSQSim simulations for the Californian fault system can replicate the seismic hazard statistics from the long-term Third Uniform California Earthquake Rupture Forecast (UCERF3; Field et al., 2014).Similarly, Milner et al. (2021) employed RSQSim to establish a nonergodic framework for probabilistic seismic-hazard analysis, while Chartier et al. (2021) utilized RSQSim as an auxiliary tool for weighting a logic tree that explores seismicity uncertainties in the Sea of Marmara.These papers demonstrate the utility of RSQSim and coupled with its increasing availability suggest that it will continue to have seismic hazard applications.Here we use the 3D geometry and slip rates of an active normal fault in offshore Aotearoa New Zealand, the Cape Egmont Fault (CEF), which is well characterized from seismic reflection data and stratigraphic wells (Nicol et al., 2005;Seebeck et al., 2020Seebeck et al., , 2021)).Simulated earthquakes generated using geophysically defined realistic fault geometries and slip rates are compared to output from models with more planar geometries and uniform slip rates.The results of this study demonstrate that adding variability to the input fault geometries and slip rates can modify the resulting RSQSim earthquake catalogs, changes which have implications for the understanding of seismic processes and hazards determined from earthquake simulators.

The Cape Egmont Fault
The active CEF is located in the offshore Southern Taranaki Basin southwest of the Taranaki Peninsula, North Island, Aotearoa New Zealand, within the continental crust of the Australian Plate (Figure 1a).The CEF is a ca.70 km long, northeast to east-northeast trending, normal fault.Contemporary extension in the Taranaki Basin is demonstrated by historical seismicity (up to M w 5.4) and onshore active faults which indicate earthquakes that ruptured the ground surface (Robinson et al., 1976;Sherburn & White, 2005;Townsend et al., 2010).Offshore, a seabed scarp (up to ca. 6 m in height) indicates seafloor rupture along the CEF in one or more moderate to large magnitude earthquakes up to M w ∼ 7.3 during the Late Holocene (Nicol et al., 2005;Nodder, 1993Nodder, , 1994;;Thingbaijam et al., 2017).
The CEF is well imaged by geophysical data and has been extensively studied (King & Thrasher, 1996;Nicol et al., 2005;Nodder, 1993Nodder, , 1994;;Seebeck et al., 2020Seebeck et al., , 2021)).The 3D structure of the CEF is resolved to a depth of ∼8 km using high-quality industry 2D and 3D seismic-reflection lines covering a region 80 km long and up to 35 km-wide (Figure 1; Seebeck et al., 2020Seebeck et al., , 2021)).The CEF is a Late Cretaceous-Palaeocene normal fault that  et al., 2021).GeoNet earthquakes from 1970 to 2020 occurring at depths less than 33 km and magnitudes greater than M L 3, including the 1853 M w 6.8 (gray-filled circle) pre-seismograph network event, the 1974 M w 5.4 relocated event (red-filled circle; Webb & Anderson, 1998), and the recent M w 3.7-4.7 Regional Moment Tensor focal solution events (blue-filled circles; Ristau, 2008).Inset map is a depth structure map of Base Pliocene unconformity (after Seebeck et al., 2021).(b and c) Interpreted seismic sections of the CEF and some horizons (after Seebeck et al., 2020).See (a) for location of sections.
was inverted (i.e., reverse reactivated) during the Miocene and has subsequently accommodated a further phase of Plio-Pleistocene to Recent normal faulting (Nicol et al., 2005;Reilly et al., 2015;Seebeck et al., 2020), which is the subject of this paper.The geometry of the CEF fault surface changes along strike and down dip, therefore providing the opportunity to test how introducing a realistic fault geometry influences simulated earthquake catalogs (compared with, e.g., a simple planar approximation to fault geometry).Northern segments have steep easterly planar dips (ca.60-70°), while southern segments appear to have more listric geometries (Seebeck et al., 2021).Prominent changes in fault strike occur at a high-angle junction between the northern (characterized by multiple synthetic and antithetic faults) and the southern segments along fault branch lines (Nicol et al., 2005;Seebeck et al., 2021).These changes in fault geometry at least partly reflect differences in the reactivation and inversion history of the fault since the Late Cretaceous (Nicol et al., 2005;Seebeck et al., 2020).The most recent phase of normal faulting on the CEF commenced at ca. 3.7 Ma and resulted in a maximum throw of ca.2.3 km on base Pliocene strata approximately in the center of the fault (Nicol et al., 2005;Seebeck et al., 2020Seebeck et al., , 2021; see inset of Figure 1a and seismic sections in Figures 1b and 1c).
Historical seismicity indicates that most earthquakes in the region of the CEF occur beneath the Late Cretaceous to Holocene sedimentary sequence, at depths >∼5 km, and the maximum depth of fault rupture is ca.20 km, above which 90% of relocated earthquake hypocenters occur (Ellis et al., 2021;Seebeck et al., 2021; Figure 2).This lower limit of seismicity in the upper crust is attributed to the transition from brittle to ductile deformation where crustal materials reach their maximum yield strength (Hauksson & Meier, 2019;Sibson, 1982; see Figure 2).Consequently, the strength of the fault increases linearly with depth in the brittle (frictional) zone until the transition to ductile deformation (e.g., aseismic creep), below which yield strength diminishes exponentially with increasing depth (Figure 2).Sibson (1983) and Scholz (1998Scholz ( , 2002)).(c) Yield strength profiles for quartz-and diorite-dominated crustal compositions for heat flows of 59 and 65 mW/m 2 (solid and dashed lines, respectively; after Hauksson & Meier, 2019) consistent with the Taranaki region (Funnell et al., 1996).Temperature-depth relationship estimated using the 60 mW/m 2 temperature model of Funnell et al. (1996). 10.1029/2023JB026746 5 of 17 The dip of the fault between ∼8 and 20 km depth is constrained by structural modeling using the geometries of Cretaceous and younger strata displaced by the fault; below 20 km depth, aseismic creep is assumed for the fault (Figure 2; Seebeck et al., 2021).Only the large-scale fault surfaces were utilized for this study with many shorter length faults (e.g., <5 km), either mapped using the available seismic reflection surveys or sub-seismic scale, omitted.In total, four main, two synthetic and two antithetic fault surfaces were included in the model (Figure 3a-the "realistic fault geometry").The corresponding slip rates along each fault surface were calculated from displaced strata dated using micropaleontological ages from exploration wells (Table 1; Seebeck et al., 2021).More specifically, slip rates can be calculated in circumstances when the long-term sedimentation rates in the hanging wall of a fault are greater than fault slip rates.The CEF fulfills this condition during the post-Miocene period (Nicol et al., 2005), with successively older syn-faulting horizons having progressively greater slip with the difference in slip between two horizons principally reflecting the horizon ages and the faultslip rate (e.g., Childs et al., 2003;Petersen et al., 1992).

RSQSim Earthquake Simulator
To simulate earthquakes, we use RSQSim, a multicycle physics-based earthquake simulator developed by Dieterich and Richards-Dinger (2010).RSQSim is a computationally efficient boundary element, event-driven, three state algorithm based on approximations of the rate-and-state friction equations, which allows the generation of synthetic earthquake catalogs over periods of hundreds of thousands of years (Richards-Dinger & Dieterich, 2012).RSQSim utilizes the following physics: (a) elements on individual fault surfaces or between faults interacting with quasistatic elastic stresses (dynamic stresses are neglected), (b) rate-and-state frictional behavior is simplified into a three-regime system in which elements are either stuck, nucleating, or sliding dynamically, and (c) dynamic sliding in which slip rate is fixed at a constant value.Faults are loaded using the backslip method, where slip rates are specified on each fault surface and based on their static stress interactions, stressing rates applied to approximately reproduce the specified slip rates (Savage, 1983).The backslip method produces stress singularities at fault tip lines for uniform slip-rate faults, leading to unrealistic seismicity nucleating at these locations (Shaw, 2019).Here we examine whether these stress singularities can be addressed by applying variable slip rates which decrease gradually to zero at fault tips.
We have used eight RSQSim models to examine how changes in fault geometry and slip-rates modify earthquake catalogs.These models comprise four different fault geometries, which are: (a) "realistic fault geometry" comprising the large-scale fault surfaces (e.g., >5 km) consisting of four main fault surfaces, including two synthetic and two antithetic to the main fault surfaces, with a mix of planar and listric dips (Figure 3a), (b) "simplified main fault surface" where the main fault trace was extrapolated to 20 km depth with constant dip of 60° (Figure 3b), (c) a "planar main fault surface" with a uniform strike and constant dip of 60° (Figure 3c), and (d) a "listric main fault surface" with uniform strike and 70° dip at shallow depths that gradually decreases to 25° at 20 km depth (Figure 3d).Each fault surface was constructed using Coreform Cubit software and comprises triangular elements with an area of ca. 1 km 2 (e.g., Figure 4a).The use of triangular elements with a size of ca. 1 km 2 or less improve the synthetic earthquake catalogs with reduced stress concentrations and more realistic ruptures (Gilchrist, 2015;Richards-Dinger & Dieterich, 2012).We chose to use 1 km 2 -area triangles for two reasons.First, this size makes our results consistent with regional or national scale applications of earthquake simulators.The fault models used in such simulations include tens or hundreds of faults-often simplified representations-and adopt a triangle area of 1 km 2 or greater to make the model computationally tractable (e.g., Shaw et al., 2022 use a triangle area of ∼1.7 km 2 ).We do not expect triangle size to influence our results significantly, but in the absence of a formal sensitivity test of triangle size, we choose our mesh size to be consistent with this intended use case.
Second, using a triangle area of 1 km 2 made it feasible to run a range of models to explore sensitivity of results to model parameters (2-5 CPU hours per model run), initially following a trial-and-error approach but also for the more formal sensitivity test (60 model runs) in the supplementary data.
Each of the four geometric fault-surface models was also run with (a) uniform and (b) variable slip rates (e.g., Figure 4 and Figures S1-S3 in Supporting Information S1) to investigate the effect of different slip-rate distributions on the resulting synthetic earthquake catalogs.In the upper ∼6 km of crust the nonuniform slip-rate models are constrained by measurements from seismic reflection lines.These lines indicate that the maximum slip-rates are approximately at the center of each fault surface in their strike dimension (e.g., symmetrical slip rate distribution).Deeper in the crust we use the seismicity rate as a proxy for slip rate and locate the maximum slip rate at ∼16 km depth, where the seismicity peaks (e.g., Eberhart-Phillips & Reyners, 2022; Sherburn & White, 2005; Figure 2).For the purposes of this study, slip rates were assumed to decrease linearly from the maximum to the fault tips, as is often observed for geological displacements (e.g., Figure 4b and Figures S1b-S3b in Supporting Information S1; Dawers et al., 1993;Manighetti et al., 2015;Nicol et al., 2005;Walsh & Watterson, 1987).The average slip rate for variable slip-rate models is approximately equal to the value for uniform slip-rate models (e.g., Figure 4a and Figures S1a-S3a in Supporting Information S1), with the earthquake simulator producing  Fault frictional properties are assumed to be uniform for all realizations, with the coefficient of friction of 0.6.Despite the fact that the rate-and-state friction parameters and the initial stress conditions are also assumed to be uniform, a sensitivity analysis was performed to identify the best combination of these parameters that produces synthetic catalogs best fitting the characteristics of natural seismicity (e.g., magnitude-rupture area relation).
The sensitivity analysis was carried out using the "realistic fault geometry" with variable slip rates (Figure S1b in Supporting Information S1) and the parameters listed in Table 2.The outputs of a total 60 RSQSim simulations are shown in Figure 5.The synthetic catalog that best correlates with the magnitude-rupture area scaling relation (Figure 5b; Gerstenberger et al., 2022;Stirling et al., 2021) and with the MFD of the natural seismicity (Figure 6) is highlighted in red in Figure 5.The combination of parameters that resulted in this catalog have rate-and-state friction parameters a = 0.001 and b = 0.004, an initial shear stress of 60 MPa and a normal stress of 120 MPa.These normal and initial shear stress values are similar to those used in previous studies employing RSQSim (e.g., Herrero-Barbero et al., 2021;Howarth et al., 2021;Shaw et al., 2022) and we use them for consistency with that previous work.We note that the normal stresses are higher than those employed by some other quasi-static and dynamic multi-cycle earthquake simulators (50 MPa; Cattania, 2019;Lapusta & Liu, 2009;Perry et al., 2020;Thomas et al., 2014).This difference may be related to the way the rate-and-state equations are simplified differently in different models; it is not important for our results because in the RSQSim simplification there is a trade-off between normal stress and the imposed rate-and-state (a-b).Finally, a dynamic stress overshoot factor of 10% and an alpha-coefficient α = 0.25 were applied in order to deal with the dynamic limitations (Gilchrist, 2015;Herrero-Barbero et al., 2021;Richards-Dinger & Dieterich, 2012).
Utilizing the aforementioned parameters, we have run eight RSQSim simulations to examine how changes in fault geometry and slip-rates modify earthquake catalogs.Each of the eight simulator models produce earthquake catalogs spanning 500 kyr producing an average of ca.43,000 earthquakes  The RSQSim input parameters, fault surface geometries and properties, together with the resulting synthetic earthquake catalogs are available in a public data repository (Delogkos et al., 2023).

Results
The eight RSQSim models display different spatial distributions of earthquake nucleation locations, MFDs, and maximum magnitudes.
The catalog of simulated earthquakes changes with increasing fault geometry and slip-rate variability.The planar fault model produces a characteristic earthquake catalog with a linear increase in the number of earthquakes from ∼M w 4.5 to M w 6 and a more modest decrease in the numbers of earthquakes from ∼M w 6 to M w 7.3 (Figures 9c and 9d).By contrast, the "realistic fault model" with variable slip rates produces an approximately Gutenberg-Richter distribution (Figures 9c and 9d).The tendency toward more Gutenberg-Richter distributions with increasing complexity in fault geometry and slip rates suggests that the size distributions for simulated earthquakes on the CEF is strongly dependent on the model input parameters, with Gutenberg-Richter scaling a fundamental characteristic of natural earthquake systems.
To quantify the differences between the MFDs of our eight model runs, we calculated overall Gutenberg-Richter b-values for each model.We also calculated and compared the slope of the MFD for each model over two different magnitude intervals: between M w 5.2-5.7 and M w 6.2-6.9 (Table 3).These magnitude intervals were chosen arbitrarily following inspection of the model MFDs, to help measure how characteristic (or not) our model MFDs are.Overall b-values were calculated between a minimum M w of 5.2 (below which steps in the MFDs are introduced by the discretization of triangles in our mesh), and a maximum M w of 6.9.The overall b-values are somewhat lower than those estimated for the Cape Egmont Fault Zone by Sherburn and White (2005) (∼0.9-1.0 for our preferred models compared with 1.3 for Sherburn and White).However, we note that their b-value is calculated only for smaller earthquakes in an onshore region north-east of our study area, and we suggest that our b-values for the offshore fault system are comparable to other fault systems in New Zealand and worldwide (e.g., Stirling et al., 2012).For all models, the slope of the MFD is steeper from M w 5.2-5.7 than M w 6.2-6.9, suggesting a variable degree of characteristic behavior; for a truly Gutenberg-Richter MFD, the slope across these two intervals would be identical.The ratio of MFD slopes (slope 5.2-5.7 /slope 6.2-6.9 ; slope ratio hereafter) for each model also gives an indication of how characteristic that model is (Table 3).A perfectly Gutenberg-Richter MFD would  have a slope ratio of 1 whereas, more characteristic MFDs produce a slope ratio >1.For all of our models, the slope ratio is >1, indicating some level of characteristic behavior.Nevertheless, for all geometries, introducing a variable slip rate brings the slope ratio closer to 1, indicating an MFD that is closer to Gutenberg-Richter and less characteristic.Similarly, the realistic geometry (for both constant and variable slip-rate distribution) yields a slope ratio that is closer to 1 than for the other geometries.Our model with variable slip rates and realistic fault geometry has the least characteristic MFD of all our models, with a slope ratio of 1.38.
We note that even with a realistic fault geometry and the variable slip-rate distribution, the MFD is still somewhat characteristic.This (albeit lessened) characteristic behavior may be a consequence of the surface smoothness of the segments that make up the realistic fault model; introducing finer-scale geometric roughness (possibly random or fractal) to fault surfaces might make the MFD more linear and less characteristic.
Changes in the size distribution of earthquakes between models are accompanied by variations in the maximum seismic moments.Comparison of the models indicates that the largest seismic moments correspond to the "listric main fault surface" for both constant and variable slip rates (Figure 10).This increase in moment arises because all fault geometries (Figure 6) extend down to 20 km depth and, therefore, the "listric main fault surface" has the largest available area to rupture (Table 1, Figure 10).
All fault geometries with uniform slip rates except the "realistic fault geometry" rupture the whole fault surface with a wide range of seismic moments (Figure 10a).The maximum values of both seismic moment and ruptured area, however, tend to be relatively smaller for the faults with variable slip rates (Figure 10b) than with uniform slip rates (Figure 10a).In all four fault geometries, the decrease in the maximum values of seismic moment is also greater than the decrease in ruptured area.The decrease in these maximum values is a consequence of the spatial variability in imposed slip rate, and specifically for this study, the gradual decrease in slip rate toward the fault surface edges.These slip decreases (tapering) restrict the area that ruptures during the largest earthquakes with these events extending to the fault edges with tapered slip.Consequently, the recurrence interval of the largest events is smaller for faults with variable slip rates than the faults with constant slip rates (see the steps at the curve of the cumulative seismic moment in Figure 7).In other words, for the same time interval (e.g., 500 kyr), a fault with variable slip rates produces more frequent and smaller major events than a fault with constant slip rates (Figure 7).Interestingly, the "realistic fault geometry" presents several multi-fault ruptures (e.g., Figure 11) that result in a relatively wide range of maximum values of both seismic moment and ruptured area (Figure 10) depending on the different combinations of ruptured fault surfaces.In contrast to the other three fault geometries, the "realistic fault geometry" with variable slip rates also presents significantly lower maximum values of ruptured area than with uniform slip rates.This is also reflected in the smaller number of fault surfaces that participate in a multi-fault rupture.In the case of, for example, variable slip rates, only two events involved synchronous rupture of more  than two fault surfaces with a minimum moment magnitude threshold of M w 6.0 on more than one segment, whereas in the case of uniform slip rates, the corresponding number of events was 24 (Figure S10 in Supporting Information S1).This is likely due to tapering of coseismic slip toward the edges of the fault surface(s) resulting in the rupture being restricted and propagating to adjacent fault surface(s).
Uniform slip rates produce stressing rate singularities at the lower and lateral fault surface boundaries resulting in higher frequencies and magnitudes of earthquakes toward these fault boundaries (e.g., Figure 7).Compared to historical earthquakes in CEF region, uniform-slip synthetic catalogs also have unrealistic hypocenter depth distributions with most events nucleating at the base of the fault (ca.20 km depth) (Figure 9a).In contrast, the variable slip rates shown in Figure 4b (and Figures S1b-S3b in Supporting Information S1) appear to correct the issues associated with singular stressing rates at the edges of the fault surfaces (Figure 7).Variable slip rates also result in a more realistic depth distribution of seismicity (e.g., Eberhart-Phillips & Reyners, 2022;Sherburn & White, 2005; Figure 2), with the majority of events in the variable slip rate model originating at depths of 15-17 km (Figure 9b).A relatively strong peak, however, is observed at these depths, reflecting the corresponding peak in slip rates (e.g., Figure 4b and Figures S1b-S3b in Supporting Information S1); a smoother transition in slip rates with depth would likely reduce this effect.Slip-rate gradients can be specified for the variable slip-rate models which would produce results comparable to the hybrid loading method, which smooths the stress singularity along the fault edge (Shaw, 2019).

Discussion
A key aim of this paper was to test whether introducing more geometric and slip-rate complexity into fault-source models would significantly impact the output of an earthquake simulator.Our analysis demonstrates that the MFDs of synthetic earthquake catalogs are critically dependent on the fault geometries and slip rate distributions input to fault-source models.The results of this paper show that using nonuniform ("realistic") fault geometries and slip rates typically produces more Gutenberg-Richter synthetic earthquake catalogs than for planar faults.For the "realistic fault geometry" of the CEF, only the large-scale fault surfaces were included in the model with many shorter length faults (e.g., <5 km), either mapped using the available seismic reflection surveys or sub-seismic scale, omitted.The effect on earthquake behavior of smaller, sub-seismic scale faulting was not examined.However, a higher frequency of events at shallower depths of the "realistic fault geometry" simulation (Figure 9b) suggests that the incorporation of the two relatively small-size antithetic fault surfaces at shallower depths toward the north-eastern part of the CEF (Opunake 1 and 2 segments, Figure 3a) can have a significant impact on the resulting earthquake catalog.This is also confirmed by running an additional RSQSim simulation of the "realistic fault geometry" without including these two antithetic faults (Figure S11 in Supporting Information S1) indicating lower frequency of events at these depths but without significantly modifying the MFD.Furthermore, we have considered very simplified slip rate distributions along the active fault surfaces without taking into account (a) elevated slip rate gradients resulting in asymmetric slip rate distributions (e.g., Manighetti et al., 2005Manighetti et al., , 2001) ) due to fault interactions (e.g., Childs et al., 1995) and (b) slip rate variations due to other geometrical complexities such as fault bends (e.g., Walker et al., 2009).Larger variations in fault geometry and slip rate distributions (such as asymmetric slip rate distributions and elliptical or bell-shaped slip rate profiles) could, however, modify the MFDs and could be investigated in future RSQSim models.3a), highlighting (a) a M w 6.9 multi-fault rupture, and (b) a M w 6.7 single segment rupture, simulated using RSQSim.Fault patches that did not participate in the rupture are shown in gray.

10.1029/2023JB026746
13 of 17 Many other fault properties could also influence the resulting simulated earthquake statistics, including slip rake (slip direction), normal and shear stress conditions, fault frictional properties and fault roughness, all of which are here assumed to be homogeneous along the fault surfaces.A recent study of Bayesian finite-fault inversion has shown that variable fault slip directions allow better determination of the nonplanar fault geometries associated with large earthquakes (Wei et al., 2023).In addition, frictional properties of large seismogenic faults may vary spatially and/or temporally (e.g., before, during, and after an earthquake).Earthquake rupture nucleation and propagation can be controlled by the spatial variation of frictional properties and normal stress conditions (e.g., Niemeijer & Vissers, 2014;Ray & Viesca, 2017).Furthermore, contrary to a recent study by Beeler (2023), previous studies suggest that fault-surface roughness can also affect earthquake rupture processes (e.g., Zielke et al., 2017).The gradient of structural maturity has also been shown to impact the earthquake behavior (Perrin et al., 2016).A recent study by Allam et al. (2019), based on long-term simulations using RSQSim at a constant backslip rate, demonstrated quantitative links between fault roughness and earthquake properties including hypocenter location, coseismic slip variability, and characteristic slip distributions.These conclusions are backed up by our sensitivity analysis (Figure 5) of a range of RSQsim parameters (uniform distribution of rate-and-state friction (a and b), initial normal (σ 0 ) and shear (τ 0 ) stress), which may be approximately equally important for controlling the MFDs of synthetic earthquake catalogs.The analysis shows that different value combinations produce catalogs with significant variations in MFDs (Figure 5a) and magnitude-rupture area scaling relations (Figure 5b).
In many natural fault systems, it can be challenging to specify how these fault properties vary in three-dimensions.
In the absence of such information, imposing (a) simplified slip-rate distributions that decrease toward the fault tip lines and (b) observed empirical variations and scaling relations of the fault-surface segmentation/roughness properties (e.g., Manighetti et al., 2021) could produce more realistic synthetic earthquake catalogs than planar faults with uniform slip rates.It is worth noting, however, that in our case study, the spatial variability of the slip rate seems to have a stronger effect on the synthetic earthquake catalogs than the fault geometry, with the only exception of the listric fault geometry that presents a significant MFD discrepancy from the other three fault geometries with variable slip rate distributions (Figure 9).On a regional or national scale, our results have important implications for the use of earthquake simulators like RSQSim in seismic hazard models.A major feature of the MFDs of California-and NZ-wide synthetic earthquake catalogs derived using RSQSim (Shaw, 2019;Shaw et al., 2022) is characteristic behavior, with a flattening of the MFD above M w ∼ 7.0.This characteristic MFD influences modeled seismic hazard from these synthetic earthquake catalogs, and differs significantly from the linear MFDs of the seismic hazard models for California and New Zealand (Field et al., 2014;Gerstenberger et al., 2022).It is possible that the characteristic MFD of the California RSQSim catalog (Shaw, 2019) is partly caused by using simplified fault surfaces.Adding realistic geometric roughness on a range of scales (either constrained by data or randomly) might make the MFD of individual faults less characteristic.In turn, this would likely lead to a less characteristic overall MFD for the whole model, which could potentially improve modeled seismic hazard results.
Collectively, published studies and this paper suggest that introducing more variability to the input parameters tends to produce more small-moderate earthquakes, smaller maximum seismic moments, and greater slip variability.One of the outcomes of more heterogeneous input parameters are synthetic earthquake catalogs suggesting more Gutenberg-Richter rather than characteristic scaling.Future studies are required to determine what level of variability of the input parameters is required to generate results that have the potential to advance understanding of earthquake processes and seismic hazards.

Conclusions
This study demonstrates that the use of realistic 3D fault geometries with nonuniform slip rates has a significant impact on the simulated earthquake catalogs produced by a physics-based earthquake simulator.
• Variable slip rates decreasing toward fault tips appear to resolve issues associated with singular stressing rates at fault surface boundaries, and also produce depth distributions of simulated seismicity better matching natural seismicity.• The realistic 3D fault geometry also produces better depth distributions of seismicity but also, in combination with variable slip rates, provides better-defined Gutenberg-Richter distributions of events than the other simplified representations of the fault.
• The maximum seismic moment decreases for the faults with variable slip rates (compared to uniform slip models) due to a decrease in maximum rupture area.
However, because the determination of detailed 3D fault geometries and slip rate spatial variations is often associated with a high degree of uncertainty, the earthquake simulator results should be used with due care in seismic hazard assessment, ensuring that key epistemic uncertainties that affect hazard outcomes have been reasonably well explored.This project has received funding from the Irish Research Council under Grant GOIPD/2020/530 and from Science Foundation Ireland co-funded under the European Regional Development Fund (Grant 13/RC/2092).It was also partly supported by a grant from the New Zealand Resilience to Natures Challenges programme "Earthquakes and Tsunami" theme (contract C05X1901).Mapping of the Cape Egmont Fault was funded by the New Zealand Ministry of Business, Innovation and Employment Endeavour Fund -Smart Ideas (contract C05X1801) under the project title "Geological modelling of active fault systems."Academic Licenses of the software Petroleum Experts MOVE provided to the University College Dublin are kindly acknowledged, as are the RCC facilities at the University of Canterbury.We also thank two anonymous reviewers for their constructive comments, Alice-Agnes Gabriel and Isabelle Manighetti for their valuable feedback, comments, and editorial handling.Open access publishing facilitated by University of Canterbury, as part of the Wiley -University of Canterbury agreement via the Council of Australian University Librarians.

Figure 1 .
Figure 1.(a) Cape Egmont Fault (CEF) with identified and inferred active faults marked by solid and dashed red lines, respectively (afterSeebeck et al., 2021).GeoNet earthquakes from 1970 to 2020 occurring at depths less than 33 km and magnitudes greater than M L 3, including the 1853 M w 6.8 (gray-filled circle) pre-seismograph network event, the 1974 M w 5.4 relocated event (red-filled circle;Webb & Anderson, 1998), and the recent M w 3.7-4.7 Regional Moment Tensor focal solution events (blue-filled circles;Ristau, 2008).Inset map is a depth structure map of Base Pliocene unconformity (afterSeebeck et al., 2021).(b and c) Interpreted seismic sections of the CEF and some horizons (afterSeebeck et al., 2020).See (a) for location of sections.

Figure 3 .
Figure 3. Different approximations of the Cape Egmont Fault (CEF) structure used to examine the impact of fault geometry on the synthetic earthquake catalogs derived from RSQSim simulator.
Figure 3) Used to Examine the Impact of Fault Geometry on the Synthetic Earthquake Catalogs Derived From the RSQSim Simulator synthetic catalogs of approximately equal cumulative seismic moments for both types of slip-rate model.The slip-direction along the fault elements is also pre-defined as pure normal faulting for all realizations.

Figure 4 .
Figure 4. Graphical illustration of constant (a) and variable (b) slip rate distributions along the "simplified main fault surface" of Cape Egmont Fault (CEF) (Figure 3b).The corresponding slip rate distributions for the other three approximations of the CEF structure are shown in Figures S1-S3 in Supporting Information S1. ).

Figure 5 .
Figure5.The total of 60 RSQSim simulations performed using the "realistic fault geometry" with variable slip rates (FigureS1bin Supporting Information S1) and different combinations of the rate-and-state friction parameters and initial stress conditions listed in Table2.The earthquake catalog that best fits the natural seismicity and its parameters were applied for the rest of the analysis is shown with red color.(a) Frequency-magnitude distributions.(b) Magnitude-rupture area scaling plot.Solid black line shows a generic scaling M = log 10 A + C with C value of 4.2 (Gerstenberger et al., 2022; Stirling et al., 2021).As reference C values, dashed lines show values of 4.1 and 4.3, and dotted lines show values of 3.7 and 4.7 respectively.(c) Moment magnitude depth distribution.(d) Density of events depth distribution.

Figure 6 .
Figure 6.Frequency-magnitude distributions for the natural seismicity in the area shown in Figure 1a (Eberhart-Phillips & Reyners, 2022) and the best-fit synthetic catalog produced by RSQSim using the "realistic fault geometry" with variable slip rates (Figure S1b in Supporting Information S1), rate-and-state friction parameters a = 0.001 and b = 0.004, a steady-state coefficient of friction μ 0 = 0.6, an initial shear stress of 60 MPa, and an initial normal stress of 120 MPa.

Figure 7 .
Figure7.Synthetic catalogs produced by RSQSim that span 500 kyr and the corresponding depth distribution of the synthetic seismicity for constant (a) and variable (b) slip rates along the "simplified main fault surface" of the Cape Egmont Fault (CEF) (e.g., Figures3b and 4).The corresponding synthetic catalogs for the other three approximations of the CEF structure are shown in Figures S4-S6 in Supporting Information S1.

Figure 8 .
Figure8.Maximum earthquake magnitude initiated at each triangular element (left) and frequency of earthquakes (right) for constant (top) and variable (bottom) slip rates along the "simplified main fault surface" of Cape Egmont Fault (CEF) (e.g., Figures3b and 4).The corresponding distributions for the other three approximations of the CEF structure are shown in Figures S7-S9 in Supporting Information S1.

Figure 9 .
Figure 9. Depth distributions (a and b) and frequency-magnitude distributions (c and d) for synthetic catalogs produced by RSQSim using different fault geometries (Figure 3) with constant and variable slip rates.

Figure 10 .
Figure 10.Ruptured area (m 2 ) versus seismic moment (N × m) for all fault geometries shown in Figure 3 with constant (a) and variable (b) slip rates.

Table 1
Input Parameters Assigned to the Different Approximations of the Cape Egmont Fault Structure (

Table 3
Characteristics of Synthetic Earthquake MFDs for Different Model Runs Figure 11.3D perspective view looking northwest of the "realistic fault geometry" of Cape Egmont Fault (Figure