A Novel Method to Determine Probabilistic Tsunami Hazard Using a Physics‐Based Synthetic Earthquake Catalog: A New Zealand Case Study

Tsunamis can have devastating consequences for coastal communities. Yet hazards from future tsunamis are difficult to quantify due to their rarity in the instrumental record. Statistical earthquake catalogs have previously been used to quantify tsunami hazards. For the first time, we use a physics‐based synthetic earthquake catalog to assess probabilistic tsunami hazard in a local region. We analyze a 30,000‐year catalog to capture a wide range of the source complexity evident in tsunamigenic earthquakes. We produce models of land‐surface displacements and subsequent tsunamis from 2,585 earthquakes with magnitude > MW 7.0. Modeled slip of the Hikurangi and Tonga‐Kermadec subduction thrusts generated maximum wave heights at the coast of up to 28 m and these earthquakes pose the greatest tsunami hazard along New Zealand's coastline. The results provide a “proof of concept” for using earthquake simulators for probabilistic tsunami hazard models. They further present a platform upon which next‐generation probabilistic tsunami inundation models may be constructed.


Introduction
In coastal regions, quantifying the size and frequency of tsunamis is important for estimating the likelihood of devastating events, like those observed in the Indian Ocean in 2004 (e.g., Kowalik et al., 2005;Satake, 2014;Wang & Liu, 2007) and Japan in 2011 (e.g., Løvholt et al., 2012;Sugawara, 2021).Globally, around 80% of tsunamis are triggered by earthquakes (Grezio et al., 2010;Röbke & Vött, 2017).Earthquake simulators coupled with surface displacement and tsunami wave models have the potential to improve probabilistic estimation of tsunami hazard.Physics-based synthetic earthquake catalogs, based on earthquake source parameters including fault geometry, friction, slip and stress, provide a new tool to help investigate the local tsunami hazard posed to global or individual regions.While similar catalogs have been used globally to investigate seismic hazard (e.g., Rafiei et al., 2022;Shaw et al., 2018), this is the first time that they have been used for tsunami hazard analysis.

Tsunami Hazard Global Research
Following the 2004 Indian Ocean and 2011 Tōhoku tsunamis, considerable research has been undertaken to assess the tsunami hazard on a global scale, as well as for specific countries (e.g., Grezio et al., 2010;Horspool Abstract Tsunamis can have devastating consequences for coastal communities.Yet hazards from future tsunamis are difficult to quantify due to their rarity in the instrumental record.Statistical earthquake catalogs have previously been used to quantify tsunami hazards.For the first time, we use a physics-based synthetic earthquake catalog to assess probabilistic tsunami hazard in a local region.We analyze a 30,000-year catalog to capture a wide range of the source complexity evident in tsunamigenic earthquakes.We produce models of land-surface displacements and subsequent tsunamis from 2,585 earthquakes with magnitude > M W 7.0.Modeled slip of the Hikurangi and Tonga-Kermadec subduction thrusts generated maximum wave heights at the coast of up to 28 m and these earthquakes pose the greatest tsunami hazard along New Zealand's coastline.The results provide a "proof of concept" for using earthquake simulators for probabilistic tsunami hazard models.They further present a platform upon which next-generation probabilistic tsunami inundation models may be constructed. Plain Language Summary Tsunamis with devastating, wide-reaching consequences are significant hazards to nearly all of the world's coastlines that needs to be analyzed and quantified.Previous tsunami hazard studies have relied on sparse data retrieved from geological and/or probabilistic earthquake models to examine this infrequent but potentially deadly hazard.Consequently, they lack crucial information about the broad range of tsunami hazards resulting from local tsunamis.Our study represents the first time that a physics-based synthetic earthquake catalog has been used to assess probabilistic tsunami hazard.It overcomes many of the challenges of previous tsunami hazard assessments.We model the tsunamis generated from 2,585 earthquakes which have magnitudes greater than 7.0 from 30,000 years of a synthetic catalog of earthquakes nucleating on both the subduction zones and crustal faults through New Zealand.Tsunami wave heights at the coast of up to 28 m were recorded.Additional analysis indicates that modeled earthquakes on the Hikurangi and Tonga-Kermadec subduction thrusts pose the greatest hazard along many parts of New Zealand 's coastline. et al., 2014;Løvholt et al., 2022;Zamora & Babeyko, 2020).A wide variety of approaches have previously been used to generate earthquake catalogs for tsunami hazard assessment.These include: Monte Carlo simulations, Bayesian procedures, and examination of historical earthquake catalogs (e.g., Davies et al., 2022;Grezio et al., 2010;Løvholt et al., 2012).However, physics-based synthetic earthquake catalogs have yet to be used.To assess the tsunami wave height at the coast, either empirical relationships for regions of interest or full wave form modeling is undertaken, initializing the tsunami model with the earthquake deformation and modeling propagation, then assessing the wave height at a depth contour or at the coast, is used.These studies allow researchers to investigate how tsunami waves propagate from the source, where tsunami waves are amplified along various coastlines, as well as the maximum wave heights that could occur (e.g., Horspool et al., 2014;Zamora & Babeyko, 2020).This is important not only for understanding the hazard that tsunamis pose to coastal regions globally, but also for identifying which of those regions are most at risk of being impacted by large tsunami waves (e.g., Løvholt et al., 2022).
Physics-based synthetic earthquake catalogs allow us to extend our understanding of the impacts of local tsunamis, by investigating the interactions between earthquakes and the tsunamis on both spatial and temporal scales.Previous tsunami hazard studies have focused on a specific source, either a specific crustal fault or subduction zone.In contrast, these new physics-based synthetic earthquake catalogs allow us to investigate individual sources of tsunami hazard, while preserving not only the stress and strain relationships of individual faults, but also the large scale interactions between elements of the fault system.These catalogs also allow us to investigate the effect of multi-fault ruptures on tsunami generation and how this impacts our understanding of the tsunami hazard.This study shows that the development of physics-based synthetic earthquake catalogs provides a crucial step forward for not only understanding fault interactions, but how these interactions influence the tsunami hazard.

Previous Tsunami Hazard Assessments in New Zealand
In New Zealand's active seismic and tectonic setting, seismogenic tsunamis have the potential to generate a significant hazard (e.g., Fry et al., 2010;Kaiser et al., 2017).Local tsunami hazard (here defined as, tsunamis reaching New Zealand's coastline within 2 hr) remains poorly understood due to the short duration of the historical tsunami record and the incompleteness of paleotsunami data, both of which only provide limited information about potential tsunami sources.For New Zealand, numerous studies have been undertaken to analyze the potential tsunami hazard (e.g., Berryman, 2005;Lane et al., 2011Lane et al., , 2013;;Power, 2013;Power et al., 2013Power et al., , 2018;;Wang et al., 2017).These studies cover tsunamis generated from local to far-field sources, triggered by tectonic, volcanic and landslide processes.Previous studies that assessed the seismically-triggered tsunami hazard have either modeled specific earthquakes or used earthquake catalogs generated from statistical methods, estimated fault magnitude ranges, or historic and paleoearthquake catalogs, to assess the probabilistic tsunami hazard (e.g., Grezio et al., 2010;Horspool et al., 2014;Papadopoulos, 2003).Moreover, few of these studies encompass New Zealand as a whole (e.g., Power, 2013;Power et al., 2018Power et al., , 2022)), but instead most have focused only on specific high-risk regions, due to the computationally intensive nature of running tsunami models (e.g., Fraser et al., 2014;Power et al., 2012;Prasetya et al., 2011).
The most recent national scale assessment of New Zealand's tsunami hazard is the 2021 National Tsunami Hazard Model (NTHM; Power et al., 2022).The 2021 NTHM is an update of the 2013 NTHM and uses a combination of empirical and hydrodynamic modeling to investigate the tsunami hazard posed by local to distant seismogenic sources.It also provides tsunami hazard curves and hazard deaggregations for major cities around New Zealand to allow for further analysis (Power et al., 2022).In this study we investigate the tsunami hazard posed to the whole of New Zealand's coastline, from local earthquake sources generated from a physics-based earthquake simulator (Shaw et al., 2022) and then compare the tsunami hazard calculated in this study with the results of the 2021 NTHM.Having an alternative method for quantifying not only future potential seismicity but also the tsunamis that could be triggered, is important to investigating and quantifying the potential of these devastating natural disasters.Physics-based synthetic earthquake catalogs, like the one presented here, not only provide an alternative understanding of the tsunami hazard, but they also allow us to investigate the impacts of cascading natural hazards that can be triggered by seismicity.

Physics-Based Synthetic Earthquake Catalogs
Physics-based synthetic earthquake catalogs are generated by multi-cycle earthquake simulations (earthquake simulators hereafter), using a specified fault system and physical models of stress and strain accumulation and release (e.g., Console et al., 2017;Khodaverdian et al., 2016;Richards-Dinger & Dieterich, 2012;Tullis 10.1029/2023JB027207 3 of 20 et al., 2012).These earthquake simulators use various combinations of fault geometry and slip rates, bulk rheology, stresses, and fault friction to generate sequences of earthquakes in regions of interest (Tullis et al., 2012).With RSQSim (Rate and State Earthquake Simulator; Richards-Dinger & Dieterich, 2012), our chosen earthquake simulator, the model self-organizes to a statistically steady-state after an initial spin-up period.Due to the long time scales that these simulators can generate events over (100,000-year time scales), synthetic earthquake catalogs allow modeled hazards to be considered in a time average, probabilistic way (e.g., Console et al., 2017;Khodaverdian et al., 2016;Milner et al., 2021;Shaw et al., 2018).These earthquake catalogs can contribute to a more comprehensive understanding of the seismic hazard than using paleoseismic and digital records alone, due to the incompleteness of those records and the ability of simulators to sample a greater complexity of earthquake sources.So far, the use of earthquake simulators in hazard investigations has mostly focused on ground motion models; in that field, they offer an alternative to traditional statistics-based techniques resolving fault rupture in time and space, while yielding statistical results that are comparable to existing seismic hazard models (e.g., Rafiei et al., 2022;Shaw et al., 2018).Importantly for tsunami hazard studies, earthquake simulators allow consideration of: (a) better spatial resolution and variety of fault ruptures than many other techniques; and (b) the potential effects of complex multi-fault ruptures, including the combined rupture of a subduction interface with upper-plate faults.
There are at least three known limitations of RSQSim (Richards-Dinger & Dieterich, 2012) -and earthquake simulators in general -that may influence the modeled tsunami hazard.First, fault models used as simulator inputs are necessarily simplified surfaces that lack the roughness, geometric complexity, and spatial slip-rate variability of real fault systems.Consequently, there are fewer barriers to rupture propagation in models than in reality, leading to more "characteristic" behavior; that is, running an earthquake simulator on a simplified fault geometry leads to the more frequent occurrence of whole-fault ruptures relative to partial ruptures than for a more realistic geometry (Delogkos et al., 2022).Second, as fault ruptures can only occur on specified faults, these simulators do not account for "background seismicity" caused by unmapped active faults.Rafiei et al. (2022) accounted for this by combining physics-based and statistical approaches in a hybrid model.However, "background seismicity" is difficult to apply to tsunami modeling.Third, it is hard to anticipate how the simplifications of earthquake physics incorporated into earthquake simulators control the occurrence of phenomena such as multi-fault earthquakes.Such earthquakes clearly occur (e.g., Hamling et al., 2017), but there is very little empirical data about the rates of multi-fault versus single-fault ruptures, so that the multi-fault behavior observed in earthquake simulators is difficult to validate for specific faults (Page, 2021;Rodriguez Padilla et al., 2022).Despite these limitations, modeling tsunami hazard using a synthetic earthquake catalog is a worthwhile exercise because it represents an advance on previous work in several areas.It is important to consider possible effects of multi-fault ruptures even if their precise rates are uncertain, and background seismicity and characteristic rupture behavior are problems facing many probabilistic seismic and tsunami hazard models.By using these catalogs for probabilistic tsunami hazard assessment, in conjunction with suitable testing, we can better understand their performance and prospect for providing a crucial step forward for physically-based tsunami hazard assessments.
Here, we present the results of using a physics-based synthetic earthquake catalog to evaluate the local tsunami hazard in New Zealand.We first outline the methodology, present the results and discuss the calculated tsunami hazard and spatial variability.We then evaluate the length of the earthquake catalog used and its effects on the tsunami hazard calculations, as well as the use of physics-based earthquake catalogs in tsunami hazard assessments, including current limitations and future recommendations.This is a preliminary case study and "proof of concept" that demonstrates the feasibility of probabilistic tsunami hazard assessments from physics-based synthetic earthquake catalogs and points to the tantalizing prospects of developing probabilistic tsunami risk models.

Generation of the Physics-Based Synthetic Earthquake Catalog
In this study we use a catalog generated by RSQSim, which calculates slip on individual patches of the faults according to rate-and state-dependent friction regimes (Richards-Dinger & Dieterich, 2012).RSQSim takes the seismic cycle and divides it into three states for each fault patch (Richards-Dinger & Dieterich, 2012).These three states are healing, nucleation, and seismic rupture.During the healing stage, the fault patch is completely locked, and no stress is transferred to other fault patches.When the fault patch transitions to the nucleation state, a time-dependent breakdown of the fault strength occurs (Richards-Dinger & Dieterich, 2012).As the slip generated during the nucleation state is negligible compared to the slip during the seismic rupture phase, the stresses from the slip in this phase on neighboring patches are ignored (Richards-Dinger & Dieterich, 2012).When the fault patch transitions to the third state, seismic rupture, the slip speed continues to increase to a peak constant speed, before reverting to the healing state (Richards-Dinger & Dieterich, 2012).During the seismic rupture stage, stresses are transferred to neighboring fault patches in the system.To generate the earthquake catalog, the slip on individual patches is combined into earthquakes according to their proximity in space and time (Richards-Dinger & Dieterich, 2012).We use a catalog for New Zealand developed using RSQSim that is similar to that presented in Shaw et al. (2022), but which includes the Kermadec Subduction Zone to the north of New Zealand (Figure 1a) and uses a more realistic fault system (Seebeck et al., 2022).
After a short transient spin-up, the system self-organizes to a complex statistically steady state.Thus initial conditions do not play a role in the statistics of the catalog, though the particular sequence does depend on initial conditions, since it is a deterministic system.To minimize the number of free parameters, we choose uniform material properties and normal stress on the fault planes.For this catalog we use rate-and-state frictional parameters a = 0.001 and b = 0.008, along with normal stress sigma = 100 MPa.Elastic moduli are uniform throughout, with the Lame parameters set to 30 GPa.The combination of fault geometry, and loading aiming to match long-term slip rates, then interact with the frictional parameters to produce sequences of complex events.The final conditions, when an event stops, updated by loading between events, acts as initial conditions for the next event.Complex events, and complex sequences of events, then emerge in a statistically steady way.Friction parameter b is set by aiming to have large events match large event scaling relations, in particular scaling relations of magnitude with area and slip with length.Changing grid resolution requires a slight adjustment in b, so higher grade resolutions requires slightly higher b, since the smallest events initiating rupture at the scale of a single cell, are not fully resolved.This is effectively a single uniform global friction parameter adjustment for a single global target grid resolution parameter.A grid resolution targeting triangles of 2 km on a side is used.Behavior of the model has been calibrated against historical seismicity and scaling relations (Shaw et al., 2022).The statistical behavior is structurally stable to changing parameters with only slight quantitative changes in results, not qualitative changes, when parameters are changed (Erickson et al., 2023;Jiang et al., 2022;Shaw et al., 2018).
The fault system used to generate the catalog is an adapted version of the 2022 New Zealand Community Fault Model (NZCFM; Seebeck et al., 2022).The 2022 NZCFM holds a number of advantages, compared to previous fault models.These advantages include that it provides three dimensional representations of the faults, as well as in depth details of the fault kinematics, and it is the most comprehensive fault model currently available for New Zealand (Seebeck et al., 2022).The simulator approximation of this fault system used here does, however, use constant downdip dips, with no downdip curvature, though dip can vary along-strike.In this catalog, the Hikurangi Subduction Margin transitions to the Tonga-Kermadec Subduction Zone at 36°S and we model the latter as far as 25°S.The NZCFM contains 880 fault sources and, in cases where faults intersect, the fault with the lower slip rate is assumed to truncate against the fault with the faster slip rate (Shaw et al., 2022).To the south of New Zealand, the Puysegur Subduction Zone remains the same as that used by Shaw et al. (2022), where it connects to the southern extent of the Alpine Fault and extends to 50°S (Figure 1a).Over the time span of the catalog, the tectonic setting, the configuration of the Australian and Pacific Plate boundary, and the convergence rates of the plate boundary all remained fixed: the same as in Shaw et al. (2022).Like Shaw et al. (2022), the earthquake catalog is comprised of single-fault and multi-fault ruptures across the fault system.Due to scarcity of paleoseismological data in New Zealand, and the shortness of both the historical and digital seismic records, it is difficult to compare the seismicity generated by the simulator to these records.Instead, we have checked consistency with available statistical data, including catalog productivity and scaling relations (Shaw et al., 2022).
For the tsunami modeling, we discard the first ∼30,000 years to remove transient effects present early in the catalog, and initially take all the events greater than M W 7.0 that occur between 30,000 years and 60,000 years of the simulation time.The transient spin-up only lasts a few thousand years, but since the tsunami calculations are expensive computationally, we only use a part of the catalog, and choose to remove an initial part of the catalog comparable to the length of the catalog used.We then filtered the catalog by the earthquakes' tsunami potential energy to remove small earthquakes, which do not produce hazardous tsunamis, following the methods of Saito ( 2019) and Titov et al. (2017).The tsunami potential energy of an earthquake is the potential energy generated by displacing the sea surface from mean sea level, which then propagates horizontally as a tsunami (Saito, 2019;Titov et al., 2017).To ensure we were identifying all earthquakes that constitute the tsunami hazard, we investigated the lower bound of tsunami potential energy to use as a cut off.We identified that events with tsunami potential energies greater than 1 × 10 11 J potentially contribute to the tsunami hazard.The exclusion of events with lower tsunami energies does not significantly alter the calculated tsunami hazard.The final dataset that we used to analyze the tsunami hazard contains 2,585 earthquakes, and has a magnitude range of M W 7.0 to M W 9.25 (Figures 1b and 1c).2,024 earthquakes were excluded from the final dataset because their calculated tsunami potential energy was below the threshold.

Creation of the Earthquake Surface Deformation Models
Surface deformation models were constructed for all the earthquakes in the dataset.The vertical component of surface deformation due to each earthquake was calculated separately, following the method of Nikkhoo and Walter (2015), optimized for GPU by Thompson (2021).This method assumes an infinite elastic half space with no topography, so that we evaluated surface deformation at a constant elevation of 0 km.We assumed a Poisson's ratio (υ) of 0.25.To reduce computation time, we only calculated surface deformation at sites closer than 100 km horizontal distance to the nearest triangular fault patch that ruptured in a given earthquake.Beyond that, we assumed vertical surface deformation to be zero.To cover the entire earthquake generating domain, the earthquake deformation models were generated on a 5 km resolution mesh.For a subset of earthquakes in the earthquake catalog generated in Shaw et al. (2022), we ran the tsunami simulations using deformation models that were generated on a 2 km grid.For similar mechanisms and magnitudes of events no differences were observed between using the 2 km earthquake deformation models and the 5 km in the tsunami models.By generating the earthquake deformation models on a 5 km grid it allowed us to save computational power and memory.These deformation models were then used to initialize the tsunami propagation models.Some of the earthquakes used in the modeling process include events that rupture the sea floor, allowing us to investigate the differences resulting from events of similar magnitude that do and do not rupture the seafloor.

Running the Tsunami Simulations
The Cornell Multi-grid Coupled Tsunami model (COMCOT; Wang & Power, 2011) is used in this study to model tsunami generation and propagation.COMCOT uses a nested grid system to model the tsunami from a low resolution grid to a high resolution grid (Wang & Power, 2011).Three nested grids were used in this study.The lowest resolution grid had a spacing of ∼5 km (area covered 150°E−200°E and 20°S-60°S), while the highest resolution grid had a spacing of ∼0.6 km (area covered 166°E−179°E and 34°S-48°S).Testing was undertaken to investigate the grid resolution which captured the details of the tsunami waves while optimizing the run time of the simulations.A grid resolution of ∼0.6 km allowed us to capture the major details of the tsunami wave interactions at the coast but minimized the time taken to run the simulations.We acknowledge that running the tsunami simulations at a higher resolution would better capture the interactions of the tsunami waves along sections of the coast with complex structure and topography.However, due to this being a case study to investigate if physics-based synthetic earthquake catalogs can be used to undertake a probabilistic tsunami hazard assessment, we have chosen to optimize the computational requirements.The earthquake deformation models discussed previously were used directly by COMCOT to generate the initial sea surface displacement triggered by the earthquakes displacing the seafloor.Following generation, the propagation of the waves was simulated using linear shallow water equations in the lowest resolution grid and nonlinear shallow water equations are used in the highest resolution grid (Wang & Power, 2011).Likewise, for the lower resolution grid, bottom friction is not applied, but for the highest resolution grid, bottom friction is applied and a Manning's coefficient (a representation of the roughness applied to base of the flow) of 0.013 is used.In COMCOT a coastline grid is set up and the total depth of the water along this boundary is recorded to determine the wave height at the coast (Wang & Power, 2011).To generate this grid a 10 arc-second resolution topography/bathymetry grid was used.We first tested the required propagation time and determined that simulating 16 hr of propagation ensured that the resulting maximum wave height at the coast was always captured over all of New Zealand.The initial sea surface displacements and the maximum wave height were recorded for analysis.

Analysis of the Simulation Results
To analyze the maximum wave height at the coast, we identify the maximum wave height through time, captured in the highest resolution grid, then isolated the grid cells along the coast, by identifying the cells in the grid that had a dry neighbor and a wet neighbor.For each of the events, the spatial maximum was assessed by identifying the location along the coast where the maximum wave height was recorded.We will refer to this as the overall maximum wave height for that event.These wave heights at the coast were then used to assess the numerical and spatial distributions of the overall maximum wave heights at the coast.We recognize that if higher resolution grids had been used, larger overall maximum wave heights at the coast would likely be observed.To analyze the tsunami hazard (i.e., the expected maximum wave height at a given location for a given return period), over varying return periods, we sorted the maximum wave heights for each coastal point from smallest to largest.From this we can calculate the tsunami hazard curve as a function of return period for that point.We can then calculate the expected maximum wave height for given return period at each coastal point.As the catalog represents 30,000-year of synthetic tsunami history, the 12th highest maximum wave height at each coastal point is an estimate of the 2,500-year return period hazard (i.e., the height that can be expected to be exceeded at least once in 2,500 years on average).Likewise, the 30th highest maximum wave height is the 1,000 years return period hazard, the 60th is the 500 years, and the 300th the 100 years.These coastal values were then used to investigate how the tsunami hazard varies across New Zealand.
Further, to understand better the impact of different faults on tsunami hazard we identified earthquakes nucleating on the Puysegur Subduction Zone, Hikurangi Subduction Margin, Tonga-Kermadec Subduction Zone, and New Zealand crustal faults.Subduction and crustal faults could have ruptured in the same earthquake, and in these cases, we assign the earthquake to a fault group based on where it nucleated (i.e., an earthquake that nucleated on the Hikurangi Subduction thrust and propagated onto crustal faults would be assigned to the Hikurangi Subduction thrust group).Once we divide the events based on their fault origin, we then proceed to analyze the results using the methods previously discussed.
We investigated the spatial variability of the tsunami hazard over four regions with high hazard around New Zealand (eastern Northland, Hawkes Bay, southeast North Island and southwest South Island).For each region we calculated the tenth, 50th and 90th percentile of the expected maximum wave height over all the coastal points in that region.These were then sorted and plotted against return periods in the same way as we did for individual coastal points to give a regional tsunami hazard curve.
Finally, to justify the length of synthetic earthquake catalog, we recalculated the 2,500-year return period hazard for different time subsets of the data.Taking the 30,000-year dataset, we split the dataset into two 15,000-year subsets and three 10,000-year subsets and recalculate the 2,500-year return period hazard for each of the subsets.
We consider the differences between the recalculated hazard maps and the 30,000-year dataset to see how the hazard varies over different subsets of the data.

Maximum Tsunami Wave Height at the Coast
The largest earthquake in the dataset is a M W 9.25 earthquake that ruptures the subduction zone from the north of the Hikurangi Subduction Margin, along the Kermadec Subduction Zone, terminating at 24°S (Figure 2a).This earthquake produces tsunami waves at the coast of up to ∼19m (Figure 2b).We compare the largest magnitude earthquake to a M W 9.22 earthquake which ruptures a similar region to investigate how variability in the seafloor deformation affects the resulting maximum tsunami waves at the coast.While the maximum tsunami wave heights are generally similar around the coast, due to the larger amount of deformation around 37°S, the wave heights recorded around the northeast coast of the North Island are slightly larger for the M W 9.25 event compared to the M W 9.22 event.This illustrates how the underlying differences in earthquake slip from the physical model translate into differences in the tsunami models.
When investigating the overall maximum wave heights at the coast of the modeled tsunamis, even with our threshold some earthquakes generated negligible overall maximum tsunami wave heights at the coast (eight events have overall maximum wave height less than 0.2 m), while the largest magnitude earthquake did not trigger the largest maximum wave height at the coast (28 m), due to the distance from the earthquake epicenter to New Zealand's coast (Figure 2a).The largest wave height at the coast was triggered by an earthquake which has a magnitude of M W 9.13 and an epicentral location of 34.32°S 181.06°E.
Of all the tsunamis analyzed, 73% (1,895 events) generated an overall maximum wave height of at least 1 m, while 15% (390 events) were at least 5 m, 3% (78 events) were at least 10 m, and 2% (52 events) were at least 15 m.From these results, we would expect that, on average, a tsunami would result in an overall maximum wave height at the coast of at least 5 m once every 77 years, and a tsunami would result in an overall maximum wave height at the coast of at least 15 m once every 580 years.We plotted the spatial distribution of the overall 8 of 20 maximum wave heights, to investigate where they were concentrated along New Zealand's coastline (Figure 3b).The spatial distributions of the overall maximum wave heights at the coast allow us to investigate which regions are most likely to be impacted by the largest waves from individual tsunamis, however, many tsunamis affect extensive regions of New Zealand's coast, as can be seen in Figures 2b and 2d.

Local Tsunami Hazard at the Coast
In addition to investigating the overall maximum wave height at the coast, we assessed the tsunami hazard at the coast for four return periods: 2,500-year, 1,000-year, 500-year, and 100-year (Figure 4).Here, we define return period hazard as the probability of the tsunami wave height at the coast of a given size being exceeded at least once in a given number a years.For 30% of the North Island's coastline and 17% of New Zealand's total coastline, the 2,500-year return period hazard is greater than 5 m (Figure 4a, yellow sections of the coast).The 1,000-year return period shows similar patterns around the North Island's coastline; however, the wave heights are smaller than the 2,500-year return period.For the 1,000-year return period hazard, for 34% of the North Island's coastline and 31% of New Zealand's total coastline, wave heights are greater than 2 m (Figure 4b, green and yellow regions).For 52% of New Zealand's total coastline the 500-year return period hazard is greater than 1 m (Figure 4c, green and yellow regions), and for 25% of New Zealand's total coastline the 100-year return period hazard is greater than 0.5 m (Figure 4d).
Figure 5 shows the comparison between the 2,500-year return period tsunami hazard posed by earthquakes in the entire catalog versus the 2,500-year return period tsunami hazard posed by earthquakes nucleating on the Tonga-Kermadec Subduction Zone, the Hikurangi Subduction Margin, and the Puysegur Subduction Zone.Further discussion of the tsunami hazard posed by the subduction zones around New Zealand is included in Section 4.1.3.

Tsunami Hazard Posed by Individual Events
We have conducted a preliminary local-source hazard assessment for New Zealand's entire coastline, using a catalog of earthquakes generated from a physics-based earthquake simulator as the triggering sources.When analyzing the coastal locations of the overall maximum wave height for each tsunami modeled, the highest concentration is along the southwest coast of the South Island between 43.9°S and 46.7°S.Thirty percent of the tsunamis have overall maximum wave heights located here (Figure 3b).This is due to the proximity of the Puysegur Subduction Zone to the south of New Zealand.We also observe that ∼20% of the tsunamis have overall maximum wave heights that occur along the east coast of the North Island between 37.44°S and 41.96°S (Figure 3b).This is due to the proximity of the Hikurangi Subduction Margin, as well as the offshore extensions of crustal faults throughout the region.Some tsunami waves generated from earthquakes nucleating along the Tonga-Kermadec Subduction Zone are amplified along sections of the east coast of the North Island.Thirteen percent of the overall maximum wave heights are located along the northeast coast of the South Island between 41.4°S and 43.14°S, which is again in part due to the proximity to the Hikurangi Subduction Margin, as well as to the numerous crustal faults around the upper South Island and lower North Island, and throughout Cook Strait (Figure 1).In contrast, no individual tsunami recorded its overall maximum wave height at the coast along the north-western coast of the North Island, nor along the central east coast of the South Island (Figure 3b).This is due to there being no tsunamigenic earthquakes occurring in those regions in the catalog during the chosen time span, despite there being faults in the 2022 NZCFM that extend offshore in those regions.If a longer time span of the RSQSim earthquake catalog were to be examined, tsunamigenic earthquakes in those regions may occur.Despite no overall maximum wave heights being recorded in these regions, they can be affected by tsunami waves, as discussed in the following sections.

Local Tsunami Hazard at Varying Return Periods
We analyze both the long-term and short-term tsunami hazard posed by tsunamis.The long-term tsunami hazard (2,500-year and 1,000-year return periods; Figures 4a and 4b), captures the tsunami hazard posed by larger earthquakes, which occur less frequently in the earthquake catalog.The north and eastern sections of the North Island are most at risk of experiencing the largest tsunami waves at the coast.The scale of the expected maximum wave heights in these regions is > 8 m (Figures 4a and 4b).This is due to the triggering mechanism for these events being dominated by large magnitude earthquake sources that nucleate along both the Hikurangi Subduction Margin and the Tonga-Kermadec Subduction Zone.These are not the only regions that are affected by waves of this scale.The south-west and north-east coasts of the South Island are also locations where waves >5 m are observed.In these locations, the dominant tsunami sources are earthquakes nucleating on the Puysegur Subduction Zone and the southern Hikurangi Subduction Margin, respectively.While the above regions are affected by the largest tsunami waves, the remainder of New Zealand's coast is modeled to be affected by tsunami waves on the order of 1-5 m.Despite these waves being of a smaller magnitude they can still have damaging consequences to New Zealand's coastal regions.
We also analyze the short-term local tsunami hazard, as it captures the tsunami hazard posed by smaller earthquakes that occur more frequently in the earthquake catalog.The maximum wave heights at the coast calculated for shorter return periods are generally around 0.5-1 m (Figures 4c and 4d).Like the 2,500-year and 1,000-year return periods, there are similarities between the 500-year and 100-year return periods in the spatial distribution of the expected maximum wave heights.The east coast of the North Island is still the region which is most affected by tsunami waves compared to other sections of New Zealand.Expected maximum wave heights are around ∼2-5 m for the 500-year return period and around ∼0.5 m for the 100-year return period.For the 500-year return period, the expected maximum wave heights are less than 2 m, with the exception of a few locations where larger waves occur, but for the 100-year return period wave heights at the coast are less than 0.5 m.Even though the expected maximum wave heights are considerably less for the 100-year return period compared to the 500-year return period, we note that even tsunami waves on the scale of 0.5 m can result in significant damage to coastal communities, marinas, and port facilities (Borrero & Goring, 2015).

Tsunami Hazard Posed by Local Subduction Zones
Using physics-based synthetic earthquake catalogs, it is easy to split the catalog based on where earthquakes are nucleating, while preserving the interactions between the faults in the system (this is discussed further in Section 4.5).By dividing the catalog based on which fault each earthquake nucleates on, we can identify not only which regions are most likely to experience the largest tsunami waves over the 2,500-year return period, but also the source of those tsunami waves.We compare the 2,500-year return period hazard previously discussed for the entire catalog, with the 2,500-year return period hazard calculated using only events that nucleate on the Puysegur Subduction Zone, the Hikurangi Subduction Margin and the Tonga-Kermadec Subduction Zone to investigate how the hazard varies (Figure 5). Figure 5b reveals that if an earthquake nucleates on the Tonga-Kermadec Subduction Zone, then over a 2,500year return period, the northern coast of the North Island is likely to experience the largest expected maximum wave heights, which are on the scale of >10 m.However, large portions of the coast, of both the North and South Island's, could be affected by waves >1 m due to the refraction of the tsunami waves around the North Island, and the interaction of the waves with the bathymetry off the coast.In contrast, if an earthquake nucleates on the Hikurangi Subduction Margin, then the southeast coast of the North Island and the northeast coast of the South Island, are most at risk of experiencing the largest wave over a 2,500-year return period (Figure 5c).These expected maximum wave heights are also on the scale of >10 m.Large sections of the east coast of the South Island are also affected by smaller scale tsunami waves generated by Hikurangi Subduction Margin earthquakes.For the hazard posed by the Puysegur Subduction Zone, the area of coastline affected by the largest tsunami waves (>5 m) is much smaller than the areas affected by the other subduction zone discussed previously (Figure 5d).However, smaller scale tsunami waves (<1 m) are modeled along the western coast of both the North and South Island as well as the east coast of both the North and South Island, south of 41°S.
The analysis presented in Section 4.1 forms only a preliminary tsunami hazard assessment for New Zealand using a physics-based synthetic earthquake catalog for the tsunami triggering sources.However, it shows that physics-based synthetic earthquake catalogs can be used to undertake tsunami hazard assessments.Our assessment includes one component of the hazard, as more distant earthquake sources are not included, nor are other triggering mechanisms such as landslides and volcanoes.Future work combining results presented here with other studies assessing the tsunami hazard posed by distant tectonic sources, will result in a more thorough assessment of the tsunami hazard posed to New Zealand.

Analysis of Tsunami Hazard Curves
To evaluate the spatial variability in the expected maximum wave heights along sections of New Zealand's coast, we constructed tsunami hazard curves.The hazard curves characterize the spatial variability in expected maximum wave heights in the entire region of interest as a function of return period.These hazard curves are generated over the entire region, rather than for a singular grid point, highlighting that for a given return period, the wave heights throughout the region can vary significantly, further allowing us to quantify the variability along different coastal sections.The percentiles in the hazard curve figures show the range of expected wave heights for a given return period within that region.Figure 6 shows the hazard curves constructed for four regions around New Zealand where the tsunami hazard is dominated by either the Tonga-Kermadec Subduction Zone, the Hikurangi Subduction Margin, or the Puysegur Subduction Zone.We have included all the possible return periods calculated, due to the length of catalog that we are analyzing.However, we recommend caution when considering the hazard calculated for return periods greater than 5,000 years and less than 100 years (gray areas in Figure 6), as there are few events in the dataset that have long return periods, and at the shorter return periods the hazard is likely underestimated due to the exclusion of distant sources which dominate the tsunami hazard posed to New Zealand over the shorter return periods (Power et al., 2018(Power et al., , 2022)).
Of the four regions investigated, the northern North Island showed the greatest spatial variation in the tsunami wave height across the return periods, with differences in the expected maximum wave height over the longer return periods being up to 10 m (Figure 6a).This is due to significant variation in the coastal structure around the region, and there being a number of sheltered bays that protect sections of coast.There are three distinct steps in this hazard curve: at return periods less than 150-year where the expected maximum wave height is less than 0.5 ± 0.5 m, return periods between 150 and 1,500-year where the expect maximum wave height is greater than 1 ± 0.5 m, and at return periods greater than 1,500-year where the expected maximum wave height is greater than 5 ± 5 m.These steps in the hazard are likely due to the characteristic return periods of large earthquakes occurring along the Tonga-Kermadec Subduction Zone, as well as variation in the location of smaller magnitude earthquakes occurring in the region which can strongly influence the hazard in this region.
The variation in the wave heights recorded around Hawke's Bay is significantly smaller than that around the northern North Island.Here, variation in the expected maximum wave heights around this segment of coast is Due to the time covered by the earthquake catalog, the hazard at return periods less than 100-year and greater than 5,000 years is considered to be less reliable (gray areas).In (b/d/f/h) the tsunami wave heights along the coast in each of the areas were isolated to generate the tsunami hazard curves.much smaller (less than ±2 m over all return periods; Figure 6c).This is due to a simpler coastline around this region, which lacks the bays to shelter or amplify wave heights.However, the 50th percentile of the maximum expected wave heights, the scale of heights is similar between Hawke's Bay and the northern North Island, due to the influence of both long return period events nucleating on the subduction thrusts, and shorter return period events nucleating on crustal faults around the region.The pattern in the hazard curves generated for the southeast coast of the North Island, is similar to that discussed for the northern North Island (Figure 6e).At return periods greater than 1,500 years, the tsunami hazard along this section of coast is ∼10 ± 5 m.The drop in the tsunami hazard around the 1,500-year return period, is due to the return period of large magnitude events nucleating along the Hikurangi Subduction Margin.At return periods less than 1,500-years, both the wave heights at the coast and the variation in the wave heights also decrease (<4 ± 0.5 m), due to crustal faults dominating the tsunami hazard.
The tsunami hazard around the southwest South Island is notably lower compared to the other regions investigated (Figure 6g).Here, the expected maximum wave heights over all the return periods investigated are <2 ± 1 m.This is because the dominant source of the hazard is the Puysegur Subduction Zone, and the earthquake nucleating on the subduction zone have a smaller magnitude, on average compared to the Tonga-Kermadec Subduction Zone and the Hikurangi Subduction Margin (M W 7.6 compared to M W 8.2 and M W 8.0 respectively).There is also little variation in the tenth percentile of the maximum expected wave heights, which are close to zero.This is due to sections of this southwest South Island coastline being completely sheltered from seismic tsunamis.
The hazard curves calculated for the four regions discussed show how spatially variable the tsunami waves can be along sections of New Zealand's coast, over different return periods, in comparison to the maximum wave height previously discussed.While the hazard curves allow us to account for the spatial variability, in the following section we investigate the variability in the earthquake catalog, to see if this could additionally impact the tsunami hazard.

Justification of Earthquake Catalog Length
To investigate whether the length of earthquake catalog we used to calculate the tsunami hazard is suitable, the earthquake magnitude-frequency distribution for longer sections of the earthquake catalog were compared to the section of catalog we used, as we acknowledge that errors in the tsunami hazard are sensitive to the length of the catalog section used.Figure 7 shows the earthquake magnitude frequency for the entire catalog, as well as different segments of the catalog.
We also conducted a Kolmogorov-Smirnov Statistical Test (Pratt et al., 1981) on the three ∼30,000 years subsets of the data (Figures 7b-7d) to investigate if these three subsets were statistically different (Table 1).The Kolmogorov-Smirnov statistical test showed that the magnitude distribution of the earthquakes in the three segments of the entire earthquake catalog where not statistically different from one another (summary of the test is included in Table 1).
Due to the similarities between the earthquakes distributions across the three segments, there is little evidence to suggest that running tsunami simulations for longer intervals of the earthquake catalog would result in significantly different tsunami simulation results from those presented here.To complement this, we also divided the earthquake catalog presented here into smaller segments (two 15,000-year subsets and three 10,000-year subsets) to compare the tsunami hazard calculated by the shorter subsets to the tsunami hazard calculated by the 30,000-year subset.We then investigated the scale of the differences between alternative tsunami hazard estimates (Figure 8).Across all the subsets investigated, the tsunami hazard calculated using the entire dataset consistently produces a higher estimate of the hazard (Figure 8).While there are variations in the differences between the different tsunami hazard models, the maximum difference tsunami hazard calculated is ∼4.0 m.In the regions where the biggest differences were recorded, the differences are ∼15% of the wave height recorded at the coast.The largest differences between the tsunami hazard calculated using the subsets are around the eastern and northern coasts of the North Island.This is due to differences in the earthquakes that are occurring along the Tonga-Kermadec Subduction Zone and the Hikurangi Subduction Margin.There are also slight difference in the hazard recorded around the southwest South Island, which is due to differences in the earthquakes occurring along the Puysegur Subduction Zone.However, the scale of all the difference compared to the wave heights around all of New Zealand's coast is small, indicating that when any of the subsets of the wider dataset are used, similar hazard results are recorded.
Despite there being differences between the entire dataset presented here, and subsets of that dataset, as these differences are relatively small compared to the tsunami wave heights recorded at the coast, we believe that using a 30,000-year segment of the entire earthquake catalog provides a good representation of the tsunami hazard.Smaller time segments of the entire earthquake catalog are likely to underestimate the tsunami hazard, while longer segments of the entire earthquake catalog are like to capture similar patterns in the seismicity that have already been captured by our 30,000 years segment.

Comparison of the Tsunami Hazard From Synthetic Earthquake Catalogs to Those From Traditional Tsunami Hazard Models
To understand how our tsunami hazard calculated using a physics-based synthetic earthquake catalog compares to previous hazard assessments, we compare our model to a subset of the 2021 NTHM (Power et al., 2022), which is generated from more traditional tsunami hazard assessment techniques.While Power et al. (2022) considers the full seismic tsunami hazard, we have extracted only the local seismic tsunami hazard from the full model to compare it with our model.Therefore, we compare the 2,500-year return period tsunami hazard from our model (Figure 9a) with the subset of the 2021 NTHM (Figure 9b; Power et al., 2022).

Interpretation
There are no differences in the earthquake magnitudes between the two subsets of the earthquake catalog There are no differences in the earthquake magnitudes between the two subsets of the earthquake catalog There are no differences in the earthquake magnitudes between the two subsets of the earthquake catalog Note.Subset A is the time period 32,000 years to 63,000 years, Subset B is the time period 63,000 years to 93,000 years, and Subset C is the time period 93,000 years to 112,000 years.A two-sided statistical test was conducted at the 95% level of significance to investigate if the earthquake magnitudes were statistically significantly different.The p-value, the conclusion of the test (i.e., the rejection or failure to reject the null hypothesis), and the interpretation of the results are included in the table.The remaining difference between the two tsunami hazard models is that the 2021 NTHM uses a statistical model of magnitude frequency for each source and uses a modeling-informed semi-empirical method to evaluate the tsunami wave heights (Power et al., 2022), while we undertake numerical modeling for each of the source events generated using RSQSim.To account for the differences in the modeling procedure and to undertake an equivalent comparison, we compare our 2,500-year return period model to the 50th percentile 2,500-year return  c-d) Difference in the tsunami hazard calculated using the two 15,000-year subsets and the 2,500-year return period hazard calculated using the entire dataset.(e-g) 2,500-year return period hazard calculated using three 10,000-year subsets.(h-j) Difference in the tsunami hazard calculated using the three 10,000-year subsets and the 2,500-year return period hazard calculated using the entire dataset.Note.The log scale of the tsunami wave height at the coast is the same for all hazard maps, and the linear scale for the differences plots is the same.
period tsunami hazard calculated for the subset of events in the 2021 NTHM (Power et al., 2022).Overall, similar patterns in the tsunami hazard are observed in the two models.For example, for the North Island the hazard identified along the east coast is higher than the west coast, and for the South Island the hazard identified along the southwest and northeast coasts are higher than the other sections of the coastline.Along the west coast of New Zealand between 37°S and 47°S the tsunami hazard calculated by the 2021 NTHM is higher than that calculated using the synthetic earthquake catalogs (Figure 9).This is also the case along the east coast of the South Island between 43°S and 47°S.In contrast, the tsunami hazard calculated using the synthetic earthquake catalogs indicates a higher level of tsunami hazard along the northeast coast of the North Island between 34°S and 38°S, and similar levels of hazard along the east coast of the North Island between 40°S and 41°S (Figure 9).
While our methods produce similar results to more traditional methods, our methods provide an alternative approach for undertaking tsunami hazard assessments.In the next section we discuss the advantages, and limitations, of using physics-based synthetic earthquake catalogs to calculate the tsunami hazard models in comparison to more traditional methods.

Evaluation of the Methods to Determine Tsunami Hazard
The use of physics-based synthetic earthquake catalogs provides a novel method for investigating the tsunami hazard posed to vulnerable coastlines.The results presented here show a "proof of concept" that these catalogs can be successfully used for undertaking probabilistic tsunami hazard assessments.Previous investigations for New Zealand have focused on the tsunami hazard from individual sources (e.g., the Tonga-Kermadec Subduction Zone; Power et al., 2012), with few studies encompassing multiple tsunami sources (Power et al., 2022).Developing an earthquake catalog that captures the seismicity on both crustal faults and subduction zones enables us to investigate the tsunami hazard posed to New Zealand by a wider range of tectonic sources.These catalogs capture not only earthquake and aftershock sequences on individual faults, but also the spatial extent and timing of ruptures on faults, as well as the interactions between faults, and multi-fault ruptures.This further extends our ability to understand the potential tsunami hazard.The methods presented here allow us to efficiently investigate the tsunami hazard posed by the entire fault system, or by individual faults, while still preserving the stress and strain interactions between the elements that are used to construct the earthquake catalog.
In addition to isolating earthquakes occurring on specific faults, the catalog can be easily split based on the timing of the earthquakes.This allows us to investigate the tsunami hazard posed due to the temporal variations on individual faults and between faults, again still preserving the stress and strain parameters on the entire system.Combined with the advent of time-dependent earthquake forecasting, our approach also presents an opportunity to investigate time-dependent tsunami hazard forecasting.The long time period covered by this catalog provides an additional advantage, due to the wide array of potential earthquake scenarios.It enables a more comprehensive hazard assessment and generates enough source information to begin the process of developing probabilistic inundation and risk models.This is particularly important for local tsunami hazard, which can have devastating impacts and very little time between the earthquake and the arrival of the first tsunami waves at the coast.Having a wider range of potential earthquake scenarios also allows for a more in-depth assessment of how the complexities of the tsunami source influence the waves impacting the coast in the near field, which are vital to include when assessing the probabilistic tsunami hazard and risk.
One of the limitations of using the physics-based synthetic earthquake catalog presented here is that it only provides local earthquake sources, so does not include the additional hazard posed by distant seismogenic sources.Due to computational requirements, distant earthquakes from around the Pacific Ocean are not generated in the catalog.Distant tectonic sources tend to dominate the tsunami hazard posed to New Zealand over short return periods (Power et al., 2018(Power et al., , 2022)), however, the local seismically-triggering sources tend to be more devastating both in the overall size of the tsunami waves at the coast (for a given magnitude of earthquake) and the generally short time between the earthquake and the arrival of the first tsunami waves at the coast.As a result, we have chosen to focus on the local tsunami hazard.Despite this, distant earthquake sources do contribute to the overall tsunami hazard posed to New Zealand, and further work needs to be undertaken to determine the best methods to combine previous work on the hazard determined from distant sources with the local hazard assessment conducted here.This will provide a more in-depth understanding of the seismically-triggered tsunami hazard in New Zealand.
The tsunami hazard assessment discussed here is limited to the types of earthquakes generated using RSQSim.This catalog does not include the dynamic processes that lead to tsunami earthquakes, which have previously occurred in New Zealand (e.g., Bell et al., 2010).Tsunami earthquakes are those that can produce larger tsunami waves than expected for their magnitude (Polet & Kanamori, 2022).The physics behind tsunami earthquakes are not fully understood, which makes them tricky to include in current physics-based earthquake models (Polet & Kanamori, 2022).However, we note that our synthetic model attempts to account for the seismic moment release of these earthquakes, treating them more as a typical "fast" earthquake.While RSQSim can model propagating seismic waves, due to the computationally costly nature of these simulations, and that the tsunami simulations do not utilize these models, we do not use an earthquake catalog generated in this way.Earthquakes in the outer-rise and subducting slab, which are accounted for in the New Zealand 2021 NTHM, are not modeled in RSQSim, as they are not part of the fault system used to simulate earthquakes.We acknowledge that these types of earthquakes contribute to the overall tsunami hazard and need to be included in future hazard assessments.Further work on the synthetic earthquake catalog generated for New Zealand to include the fault systems where these earthquakes occur, or the development of methods to incorporate tsunami from alternative sources into the analysis, will aid in developing a more complete and realistic tsunami hazard assessment for New Zealand.Adapting the earthquake simulator, and the surface deformation models, to account for earthquakes that result in larger surface displacements due to amplification effects in soft sediment layers or complex splay-fault structures, will allow for a more comprehensive assessment of the tsunami hazard.Any biases in the RSQSim catalog will result in bias in the tsunami hazard and further research is being undertaken to quantify how well these catalogs correctly capture the full range of earthquakes.
Another limitation of using these catalogs is that the relatively characteristic earthquakes present in the catalog may not accurately reflect the range of earthquake ruptures possible, as can be seen in Figure 7a in the banded appearance of the earthquake magnitudes.This is due to the lack of roughness on the faults and the lack of variability in the stress components.As a result, we may be missing some of the variability of earthquake triggering sources in the synthetic catalog.However, more testing is required to investigate if the physics-based model provides a better estimate of the seismicity compared to previous models that use assumed magnitude-frequency scaling relationships for individual faults.Alongside this, further testing is required to determine whether the simulated tsunamis are representative of real tsunamis that are generated by earthquakes of comparable magnitude and location.For the subduction zones, specifically the Tonga-Kermadec subduction zone, the homogeneous nature of the subduction zone along strike may not accurately reflect the physical nature of the region.Due to the lack of geophysical data sampling along the Tonga-Kermadec Subduction Zone, the coupling of the subduction zone as well as the heterogeneity of the subduction zone is poorly understood.In addition, due to the short digital earthquake record in the region, it cannot be confirmed or contradicted that the subduction zone models used to generate the earthquake catalog analyzed here adequately reflects the long term properties of the subduction zones.As a result, the decision was made to used simplified models of the Hikurangi Subduction Margin and the Tonga-Kermadec Subduction Zone.Further work is needed to add variability along both the subduction zones and the faults, to encompass the full range of the faults.We also acknowledge that the dataset used to analyze the tsunami hazard here does not produce a standard Gutenberg-Richter relationship for earthquakes greater than M W 7 (Figure 1b), despite the overall dataset following an approximate Gutenberg-Richter relationship.This may in part be because we did not include the earthquakes that had low tsunami potential energy in our dataset.
While there are challenges to overcome, physics-based synthetic earthquake catalogs are a useful step forward in assessing the tsunami hazard posed to New Zealand.Future work to adapt the current catalog will provide a wider range of possible tsunami triggering scenarios, increasing the comprehensive nature of this form of tsunami hazard assessment.Moreover, by combining the local tsunami hazard present here with far-field hazard assessment, a more in-depth understanding of tsunami hazard posed to New Zealand's coastline can be determined.
The study presented here shows that physics-based synthetic earthquake catalogs can provide crucial information toward tsunami hazard analysis that can be missing when other methods are implemented.The work flow developed using this catalog, can be easily replicated using other synthetic earthquake catalogs that address the limitations we have outlined.Overcoming the limitations of the earthquake catalog, would enable the tsunami hazard to be calculated at a higher resolution so that finer details of the tsunami interactions at the coast can be captured.Additional work to investigate the sensitivity of the tsunami hazard to changes in the parameters of the earthquake simulator, will provide alternative pathways of analysis, and extend further our understanding of the tsunami hazard that can be provided by the earthquake simulators.

Conclusion
Physics-based synthetic earthquake catalogs provide a new tool for investigating tsunami hazard.We have developed a method for effectively using an earthquake catalog generated by RSQSim to conduct a preliminary tsunami hazard assessment using New Zealand as a case study.Preliminary results indicate that some earthquakes produce overall maximum tsunami wave heights at the coast of up to 28 m, and the northern and eastern coasts of the North Island are most at risk of being impacted by the largest tsunami waves, both from individual events and over 2,500-year and 1,000-year hazard return periods.Future work to improve the earthquake catalogs, and extend the hazard analysis to distant earthquake sources, will help refine the hazard analysis presented here.Further testing of the realism of the simulated earthquakes and tsunamis, as well as understanding the sensitivity of the models to input parameters and data will provide additional robustness to the hazard assessment.Extending from the work presented here, these catalogs can be used to undertake inundation modeling to underpin future generations of probabilistic tsunami modeling.

Figure 1 .
Figure 1.Physics-based synthetic earthquake catalog used for tsunami analysis.(a) Key tectonic locations through New Zealand, AUS = Australian Plate, PAC = Pacific Plate, NI = North Island, SI = South Island, AF = Alpine Fault, CS = Cook Strait, Puysegur = Puysegur Subduction Zone, Hikurangi = Hikurangi Subduction Margin, Kermadec = Tonga-Kermadec Subduction Zone.Key populations centers Auckland (a), Wellington (W) and Christchurch (c) are shown by the red letters and squares.(b) Histogram of earthquake magnitudes in the dataset.(c) Earthquake epicentral locations and magnitudes of the earthquakes in the dataset.

Figure 2 .
Figure 2. Analysis of two M W 9.2 earthquakes and the resultant tsunami.The black circle depicts the earthquake's epicenter.(a/c) Seafloor deformation/initial water surface displacements for the respective earthquakes.(b/d) The maximum wave height of the tsunami recorded on the highest resolution grid (∼0.6 km grid spacing).(e) Difference in the modeled maximum wave height of the tsunamis generated by the M W 9.25 and M W 9.22 earthquakes displayed in (b) and (d).Note the log scale used to show the maximum wave height of the tsunami, and the shortened linear scale to show the differences in the maximum tsunami wave heights.

Figure 3 .
Figure 3. Analysis of the overall maximum tsunami wave heights at the coast (i.e., the maximum both spatially and temporally for a given tsunami).(a) Histogram of overall maximum tsunami wave heights at the coast generated from the dataset.(b) Spatial distribution of the overall maximum tsunami wave height at the coast around New Zealand, generated from the dataset.Note.Each point in (b) represents where along the coast the overall maximum wave height for an individual event occurred; but for some locations multiple tsunamis register the overall maximum wave height at the coast at the same location, resulting in an overlap of data points.A log scale is used to depict the overall maximum wave heights.

Figure 4 .
Figure 4. Tsunami hazard from local earthquake sources posed to New Zealand over differing return periods.(a) 2,500-year (b) 1,000-year (c) 500-year (d) 100-year.Key population centers, Auckland, Wellington, and Christchurch, are shown by the black squares.Note.A log scale is used to show the expected maximum tsunami wave heights at the coast.

Figure 5 .
Figure 5.Comparison of the 2,500-year return period tsunami hazard posed by different earthquake sources.Maps show expected maximum wave height for: (a) All earthquakes in the dataset.(b) Earthquakes nucleating on the Tonga-Kermadec Subduction Zone.(c) Earthquakes nucleating on the Hikurangi Subduction Margin.(d) Earthquakes nucleating on the Puysegur Subduction Zone.Key population centers, Auckland, Wellington and Christchurch, are shown by the black squares.Note.A log scale has been used to display the expected maximum tsunami wave heights at the coast.

Figure 6 .
Figure 6.Tsunami hazard curves of the expected maximum wave heights for the indicated regions showing the 90th percentile, the 50th percentile and the 10th percentile hazard within each region.(a) Tsunami hazard curve calculated for the northern North Island.(b) Area used to calculate the hazard curves in (a).(c) Tsunami hazard curve calculated for the Hawke's Bay Region.(d) Area used to calculate the hazard curves in (c).(e) Tsunami hazard curve calculated for the southeast North Island.(f) Area used to calculate the hazard curves in (e).(g) Tsunami hazard curve calculated for the southwest South Island.(h) Area used to calculate the hazard curves in (g).Key population centers, Auckland, Wellington and Christchurch, are shown by the red squares.Note.Due to the time covered by the earthquake catalog, the hazard at return periods less than 100-year and greater than 5,000 years is considered to be less reliable (gray areas).In (b/d/f/h) the tsunami wave heights along the coast in each of the areas were isolated to generate the tsunami hazard curves.
the null hypothesis Fail to reject the null hypothesis Fail to reject the null hypothesis

Figure 7 .
Figure 7. Magnitude frequency distribution of all the earthquakes in the entire earthquake catalog.The teal graphs show the time frame which earthquakes in this study were identified from.(a) Earthquake magnitude through time for the entire catalog.(b)-(d) Earthquake magnitude frequency distributions for two ∼30,000-year blocks and one ∼20,000-year block of time within the synthetic earthquake catalog.Histograms shows the counts of earthquakes within the magnitude band, and the line plots show the cumulative frequency of events greater than a given magnitude.

Figure 8 .
Figure8.Comparison of the 2,500-year return period hazard calculated for different subsets of the dataset used in this analysis.(a-b) 2,500-year return period hazard calculated using two 15,000-year subsets.(c-d) Difference in the tsunami hazard calculated using the two 15,000-year subsets and the 2,500-year return period hazard calculated using the entire dataset.(e-g) 2,500-year return period hazard calculated using three 10,000-year subsets.(h-j) Difference in the tsunami hazard calculated using the three 10,000-year subsets and the 2,500-year return period hazard calculated using the entire dataset.Note.The log scale of the tsunami wave height at the coast is the same for all hazard maps, and the linear scale for the differences plots is the same.

Figure 9 .
Figure 9.Comparison between tsunami hazard calculated by different models.(a) 2,500-year return period tsunami hazard calculated in this study.(b) 50th percentile 2,500-year return period tsunami hazard calculated using a subset of the 2021 NTHM sources, containing only local earthquake sources, which are equivalent to the sources used to calculate the tsunami hazard in this study(Power et al., 2022).(c) Differences between the 2,500-year return period hazard calculated using our physically-based model (a) and the 50th percentile 2,500-year return period hazard calculated using a subset of the 2021 NTHM (b;Power et al., 2022).

Table 1
Kolmogorov-Smirnov Statistical TestRun on the Earthquake Magnitude in Three Subsets of the Entire Earthquake Catalog