Underground hydrogen storage (UHS) is a promising solution for large-scale, long-term energy storage by converting surplus renewable electricity into hydrogen for later use. Hydrogen recovery efficiency is an essential indicator for evaluating the performance of UHS. However, the impacts of permeability anisotropy, including horizontal permeability anisotropy (kh,max/kh,min) and horizontal-to-vertical permeability anisotropy (kh/kv), on hydrogen recovery efficiency have not been thoroughly investigated. This study addresses this gap with numerical simulations that couple a simplified box model and a field-realistic model based on the Ahuroa gas storage site. Our results show that kh,max/kh,min affects hydrogen recovery efficiency, and this effect becomes more pronounced at higher permeability. The difference in hydrogen recovery efficiency between kh,max/kh,min=1 and 10 expands from 1.3 percentage points at 100 mD to 4.3 percentage points at 200 mD. Increasing kh/kv can reduce the effect of kh,max/kh,min. In reservoirs with horizontal anisotropy, orienting a horizontal well along the minimum rather than the maximum permeability direction can reduce hydrogen recovery efficiency by up to 8 percentage points (from 78 % to 70 %), with all other conditions held constant. These findings from the box model were subsequently validated through simulations of a realistic geological site (Ahuroa natural gas storage site). This study indicates that permeability anisotropy is an important consideration in UHS site assessment. Ignoring it can result in errors in hydrogen recovery efficiency evaluation, leading to inaccurate predictions of recoverable hydrogen volumes.

With the rapid expansion of renewable energy deployment, underground hydrogen storage (UHS) has emerged as a promising large-scale storage option to smooth seasonal fluctuations in electricity supply. However, current geological assessments for UHS primarily focus on reservoir properties such as porosity and permeability, while overlooking the influence of interlayers. To address this research gap, this study systematically investigates how interlayer characteristics, including permeability, thickness, and geometry, affect hydrogen recovery efficiency and the extent of unrecoverable hydrogen. Numerical simulations were conducted using OpenGoSim (OGS) software, beginning with a box model and extending to a realistic geological model based on the Ahuroa gas storage site.Simulation results reveal that, under fully perforated conditions, lower interlayer permeability impedes upward hydrogen migration, thereby reducing the impact of gravity override and enhancing hydrogen recovery efficiency. With decreasing interlayer permeability, hydrogen transport into the interlayer shifts from advection-dominated transport to diffusion-dominated transport, resulting in greater hydrogen volume than predicted by a model that neglects molecular diffusion. The study further demonstrates that optimizing well configuration should consider interlayer properties to maximize recovery efficiency. Consistency between the box model and the realistic geological model supports the generality of the findings.By integrating interlayer characteristics into site evaluation, this work enhances the accuracy of hydrogen recovery efficiency evaluation and advances UHS assessment.

Faults in the brittle upper crust are thought to primarily grow due to repeated earthquakes. To understand better fault growth during incremental slip events we analyse geometric and displacement data for timescales of individual earthquakes (since 1840 AD) to millions of years on New Zealand active faults. The active faults studied are from connected networks with a range of orientations, lengths (1–200 km), displacement rates (0.1–27 mm/yr) and slip types. Our data indicate that individual earthquakes produce slip on multiple faults with variable sizes, orientations and slip types; in some cases these earthquakes cross tectonic domain boundaries. Earthquakes generally produce partial rupture of reactivated bedrock faults and show little evidence of tip propagation, characteristics most closely resembling the constant-length fault growth model, with growth primarily achieved by increases in cumulative slip. Earthquake slip profiles display a range of shapes with one or more maxima. High gradients along faults and approaching fault tips reflect coseismic slip transfer to nearby faults. These high slip gradients are consistent with stress interactions and kinematic coherence between faults during individual earthquakes (i.e., timescales of seconds to minutes). Coseismic increments of slip increase with rupture length and are described by a power function of ~0.5, while the power function for cumulative displacement and final length is ≥1. An important consequence of these divergent power functions is that larger faults broadly accrue their finite displacements in more earthquakes than smaller faults. The increase in earthquake number with fault size is achieved by a combination of shorter recurrence intervals and longer faulting histories for larger faults. We believe that our New Zealand observations have global application.

Fault geometry and the connectivity between faults at depth are both important controls on the nucleation, propagation and arrest of earthquake rupture, so modelling these parameters accurately is essential to models of the earthquake cycle. However, simulations involving complex three-dimensional (3D) fault systems rarely explore the sensitivity of results to uncertainties in geometry and connectivity — either in terms of modelled earthquake characteristics or impacts such as ground shaking and surface deformation. In many cases, geometry-related sensitivity testing is limited because it is challenging to construct a suite of alternative fault models that span the range of plausible fault geometries, intersections and connections; such alternative models are especially difficult to construct for systems where faults truncate or cross-cut each other at depth. We present a new, semi-automated method that simplifies creation of 3D models of networks of tens or hundreds of faults, combining open-source python tools with the meshing capabilities of Leapfrog(TM) software. The new workflow reduces the time to create a fault model of 113 faults in central Aotearoa New Zealand by ~80%, from 25 hours to 5 hours of human input. This improvement significantly decreases the effort required to create multiple alternative fault geometries, making detailed sensitivity analyses more feasible. The applicability of the workflow is demonstrated for the creation of three alternative models of fault geometries for central Aotearoa New Zealand.

Fault-seal algorithms were primarily developed from outcrops where clay-rich fault-rock is mainly derived from shale beds. We analyse the continuity of shale-smears and consider their impact on fault seal for interbedded sand-shale sequences in New Zealand. Our data are from six coastal outcrops of poorly lithified strata (burial depths similar to <= 1.5 km) displaced by a random sample of 194 small normal faults with displacements of 0.02-1.4 m. Each fault is 100% exposed in cross section and displaces a single shale bed by more than the bed thickness. The faulted shale beds display a range of geometries from no smear, to discontinuous smear and continuous smear. The relative frequencies of these three smear types vary between the six outcrop localities, with up to 60% of shale beds at each locality showing no smear and discontinuous smears typically covering <50% of the fault-trace length between shale-bed cutoffs. First-order changes in shale-smear continuity between sample localities reflect differences in shale composition and competence, while locally the geometries of individual smears can be controlled by the number and displacements of slip surfaces within fault-zones. The absence of shale-smear on many beds decreases fault-seal potential and could be accounted for in shale-smear algorithms.

Tsunamis have the potential to cause catastrophic damage to coastal communities. In Aotearoa New Zealand, where 3.5 million people reside within 5 km of the coast, the threat of experiencing a tsunami within their lifetime is a stark reality. Although these events are infrequent with recurrence intervals of hundreds of years, New Zealand faces an elevated risk due to its location within the tectonically active Pacific, where over 80% of the world's tsunamis occur. The region has experienced over seventy tsunamis in the past two-hundred years, with five of these causing devastating impacts to coastal communities and leaving an indelible mark on the landscape due to wave amplitudes surpassing 5 m at the coast. Recent studies, while crucial, have predominantly focused on assessing the tsunami hazard from local sources, recognising their immediate threat. However, to comprehensively assess the overall tsunami hazard to Aotearoa, we must fully account for the regional and distant sources also. This is informed by the harsh reality that some events, such as the 1877 Northern Chile and the 2004 Indian Ocean tsunamis, have inflicted staggering death tolls in distant locations, emphasising their paramount significance in our hazard assessment efforts. In this talk, I will present our innovative hybrid tsunami hazard model designed for Aotearoa New Zealand. We use observations of accumulated earthquake slip on active faults in the Pacific alongside established earthquake laws to ensure that we capture a wide variability of seismogenic tsunami sources to complement the limited historical and instrumental records. Due to recent computational advancements, we can now calculate the seafloor deformation generated from hundreds of synthetic tsunami sources across twenty subduction zones and simulate the tsunami wave propagation to the coast of New Zealand. For each source, we can estimate the wave amplitudes and timing of potential tsunamis and use these metrics to calculate the hazard that these regional and distant sources pose over common return periods. Each part of the model, from the source characteristics to the wave propagation has been independently tested and benchmarked with recorded events to ensure the rigor of the research. Our hybrid approach of blending observation-driven, physics-based, and probabilistic methodologies offers a comprehensive approach to assessing the full range of earthquakes that could cause a tsunami at the shores of New Zealand. Our work, alongside the recent research carried out on the local tsunami sources will accelerate Aotearoa New Zealand’s natural hazard resilience from Pacific earthquake-generated tsunami sources and will pave the way for other tsunami mechanisms to be incorporated into the model analysis, an urgent need given that these hazardous events do not occur independently. We look forward to having the opportunity to share our Aotearoa New Zealand tsunami hazard model with the wider tsunami community in Europe and discuss pathways that our combined research could follow to help build safer and more resilient coastal communities, globally.

Large ( >= Mw 6.5) earthquakes recorded in active fault systems are commonly clustered in space and time, which presents challenges for time-dependent seismic hazard modeling. We investigate the spatial and temporal clustering of earthquakes in the last 5500 yr on upperplate faults (Wairarapa, Wellington, and & Omacr;h & amacr;riu) and the subduction interface in the southern Hikurangi margin in Aotearoa-New Zealand. We recalibrated radiocarbon ages and reinterpreted some earthquake timing interpretations from 37 on-land sites (trenches) to produce revised earthquake timings and recurrence intervals on three upper-plate faults. We compare these ages with the timings of great earthquakes ( >= Mw 8) on the Hikurangi subduction interface and the 1848 Marlborough and 1855 Wairarapa historical surface-rupturing earthquakes. Temporally clustered surface-rupturing earthquakes occurred on two or more upper-plate faults at 270-90, 880-520, 2300-1825, 3640-2810, and 5170-4855 cal. B.P. The youngest four of these earthquakes overlap in age with the timing of ruptures on the southern Hikurangi subduction interface. A further two subduction interface earthquakes at 515-475 and 1505-1250 cal. B.P. do not temporally overlap with the upper-plate earthquakes studied. Over half of the earthquakes sampled on the subduction interface are clustered in time with upper-plate earthquakes on the Wairarapa, Wellington and/or & Omacr;h & amacr;riu faults. The observed spatial and temporal clustering of large earthquakes could reflect co-rupture of multiple faults and/or sequences of earthquakes closely spaced in time. The clustering is consistent with geometric intersection and/or stress interactions between upper-plate faults and the subduction interface.

Underground hydrogen storage (UHS) in depleted gas reservoirs is a possible solution for large-scale seasonal energy storage. A key challenge with UHS lies in an elevated gas-water contact zone due to prior hydrocarbon exploitation. This condition increases the risk of water invasion during hydrogen injection and withdrawal operations, adversely affecting overall storage performance. However, the mechanisms controlling this phenomenon and effective mitigation strategies remain poorly understood.This study aims to optimize UHS operational strategies by investigating the impact of operational parameters on water invasion behaviour using a model based on an existing underground gas storage (UGS) facility. In this study, numerical simulations were conducted to evaluate the influence of various operational parameters, including cushion gas injection rate and volume, along with hydrogen injection and withdrawal rates, durations, and volumes. Our results demonstrate several key findings. First, hydrogen's high compressibility reduces water invasion effects compared to natural gas when injecting and withdrawing equivalent volumes. Second, extending the duration of withdrawal stage and increasing withdrawal rate significantly increases water production rate. Third, cushion gas volume strongly influences reservoir pressure, with larger volumes reducing water production due to a stabilized pressure baseline. Fourth, the amount of gas injected and withdrawn in each cycle positively correlates with water production rate, while increasing injection and withdrawal frequency mitigates water invasion effects. Finally, cushion gas and hydrogen injection rates have less impact on water invasion.These findings provide a basis for optimizing UHS design parameters in depleted reservoirs by specifying injection and withdrawal schemes and cushion gas management, intending to mitigate water invasion risks.

In outcrops, the hanging-wall and/or footwall structure around a fault are often exposed, while the underlying fault is poorly resolved. In these cases, it is desirable to estimate the location and shape of the fault at depth, especially if it belongs to an active fault system prone to large earthquakes. The Mw 7.8 Kaikōura earthquake occurred two minutes after midnight on 14th November 2016, causing at least 17 faults in the northeast South Island of New Zealand to rupture, including a number of faults that had not been previously mapped. One of these smaller new faults is the Leader Fault, which at the surface displaces Mesozoic interbedded greywacke and argillite. In outcrop, the fault rupture caused an over 3 m high, 20-30 m wide, and over 120 m long hanging-wall fold to appear at the surface.In September 2022, we used a differential global navigation satellite system to map the topography of the fold. We collected a total of 1493 points over a map area of 4526 m², i.e. an average point density of ca. 1 point per 3 m². The data were meshed into a three-dimensional triangular surface, which was then sectioned into ten cross-sections, each 10 m apart and perpendicular to the fold axes. We present fault-prediction modelling of two of these sections. In the Movetm software (Petroleum Experts), we used two methods of fault prediction; constant heave and constant slip. Both methods require implicit information about the hanging-wall shape, the position of the fault at the surface and the “regional”, i.e. the position of the hanging wall before deformation. Before the modelling, all this information was known apriori; i.e. we mapped the shape of the ground surface, we knew the fault to outcrop at the break of slope at the front of the leading edge, and the regional is an extension of the undeformed footwall. Both modelling techniques require a seed, i.e., a small portion of fault at the surface with a certain angle of dip. We use a horizontal and a 60° dipping seed.We can estimate the fault geometry down to a depth of 20-25 m. For both sections, we predict the fault is steep, greater than 60°. Using a flat seed gives a slightly listric fault geometry, but in any case, the fault is steep down to 20 m depth before flattening out slightly. Compared to a small (15 cm) outcrop of the fault plane (dipping 75° WNW) at the surface at the northern end of the outcrop, the best matches are given by modelling with constant slip. The steep fault geometry is governed by the basement rock that has steep bedding that also dips ca. 70° WNW.

The timing and size of successive prehistoric earthquakes on individual active faults are key for understanding seismic processes and time‐dependent seismic hazards. Here, we analyze interevent and elapsed times for 890 large prehistoric and historic earthquakes on 210 normal, reverse and strike‐slip faults from five active tectonic regions globally (Japan, Greece, New Zealand and the California & Basin‐and‐Range provinces in the US). Most faults (∼80%) have mean interevent times greater than the elapsed time (open‐interval) since their last recorded earthquake. We also find that 85%–100% of closed interevent times, defined by 64 historic ruptures and their penultimate events on these faults, occurred within a factor of two of their mean recurrence‐interval, with 75% less than the mean. These observations hold for a variety of tectonic settings and fault parameters, with faster slip‐rate faults (>10 mm/a) being consistently more “advanced” in their seismic‐cycle than slower moving faults. The entire global population of closed interevent‐times is consistent with a Weibull probability density function (PDF), while stochastic modeling tailored to closed‐interval parameters indicates that open recurrence‐interval data sets are best “predicted” by positively skewed elapsed time distributions (52%–78% overlap integral) for all regions, except California. Thus, the rarity of elapsed times exceeding mean interevent‐times on individual faults may be due to skewed recurrence PDFs (i.e., Brownian Passage Time, lognormal, etc.), in which the median and mode are less than its mean, while California is an outlier potentially because its open‐intervals derive from a single geometrically interconnected fast‐moving (>10 mm/a) fault system that is presently experiencing an earthquake‐hiatus.

Greece is Europe’s most seismically active country, as it is being deformed by an active subduction-system and one of the world’s fastest-spreading continental rifts. Onshore active faults pose seismic-hazard that cannot be reliably assessed in the absence of a comprehensive map of potential earthquake sources. Here, we use high-resolution Digital Elevation Models (DEMs), in conjunction with hillshades and slope-models, to map and characterise faults in Greece at a scale of 1:25000. The Active Faults Greece (AFG) database records 3815 fault-traces assigned to 892 interpreted faults. Of these traces, 53% were mapped here for the first time, with their geometries and slip-sense constrained by displacement of landscape features. AFG includes >2000 active and 1632 probably-active traces, while 35 traces result from historic surface-ruptures. Many faults (57%) exhibit strong depositional-control (DC) on sedimentation patterns, with active faults featuring approximately equal numbers of sharp (32%), moderate (29%) and rounded (29%) scarps. AFG is the first fault database in Greece generated using nationwide interpretation of geomorphology, with applications in paleoseismology, seismic-hazard-assessment, mineral-resources exploration and resilience-planning.

Underground storage of green hydrogen in depleted gas fields could provide Aotearoa New Zealand (ANZ) with a storage option critical for meeting peak energy demands and realising green hydrogen ambitions. During early de-risking of specific sites, it is important to develop an accurate geological model to test whether the reservoir has the desired containment, volume and hydrogen deliverability. However, where seismic reflection lines and well data are limited and/or the storage system is structurally complex, the resulting geological models may be non-unique. Therefore, injection and withdrawal simulations using different structural end members is critical to constrain how a hydrogen plume may flow within (and out of) the container and interact with existing reservoir fluids. Here we present workflows for modelling a multi-year injection and withdrawal cycle of hydrogen into a depleted gas field. We use data from the Tariki Sandstone Member of the Ahuroa field in the Taranaki Basin, currently used to store natural gas in ANZ. This reservoir is located 2 km deep at the crest of an anticline above a major thrust fault, with marine mudstones forming the top seal and low-permeability fault rock the lateral seal. With only mixed quality 2D seismic reflection lines and a tight well cluster, the precise geometry of the thrust fault and its relations to smaller secondary faults is poorly constrained. To capture this uncertainty in our simulations, we have developed two 3D geological models of the Ahuroa field in Leapfrog Energy software. We use these geological models to conduct dynamic simulation of hydrogen injection and withdrawal using the massively-parallel simulator PFLOTRAN-OGS. We develop simulations that allow us to, over a 10-year cycle, test for closure or spill into adjacent fields, and predict the amount of mixing with remnant natural gas and formation water. During the simulations, we see major differences between the two geological models related to cushion injection and working H2 volumes, rates of water production and impurities due to natural gas. Additionally, one model has high risks of unrecoverable H2 gas loss when over-pressurised. Finally, we reimport the results back into Leapfrog for visualisation of the behaviour of the two hydrogen plumes over time.

Underground hydrogen storage (UHS) in depleted reservoirs presents a promising solution for managing seasonal variations in renewable energy during the global energy transition. However, the impact of reservoir heterogeneity, particularly permeability anisotropy and interlayer characteristics, on hydrogen recovery efficiency remains insufficiently understood. To bridge this knowledge gap and improve storage site selection accuracy, we developed a systematic box model to evaluate the effects of reservoir heterogeneity and validated our findings using New Zealand's Ahuroa gas storage field.Our investigation revealed that permeability anisotropy affects hydrogen recovery efficiency, with variations depending on well patterns. For well patterns with vertical wells only, both lateral (kx/ky) and horizontal-to-vertical (kh/kv) permeability anisotropy enhanced hydrogen recovery efficiency. For combined vertical and horizontal well patterns, the effect varied by anisotropy type. Lateral (kx/ky) anisotropy enhanced efficiency when horizontal wells aligned with the maximum permeability direction. In contrast, when horizontal wells aligned with the minimum permeability direction, kh/kv anisotropy exhibited an optimal ratio, beyond which efficiency began to decline. Analysis of interlayer effects revealed that reducing permeability from 1 mD to 10-3 mD led to an enhancement in hydrogen recovery efficiency, increasing from 61% to 75%. Additionally, our investigation demonstrated that the presence of interlayer pinch-outs and discontinuities along vertical hydrogen migration pathways reduced hydrogen recovery efficiency. A realistic geological model corroborated the box model findings: hydrogen recovery efficiency improved from 66.6% to 77.9% as the kx/ky ratio increased from 1 to 10, and from 66.6% to 76.2% when the kh/kv ratio increased similarly. Furthermore, inaccurate estimation of interlayer permeability could result in an 11.3% deviation in hydrogen recovery predictions.These results underscore the importance of accurately characterizing reservoir heterogeneity, including permeability anisotropy and interlayer properties, to ensure reliable hydrogen recovery predictions and improve site selection for UHS.

Hydrogen is projected to account for at least 10% of the global energy system in 20 years and is a critical component of the future zero-emissions energy system. Underground storage of green hydrogen in Aotearoa New Zealand (ANZ) will take advantage of intermittent surplus of renewable electricity at low cost, balance seasonal fluctuations in energy supply and demand, and provide a strategic reserve of energy. This poster is part of a larger research programme primarily focused on investigating the potential for underground hydrogen storage (UHS) in Taranaki, ANZ. Here, we explore the potential for UHS in porous rock formations of depleted gas reservoirs with particular focus on the role of seal integrity for storage. In this project the overarching goal is to improve understanding of whether mudstone seal strata have the potential to prevent leakage of hydrogen from Taranaki reservoirs. The primary focus is to characterise the geometries of fault and fracture systems in seal strata, their impact on its bulk permeability and to identify the pressure conditions required to promote the loss of seal integrity. In this poster we use Formation Micro Imagery (FMI) together with stratigraphic and fault/fracture mapping of core from petroleum wells to identify fracture densities, orientations and properties in both seal and reservoir rocks. Interpretations of seismic reflection lines in Taranaki and analogous outcrop observations are used to understand the geometries and permeability properties of fault zones. Preliminary results indicate that fractures are present in both reservoir and seal rocks. The densities of fractures increase with proximity to regional fold hinges and faults, and with increasing carbonate content. Questions remain about under what conditions fractures are open and capable of transmitting hydrogen. The poster outlines preliminary results, proposed research pathways and invites discussion.

The New Zealand Community Fault Model (NZ CFM) is a publicly available representation of New Zealand fault zones that have the potential to produce damaging earthquakes. Compiled through collaborative engagement between New Zealand earthquake-science experts, this first edition (version 1.0) of the NZ CFM builds upon previous compilations of earthquake-source active fault models with the addition of new and modified information. Developed primarily to support an update of the New Zealand National Seismic Hazard Model, the NZ CFM comprises two principal components. The first dataset is a two-dimensional map representation of the surface traces of 880 generalised fault zones. Each fault zone is assigned specific geometric and kinematic attributes, including uncertainties, supplemented with a subjective quality ranking focused primarily on the confidence in assigned slip rates. The second component is a three-dimensional representation of the fault zones as triangulated mesh surfaces that are projected down-dip from the two-dimensional mapped traces to a geophysically-defined maximum fault rupture depth. This article summarises the compilation and parameterisation of the NZ CFM, along with background on its relation to predecessor datasets, and forward applications to probabilistic seismic hazard assessment and physics-based earthquake models currently being developed for Aotearoa New Zealand.

We present a mid-Late Cretaceous to present day tectonic reconstruction model for Aotearoa-New Zealand. Our GPlates model comprises 50 rigid crustal blocks grouped into regions with common deformation histories set within a well-defined Australia-Pacific-Antarctica plate circuit tied to a published global paleomagnetic absolute reference frame. Within the model, four distinct periods of deformation are recognised from both near- and far-field observations. A key model assumption is the continuity of basement terranes between North and South Zealandia prior to Middle Eocene rifting and Late Oligocene initiation of transform motion. To complement the rigid crustal block model, continuously closing polygons show a ∼25% decrease in plate boundary area since the Middle Eocene that has been compensated for by corresponding increases in crustal thickness. A kinematic fault-propagation fold model demonstrates the plausibility of post-Late Oligocene asymmetric oroclinal folding on both sides of the evolving transform boundary. ‘Missing’ areas of map section common to previous tectonic reconstructions can be reconciled through contraction, elongation and vertical-axis rotation of continental crust within the deformation zone ahead of northward propagation of the Alpine Fault. This tectonic reconstruction model provides an open, accessible, and testable foundation for current and future paleogeographic and tectonic studies across Zealandia.

Evaluating fault segmentation is important for our understanding of seismic hazard assessment and fault growth. However, it is still unclear what controls if reverse fault earthquakes will rupture across segment boundaries. Here, we combine fault mapping and trench data from the low slip rate (0.04-0.15 mm/yr) multi-segment Nevis-Cardrona Fault (NCF) in the South Island of Aotearoa New Zealand to assess if it has ruptured in single or multi-segment earthquakes during the late Quaternary. Two new trenches on its Nevis segment provide stratigraphic evidence for two surface rupturing earthquakes, which through Optically Stimulated Luminscence dating and OxCal modelling, are constrained to have occurred at 28.9 +12.9 -9.1 ka and 12.8 ± 4.9 ka. The most recent timing is only weakly correlated to surface rupture timings from two trenches along the NCF's NW Cardrona segment. Furthermore, the 2 ± 1 m Nevis segment single event displacements we estimate would be unusually low for a ~85 km long NCF multi-segment rupture. We therefore surmise that late Quaternary NCF surface rupturing earthquakes did not rupture through ~30-50° bends that link these segments. Our trench data and fault mapping also indicate lower slip rates on the Nevis segment than previous studies (0.04-0.1 mm/yr vs 0.4 mm/yr).

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

Subduction zones have the greatest potential to generate large earthquakes and tsunamis. However, when undertaking Probabilistic Tsunami Hazard Assessments (PTHAs), subduction zones are a significant source of epistemic uncertainty. Therefore, understanding how the spatial distribution of elastic strain accumulation on the subduction interface influences the tsunami hazard is important for providing comprehensive hazard assessments, as well as quantifying uncertainty. This is especially important if the spatial locking distribution is undefined, and if it changes through time. Physics-based earthquake simulators allow different interpretations of the subduction interface locking distribution to be modelled, and how this influences the long-term seismicity, and the tsunami hazard, can be explored. Using three physics-based synthetic earthquake catalogues, generated by the earthquake simulator RSQSim, we analysed the tsunami hazard in Aotearoa/New Zealand. Three alternative representations of the subduction interface locking distribution along the Hikurangi Subduction Margin and the Tonga-Kermadec Subduction Zone were specified in the simulator to generate the catalogues. We modelled the tsunamis generated by $M_W\, \gt $8.0 earthquakes from each of the catalogues and undertook PTHAs. These assessments show that patches of high slip-deficit, both along strike and dip of the subduction interface, increase the tsunami hazard at the coast. Locking along the shallowest segments of the subduction interface also significantly increases the tsunami hazard. Our results show that careful consideration of the locking distribution in physical models is necessary before using them for PTHAs. They also show that by analysing multiple physical models of subduction zones, uncertainty in hazard assessments caused by the unresolved interface properties can also begin to be quantified.