Paleotsunamis on the Southern Hikurangi Subduction Zone, New Zealand, Show Regular Recurrence of Large Subduction Earthquakes Open Access

At transitional zones between subduction and strike‐slip plate boundaries, both fault systems can generate large earthquakes that produce similar geological signatures of vertical coastal deformation and tsunami (Clark et al., 2017). Disentangling the two earthquake sources is important for enabling the examination of fault interactions and earthquake triggering between the two fault systems (Goldfinger et al., 2008). Such interactions are particularly important at populated plate boundary transition zones such as the northern San Andreas fault–Cascadia subduction zone and the Marlborough fault system (MFS)–Hikurangi subduction zone (HSZ) where clusters of large earthquakes could extend the duration of an earthquake and tsunami sequence, impacting how society prepares and recovers. Although multiple historic examples have shown that large subduction earthquakes can trigger earthquakes on nearby upper‐plate faults (Melnick et al., 2012; Gomberg and Sherrod, 2014), these examples neither span transitional plate boundaries nor has there been documentation of large upper‐plate fault earthquakes triggering subduction earthquakes.

At the southern HSZ, New Zealand, the possibility of a large crustal earthquake triggering a subduction earthquake came to prominence following the 2016 $Mw$ 7.8 Kaikōura earthquake. Rupture of multiple faults in the MFS and possibly the subduction interface was followed by afterslip and slow slip on the subduction interface (Wallace et al., 2018; Wang et al., 2018). Increased coulomb stress on the locked portion of the interface beneath the southern North Island from this sequence of events may have increased the likelihood of a subduction earthquake (Kaneko et al., 2018). However, evaluating the short‐term likelihood of a subduction earthquake is challenging without the long‐term geological context of the frequency of subduction earthquakes and their relationship to the timing of large crustal earthquakes. Previous work at Mataora‐Wairau (M‐W) Lagoon identified geological evidence of two subduction earthquakes at $∼500$ and $∼850yr$ B.P. (Fig. 1b, Clark et al., 2015); and, although this data suggests a relatively short interval between earthquakes (405–280 yr), it is insufficient to ascertain the pattern of large earthquake occurrences on the subduction interface nor is the record of the two events long enough to examine a relationship with large earthquakes on upper‐plate faults that mostly have recurrence intervals (RIs) longer than 500 yr. Here, we present evidence for two paleotsunamis from Lake Grassmere that predate the earthquakes documented at M‐W Lagoon, extending the record of subduction zone ruptures in the region to 2000 yr. For the first time, we are able to calculate the RI of large subduction earthquakes on the southern HSZ and evaluate fault system interactions at this transitional plate boundary.

Lake Grassmere is situated in the northeastern South Island, within a transitional zone between the MFS and the HSZ (Fig. 1b). The subduction interface lies at a depth of 25 km below the site, and active faults lie in close proximity to north and south of the site. A number of historic earthquakes have caused vertical deformation at the study site. In the 2016 Kaikōura earthquake, rupture of Needles fault, located 12 km southeast of Lake Grassmere, caused 0.4 m of uplift at the study site, but there was no evidence of tsunami run‐up. The 1855 $Mw$ 8.2 Wairarapa earthquake, on the Wairarapa fault, its offshore fault extensions, and possibly the deep part of the subduction interface, caused widespread vertical deformation and tsunami run‐up around the southern North Island (Beavan and Darby, 2005). No tsunami run‐up was recorded at Lake Grassmere, but there is circumstantial historic evidence of uplift of Lake Grassmere in 1855 (Grapes and Holdgate, 2014).

Lake Grassmere is currently an artificially closed coastal lake used for salt production. The lake has a Holocene history of progressive isolation from the open sea, as it has evolved from an open bay (⁠$∼4000yr$ B.P.) to a barrier‐enclosed lagoon (⁠$∼2000yr$ B.P., Ota et al., 1995). Nontectonic sea‐level changes have been negligible (⁠$∼1m$ fall) in this region over the past 2000 yr (Clement et al., 2016).

To investigate the subsurface stratigraphy of the northern margins of Lake Grassmere, we obtained stratigraphic data from 30 gouge auger cores and 10 piston cores on transects normal and parallel to the modern lake margin (Fig. 1d). Unsplit cores were scanned using medical x‐ray computed tomography, then the cores were split and described using standard Troels‐Smith (1955) descriptors. A master core (core 6) and sections of other cores were sampled for particle size, foraminifera, and mollusk assemblage analyses (see supplemental material available to this article). Particle size distributions were obtained using a Beckman Coulter L13 320 laser granulometer. Core 6 was sampled at 2 cm intervals for foraminifera and diatoms, supported by samples from five additional cores. Diatoms were absent from all samples. Foraminifera samples were sieved to $63μm$⁠, and a count of 100 tests per sample was obtained where possible. Mollusks in seven cores were counted and identified, and then categorized into high tidal mudflat, intertidal estuary, or rocky intertidal preferred habitats (supplemental material). Well‐preserved, dominantly whole juvenile bivalves were selected for radiocarbon dating. Radiocarbon age modeling was conducted in OxCal 4.3.2 (Bronk Ramsey, 2009) using Marine20 curve (Heaton et al., 2020) with a local marine reservoir correction of $−162±24$ (from Clark et al., 2019, adjusted for the Marine 2020 calibration; code in supplemental material).

To evaluate potential tsunami sources, tsunami modeling was undertaken using Cornell Multi‐grid Coupled Tsunami model (Wang and Power, 2011). In the modeled scenarios, coseismic ground surface and seafloor displacements were computed with the elastic finite‐fault method of Okada (1985). For the Cloudy fault and Wairarapa fault, fault geometries were taken from the National Seismic Hazard Model (Stirling et al., 2012). For HSZ scenarios, 20 HSZ slip variation scenarios from Wang et al. (2016) were used to model the tsunami inundation in Lake Grasmere. A slip‐rate deficit distribution (from Wallace et al., 2012) was applied as a weighting mechanism to redistribute the randomly generated slip distributions to reflect the coupling status of the HSZ interface. The present‐day topography of the Lake Grassmere area was derived from high‐resolution photogrammetry models, and adjustments were made to the barrier height to account for paleotopography at the time of tsunami occurrence (details in the supplemental material).

We combined the Lake Grassmere tsunami record and M‐W Lagoon earthquake record (Clark et al., 2015) to calculate the mean RI and coefficient of variation (CoV) for earthquakes on the southern HSZ. We determined the 95% confidence limits on the RI and the mean CoV for the RI using a bootstrapping approach that produced 10,000 earthquake series by drawing from the earthquake age probability distributions with elimination of age reversals. Conditional probabilities for rupture of the southern HSZ were estimated using the method of Rhoades and Van Dissen (2003).

Five sedimentary units are continuous across the study site (Figs. 2 and 3). Units 1, 3, and 5 are dominantly fine‐grained silt, and units 2 and 4 are coarse shell hash within sand (Fig. 3). All units contain a near‐monospecific foraminiferal assemblage of Ammonia aoteana (80%–100%), which has a living range from low tidal to subtidal. Haynesina despressula occurs at $∼6%$ in unit 3 and above, possibly indicating a slightly higher salinity than unit 1 (Fig. S3). Units 1, 3, and 5 are interpreted as subtidal lagoon sediments based on their foraminifera assemblages.

The physical and biological characteristics of units 2 and 4 are highly anomalous within the quiescent, subtidal paleoenvironment of units 1, 3, and 5. We interpret these units as tsunami deposits due to their anomalous mollusk assemblages, high‐energy mode of deposition inferred from coarser grain size, and broad spatial extent. Units 2 and 4 are dominated by a densely packed hash of Austrovenus stutchburyi and Nucula hartvigiana. A transport process is necessary to account for the juxtaposition of this predominantly intertidal bivalve assemblage within the subtidal environment of the depositional site (Figs. 2 and 3). The shell hashes display a range in bivalve maturity from juveniles to mature adults, and we interpret the abundance of articulated bivalves to signify mass transport of live individuals during tsunami inundation, like Kitamura et al (2018). Rapid burial and a lack of reworking are indicated by the excellent preservation of all species, showing no signs of abrasion, encrustation, or rounding. The absence of mollusks in bounding units 1, 3, and 5 and lack of burrowing observed across the contacts further suggests that the shell hashes are transported deposits, rather than in situ accumulations of shells. Species diversity increases with distance inland, and we infer that taxa of rocky intertidal and high‐tidal mudflat habitats at the inner edge of the lagoon were mixed into the deposit (Fig. 3a–c). A striking feature of units 2 and 4 is the lateral continuity over 1.7 km inland, most likely beyond the reach of storm surges and incompatible with localized shell accumulations commonly seen at lagoon margins and in tidal channel lags elsewhere. Units 2 and 4 also display more conventional characteristics of paleotsunami deposits such as a coarser grain size and very sharp, erosional lower contacts, further indicating high‐energy deposition (Figs. 2 and 3).

Paleotsunamis are typically dated with bracketing ages above and below the tsunami deposit (Kelsey and Witter, 2020), but the bounding units of tsunamis 1 and 2 contained no reliable terrestrial‐ or marine‐dateable material at Lake Grassmere. Instead, we preferentially dated articulated, juvenile bivalves from within the tsunami deposits. Our age model is based on an assumption that the youngest ages in the deposit represent a death assemblage caused by the tsunami, and the tsunami age must overlap with and slightly postdate the youngest age cluster. We utilize “Tau” boundaries within the sequence priors of the age model to favor the youngest ages. We obtain age ranges of 2145–1837 yr B.P. for tsunami 1 and 1505–1283 yr B.P. for tsunami 2 (Fig. 3d).

There is circumstantial evidence that tsunami 2 was accompanied by coseismic uplift at Lake Grassmere. The finely laminated silt and clay of unit 1 imply a quiescent environment sheltered from tidal currents and wind waves. The nonlaminated, more variable, and coarser sediment of unit 3 indicates permanent change in depositional environment following the earthquake that generated tsunami 1, perhaps with increased tidal and/or wind‐wave mixing, compatible with uplift of the lagoon and shallowing of the water depth.

We consider possible distant and local tsunami sources that could inundate Lake Grassmere. Transoceanic tsunamis are discounted because large South America‐sourced tsunamis (the most significant distant source for New Zealand) in 1960 and 1868 did not significantly impact the Marlborough coast (Downes et al., 2017). Submarine landslides in Cook Strait are a possible tsunami source, but modeling suggests a low hazard from these sources for the southern Marlborough coastline (Lane et al., 2016). Most of the active faults within Cook Strait are dominantly strike slip, and marine seismic profiles show negligible vertical offset on the submarine scarps (e.g., Wairau fault and Awatere fault); therefore, they are unlikely to produce sufficient seafloor displacement to generate large tsunamis, but there are some exceptions to this that we have tested. We consider rupture of the southern Hikurangi subduction interface to be the most likely source of a large tsunami at Lake Grassmere.

To test the hypothesis that paleotsunamis observed at Lake Grassmere could only be generated by subduction earthquakes, we modeled tsunamis from three different fault sources and compared modeled inundation with the extent of the paleotsunami sediments at our study site (Fig. 4). Scenario A models the Cloudy fault (Fig. 4a), because it is a dip‐slip fault with optimal orientation for propagation of a tsunami toward Lake Grassmere. Scenario B models a tsunami generated by the Wairarapa fault and its offshore extension (the Wharekahau thrust, Fig. 4b) in a replication of the A.D. 1855 $Mw$ 8.2 earthquake that generated a tsunami with run‐ups of 10 m in the southern North Island. Relatively small‐wave heights and a lack of inundation at Lake Grassmere suggest that the Cloudy and Wairarapa faults are not plausible matches for tsunamis 1 and 2 (Fig. 4a,b).

The third tsunami source is the subduction interface on which tsunami generation is strongly influenced by the slip distribution, particularly how much and how far slip propagates toward the trench. We modeled 20 different slip distributions for an $Mw$ 8.9 earthquake on the HSZ, and the majority of these produce inundation at paleo‐Lake Grassmere (supplemental material). The preferred model shown in Figure 4c has a slip distribution consistent with geological observations of net Holocene uplift at Lake Grassmere and subsidence at M‐W Lagoons (Ota et al., 1995; Clark et al., 2019). Consistency with long‐term geological observations of vertical deformation may indicate preferred subduction interface rupture patch, but each earthquake probably has a different slip distribution, and interseismic recovery of vertical deformation could mean that the long‐term signal does not reflect single‐event displacements. Our preferred subduction interface model shows wave heights of up to 10 m at the coast of Lake Grassmere and flow depths of $≤5m$ within the lake basin. Of our 20 different models of an $Mw$ 8.9, all show some inundation of Lake Grassmere, with 10 models showing $>2.5m$ of flow depth and one showing $>5m$ flow depth. When the preferred tsunami model (Fig. 4c) is scaled down to an $Mw$ 8.5 earthquake, it does not show significant inundation at Lake Grassmere. The results of the tsunami modeling support the hypothesis that the most probable source of paleotsunamis at Lake Grassmere is rupture of the HSZ; moreover, these past earthquakes were probably $Mw>8.5$⁠.

Combining evidence from Lake Grassmere and M‐W Lagoon generates a record of four subduction earthquakes in the past 2000 yr (Fig. 5a). Ideally, evidence for each past earthquake would be replicated at many sites along the margin. However, local‐scale correlation of subsidence and tsunami records on the southern HSZ is not evident due to rapid changes in preservation potential of the sites through time. Uplift of Lake Grassmere in the hanging wall of the Needles fault has progressively isolated the site from the sea, raised the barrier, and lessened its sensitivity and accommodation space for recording tsunamis, so the site does not record the youngest two earthquakes seen at M‐W Lagoon. Conversely, prior to $∼1500yr$ B.P., the M‐W Lagoon was an intertidal embayment (Hayward et al., 2010), and less sensitive to recording subsidence and tsunami. Evidence of the $∼2145–1837yr$ B.P. tsunami is possibly recorded at M‐W Lagoon, but it is poorly dated (King et al., 2017). The difficulty in getting long subduction earthquake records at a single site is a characteristic of transitional plate boundaries where permanent vertical deformation on crustal faults alters the site sensitivity to sea‐level change and tsunami inundation through time, in addition to well‐documented changes in site sensitivity due to nontectonic sea‐level change (Dura et al., 2016). Despite this, the last three subduction earthquakes on the southern HSZ show far‐field correlations with subsidence at Ahuriri Lagoon on the central HSZ and marine terrace uplift on the Wairarapa coast (Fig. 1, Clark et al., 2019) that support a subduction earthquake signal.

Integration of the Lake Grassmere and M‐W Lagoon subduction earthquake records enables the first calculation of large earthquake recurrence statistics for the locked patch of the southern HSZ. We find an RI of 500 yr (335–655 yr, 95% confidence interval) and a CoV of 0.27 (0.0–0.47, 95% confidence interval, Fig. 5c,d). Physical evidence of tsunamis accompanying at least three of the last four subduction earthquakes confirms a high tsunami hazard for the central New Zealand region and points to a high likelihood that ruptures extend well into Cook Strait, notably increasing tsunami hazard for northeastern South Island, including Christchurch (Power et al., 2016). Considering the timing of the most recent southern Hikurangi subduction earthquake (⁠$∼520–470yr$ B.P.) alongside our calculated RI and low CoV, we infer that the southern HSZ is late in its seismic cycle. We calculate a conditional probability of 26% for a southern HSZ earthquake in the next 50 yr. The proximity of the subduction interface to New Zealand’s capital city, Wellington, and the high likelihood of a tsunami across the Cook Strait region means that a subduction earthquake will probably have significant environmental and societal impacts.

The southern HSZ is one of the very few transitional plate boundaries globally for which robust paleoseismic records exist for surface‐rupturing earthquakes on the majority of upper‐plate faults over the past 2000 yr (Fig. 5a). When combined with this new subduction earthquake record, it facilitates a rare opportunity to examine the interaction between the subduction interface and the upper‐plate faults at a transition zone. Temporal association between subduction interface and strike‐slip rupture has been inferred at the transition between the Cascadia subduction zone and the San Andreas fault (Goldfinger et al., 2008), but we have not previously had the length of record required to examine this relationship on the southern HSZ.

Comparison of earthquake ages shows that there is no consistent relationship between upper‐plate fault earthquakes and subduction zone ruptures. All the past subduction earthquakes (gray shading, Fig. 5a) temporally coincide with at least one upper‐plate fault earthquake. However, this result is not surprising, because the low age resolution on many paleoearthquakes means that there are almost no time periods without an upper‐plate fault earthquake, particularly in the record $>1000yr$ B.P. One notable cluster of earthquakes is at 880–800 yr B.P. when a subduction earthquake coincides with large earthquakes on the Wellington, Wairarapa and Awatere faults, and the earthquake at 2145–1837 yr B.P. overlaps with four upper‐plate fault earthquakes, but all have large uncertainties. In contrast, the subduction earthquakes of 1509–1253 and 520–470 yr B.P. overlap with only one upper‐plate fault paleoearthquake each.

The lack of a clear relationship between upper‐plate fault earthquake and subduction earthquake timing over the past 2000 yr could be attributed to poor age resolution that inhibits detection of a pattern, or a true reflection that the two fault systems operate largely independently of another. Although recent upper‐plate fault earthquakes (2016 Kaikoura and 1855 Wairarapa) may have had a component of interface slip, the slip was at the fringes of the currently locked patch of the interface, so they are not comparable with the large paleoearthquakes that we infer rupture the locked patch of the southern HSZ (Fig. 4d).

The low CoV of the southern Hikurangi interface is consistent with quasiperiodic earthquake behavior driven by steady plate motion moment accumulation. Upper‐plate fault earthquakes on the southern HSZ appear to display a more clustered occurrence through time, but there is no consistent correlation with subduction earthquakes. This finding implies that the recent 2016 Kaikōura earthquake probably had little effect on the seismic cycle of the subduction interface because the subduction interface appears to maintain quasiperiodic earthquake behavior irrespective of upper‐plate fault activity.

The supplemental material for this article includes full core descriptions and biostratigraphy, radiocarbon dates and age models, and additional information on tsunami modeling and calculation of conditional probabilities. The sediment cores are archived at GNS Science. The digital surface model of Lake Grassmere is stored in the GNS Science digital elevation collection and available upon request. All other data used in this article came from published sources listed in the References.

The authors declare no competing interests.

This research was supported by the It’s Our Fault Programme (funded by Earthquake Commission, Wellington City Council and Greater Wellington Regional Council), and the GNS Science Understanding Zealandia and Hazard and Risk Management Strategic Science Investment Fund Programmes funded by the Ministry of Business, Innovation and Employment. The authors thank Peter Meihana (Ngati Kuia and Rangitāne o Wairau), landowners Mr & Mrs Spence and Peter Yealands, and Andrew Howell, Alan Beu, and staff of the Rafter Radiocarbon Lab for their scientific advice. Two anonymous reviewers improved this article.