# 14 Years of Slip on the Hikurangi Subduction Zone

Sylvain Michel, Romain Jolivet, Emilie Klein, Louise Maubant

► To cite this version:

Sylvain Michel, Romain Jolivet, Emilie Klein, Louise Maubant. 14 Years of Slip on the Hikurangi Subduction Zone. Journal of Geophysical Research: Solid Earth, 2025, 130 (7), pp.1-26. 10.1029/2024jb030865. hal-05472282

HAL Id: hal-05472282

https://asnr.hal.science/hal-05472282v1

Submitted on 22 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2024JB030865

**Key Points:**

* We model 14 years of slip along the Hikurangi subduction megathrust, devoid of temporal smoothing, between 2009 and 2023 based on GNSS data

* We extract a catalog of 27 slow slip events (SSEs) and build maps of SSEs cumulative slip, inter-SSE coupling, and long-term coupling

* We observe a change in the behavior of SSEs after the 2016 Kaikoura earthquake and discuss the interplay between coupling and SSE

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

S. Michel,

sylvain_michel@live.fr

**Citation:**

Michel, S., Jolivet, R., Klein, E., & Maubant, L. (2025). 14 years of slip on the Hikurangi subduction zone. *Journal of Geophysical Research: Solid Earth*, 130, e2024JB030865. https://doi.org/10.1029/2024JB030865

Received 27 NOV 2024

Accepted 25 JUN 2025

**Author Contributions:**

**Conceptualization:** Sylvain Michel, Romain Jolivet

**Data curation:** Emilie Klein

**Formal analysis:** Sylvain Michel, Romain Jolivet, Emilie Klein, Louise Maubant

**Funding acquisition:** Romain Jolivet

**Investigation:** Sylvain Michel, Emilie Klein

**Methodology:** Sylvain Michel, Romain Jolivet, Emilie Klein

**Supervision:** Romain Jolivet

**Validation:** Sylvain Michel, Romain Jolivet, Emilie Klein, Louise Maubant

**Visualization:** Sylvain Michel, Emilie Klein

© 2025. The Author(s). This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

## 14 Years of Slip on the Hikurangi Subduction Zone

Sylvain Michel¹,²,³ ID, Romain Jolivet²,⁴ ID, Emilie Klein², and Louise Maubant⁵,⁶ ID

¹Institut de Radioprotection et de sûreté Nucléaire, Fontenay-aux-Roses, France, ²Département de Géosciences, Laboratoire de Géologie, École Normale Supérieure, PSL Université, CNRS UMR, Paris, France, ³Now at Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Valbonne, France, ⁴Institut Universitaire de France, Paris, France, ⁵Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA, ⁶Now at Earthquake Research Institute, University of Tokyo, Tokyo, Japan

**Abstract** Slow slip events (SSEs) on the Hikurangi subduction megathrust in New Zealand have been studied extensively. However, we lack a comprehensive and consistent catalog over a long time period since 2010. Recent studies have focused on short time periods and rarely provide a time dependent image of slip during SSEs, let alone between SSEs. In addition, the use of different methods, parametrizations and data sets does not allow to compare SSEs at different periods of time. Our study aims to model the spatio-temporal history of slip of the subduction, devoid of temporal smoothing, over the entire 2009–2023 time period based on GNSS time series. We combine an Independent Component Analysis decomposition with an inversion scheme to derive a history of accelerations and slowdowns of aseismic fault slip on the megathrust. We then combine a SSE detection procedure with an inter-SSE linear trend correction to reference this slip time series to the inter-SSE loading rate and build a catalog of 27 SSEs. We observe a change in the behavior of SSEs after the 2016 Kaikoura earthquake with respect to the pre-Kaikoura period, highlighted by changes in the slip frequency content. We finally evaluate using a Bayesian approach the subduction inter-SSE coupling map, the SSEs cumulative slip between 2009 and 2023, to constrain the distribution of locked regions over 14 years and conclude that moment deficit accumulates mostly under Wellington peninsula at a rate of $4.9 \pm 0.5 \ 10^{18}$ N.m/yr.

**Plain Language Summary** Slow slip events (SSEs) are just like earthquakes, releasing tectonic stress, but slowly, without radiating destructive seismic waves. SSEs along New Zealand's Hikurangi subduction zone have been studied extensively but in separate studies with a wide variety of methods, making it hard to track evolutions in the behavior of these intriguing events over the last two decades. Our study aims to model the history of slip on the subduction zone over a long time period, from 2009 to 2023, using GPS data. We combine signal decomposition methods with inverse problem theory on GPS data to build a consistent catalog of 27 SSEs. Our findings indicate a change in SSE behavior following the 2016 Kaikoura earthquake with noticeable variations in slip patterns. We also map areas where the subduction zone is locked between periods of SSEs and also over the 14 years studied. We conclude that the energy build up below the Wellington peninsula since the 1430 earthquake could possibly generate a magnitude 8.1 if it occurred today.

## 1. Introduction

The Hikurangi subduction zone is characterized by the sliding of the Pacific oceanic plate, and more specifically the Hikurangi plateau, under the North Island of New Zealand and the northern part of South Island (Figure 1). The oblique convergence between the Pacific and Australian plates reaches about 45 mm/yr. In addition, internal crustal deformation within North Island induces an additional rotation of the Island with horizontal GNSS displacement rates reaching up to ~20 mm/yr. This rotation results in a motion of convergence roughly perpendicular to the subduction trench, with rates ranging from ~10 to ~60 mm/yr from south to north (Wallace et al., 2012, 2004). The direction of convergence on the subduction interface changes also with depth. Along the deeper portion of the fault (>20 km depth), relative motion is parallel to the Pacific-Australia plate convergence direction while along the shallow portion (<20 km depth), crustal rotation leads to trench-normal convergence.

This change in relative motion appears clearly when analyzing Slow Slip Events (SSEs). Those events, first detected on the Cascadia and Japan subduction zones (Dragert et al., 2001; Hirose et al., 1999), correspond to aseismic slip events on the subduction interface which can last days to years (Behr & Bürgmann, 2021). During these events, slip is not fast enough to radiate seismic waves as earthquakes do, although SSEs can be accompanied by microseismicity and/or tremors (Rogers & Dragert, 2003). Numerous SSEs have been documented in
Writing – original draft: Sylvain Michel

Writing – review & editing:

Romain Jolivet, Emilie Klein,
Louise Maubant

<figure><img src="image_3.png" /></figure>

**Figure 1.** Regional setting of the Hikurangi subduction zone. The coupling map (yellow to red shading) is estimated from our study. Green contours indicate the cumulative slip patterns of slow slip events (100 mm intervals) during the 2009–2023 period as also evaluated in this study. Black dots indicate GNSS stations (https://data.geonet.org.nz/gnss/). The blue diamonds in the south and north locate the stations LEVN and CNST, respectively, for which the time series are shown in Figure 2. Gray thin lines indicate faults (Litchfield et al., 2014). The dark blue line in the South-west shows the 2016 Kaikoura earthquake trace (Zinke et al., 2019). White stars represent the epicenters of $M \ge 6.5$ earthquakes above 80 km depth (ANSS catalog; U.S. Geological Survey, 2017). The black arrow indicates the direction of convergence between the Pacific and Australian plates (Wallace et al., 2012). The bottom left inset shows the moment deficit rate PDF of the subduction.

the Hikurangi subduction zone (e.g., Bartlow et al., 2014; Wallace & Beavan, 2010; Wallace, Beavan, et al., 2012; Wallace et al., 2016), where they can be separated in two distinct families: (a) the deep SSEs (>20 km depth) last several years and occur every 5–10 years, and (b) the shallow SSEs (<20 km depth) last days to weeks with a recurrence of the order of a few months to years (Wallace, 2020).

In New Zealand, while micro-seismicity and tremors have often been observed during SSEs, their presence is far from being systematic (Aden-Antoniów et al., 2024; Fasola et al., 2023; Jacobs et al., 2016; Nishikawa et al., 2021; Romanet & Ide, 2019). Additionally, multiple interactions between aseismic slip on the subduction interface and large seismic events ($M_w > 6$) have also been documented. In 2013, a $M_w$ 6.3 earthquake is thought to have modified the propagation of a large SSE (Wallace et al., 2014). This same SSE in 2013 might have promoted a sequence of crustal earthquakes, which included two $M_w$ 6.5 in the Cook Strait, in between the North and South Islands. The 2016 $M_w$ 7.8 Kaikoura earthquake has ruptured multiple crustal fault segments but also a portion of the subduction interface. In addition to post-seismic slip on the subduction interface, the 2016 earthquake triggered in addition shallow SSEs more than 200 km away (e.g., Hamling, 2020; Wallace et al., 2018).
<table><caption>**Table 1**

*General Information on the Time Series Analyzed*</caption><thead><tr><th>Data set name</th><th>Time period analyzed</th><th>Time period over which time series are detrended</th></tr></thead><tbody><tr><td>A</td><td>2009 to Kaikoura earthquake (~2016.8)</td><td>2009–2013</td></tr><tr><td>B</td><td>Kaikoura earthquake to 2023</td><td>2009 to Kaikoura earthquake (The data set after the Kaikoura earthquake is thus detrended based on data before this event)</td></tr></tbody></table>

SSEs have been extensively studied in the Hikurangi subduction zone. However, the last comprehensive and consistent catalog over a long time period was established in 2010 (e.g., Wallace & Beavan, 2010; a description of SSEs during the 2002–2010 period). Recent studies focused on short time periods (e.g., Bartlow et al., 2014; Wallace, Beavan, et al., 2012; Wallace et al., 2018; Yan et al., 2023). Few studies also provide the time dependent history of slip, which is, in general, impacted by temporal smoothing (e.g., Bartlow et al., 2014; Yan et al., 2023). In this study, we derive the spatio-temporal history of slip on the subduction interface based on GNSS data from 2009 to 2023, a period during which the GNSS network is stable and in the continuity of the study of Wallace & Beavan (2010). We avoid using temporal smoothing by using a signal decomposition technique called variational Bayesian Independent Component Analysis on GNSS time series (vbICA; Choudrey, 2002; Choudrey & Roberts, 2003) (see Section 3). vbICA has been successfully used in other subduction zones and tectonic regimes to reconstruct surface displacements of GNSS stations due to SSEs (e.g., Michel et al., 2019a, 2022). We use the recovered surface deformation to infer concomitantly the temporally unsmoothed spatio-temporal history of slip, a catalog of SSEs and evaluate the velocity field of the GNSS stations related to the fault loading in between SSEs (See Section 4).

In order to provide a full budget of slip along the megathrust, we estimate maps of inter-SSE coupling (i.e., the distribution of the locked regions on the fault between SSEs), SSEs final cumulative slip, and long-term coupling (i.e., the distribution of the locked regions over the 14 years of observation) using a Bayesian approach (Sections 5–7). Such approach enables us to carefully evaluate associated uncertainties. We finally list the potential limitations and biases of our study (Section 7) before discussing the implications in terms of seismogenic potential of the Hikurangi subduction zone. Our study hence provides one of the first continuous and methodologically consistent view of the dynamics of slow tectonic deformation along the Hikurangi subduction zone over the 2009–2023 period.

## 2. Data

We process GNSS position time series using the GAMIT/GLOBK software relying on the classical double-difference MIT methodology (Herring et al., 2010a, 2010b). We include New Zealand data (RINEX were downloaded from https://data.geonet.org.nz/gnss/; GNS Science, 1995) and IGS data from stations located in Australia and Melanesia in order to establish a stable reference. We process in total 177 stations divided into 4 sub-networks of slightly more than 50 stations each and half a dozen stations in common between sub-networks. We follow the same methodology as described in detail in Klein et al., 2023. Position time series are mapped in ITRF2014 using PYACS (Altamimi et al., 2017).

We use the resulting daily east, north and up position time series between 2009 and 2023. Time series are split into two periods separated by the 2016 Kaikoura earthquake, respectively named data set or period A and B (Table 1; Figure 2). Stations with more than 20% of missing data within those two periods are removed. We also remove stations below latitude 42°S and stations close to the Kaikoura earthquake (<30 km from its rupture trace) to reduce the number of time series in which the signal is dominated by this event. A total of 131 and 130 stations are selected for the two periods, respectively (see list in Table S1 in Supporting Information S1).

We first correct the time series from instrumental and earthquake-related offsets, including two M6.5 earthquakes in 2013 and a M7.3 in 2021 (see Text S1 in Supporting Information S1 for details; Okuwaki et al., 2021; U.S. Geological Survey, 2017). We then evaluate the secular trends on each time series from data set A using a trajectory model (see Text S2 and S3 in Supporting Information S1). Trends for data set A are estimated over the period 2009–2013, before a known deep SSE which started at the beginning of 2013 and dominates the signal (Wallace, 2020). Since post-seismic motion severely affects time series after 2016, we cannot safely evaluate trends on data set B. We, therefore, detrend data set B with the secular trends estimated on data set A from 2009 to
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Illustration of the data sets and the detrending strategy (Section 2). (a) Raw east position time series of station LEVN. The data is separated in two data sets, A and B for the data before and after the Kaikoura earthquake, respectively. Each data set is detrended based on data from a specific period as indicated by the red and blue arrows. Note that data set B is detrended using data before the Kaikoura earthquake. (b) LEVN detrended east position time series for data set A (in red) and B (in blue). After detrending, data set A only considers the data before the Kaikoura earthquake (dark red). The time series after the earthquake are disregarded (pink). Instead, after detrending, data set B only considers the data after the Kaikoura earthquake (dark blue) and disregards the data before the earthquake (light blue).</figcaption></figure>

2016. This difference in detrending strategy is based on multiple tests realized from the results of the ICA decomposition (see Section 3), on how well the ICA seemed to visually manage to separate the different signals. In both cases, we do not interpret the trends as the long-term secular displacement rate of the GNSS stations and re-evaluate this long-term rate after detection of the SSEs (see Section 4.2). Those periods have only been chosen to facilitate the ICA decomposition. Examples of trended and detrended time series for two stations are shown in Figure 3.

## 3. Signal Decomposition

### 3.1. Method

To separate tectonic signals from that from other sources of deformation affecting GNSS time series (e.g., seasonal signals potentially due to hydrological loading), we apply a variational Bayesian Independent Component Analysis (vbICA; Choudrey, 2002; Choudrey & Roberts, 2003; Gualandi et al., 2016; Gualandi & Liu, 2021; Larochelle et al., 2018; Maubant et al., 2020; Michel et al., 2019b, 2022). This method separates time series into statistically independent components, each of them describing the spatial pattern and temporal evolution of a source of deformation. We perform the decomposition assuming the temporal evolution of each source is independent (in opposition to their spatial pattern). Resulting independent sources are additionally assumed to be static in space. Consequently, non-stationary sources, such as propagating SSEs, cannot be represented by a single component but by the combination of multiple sources (Larochelle et al., 2018; Michel et al., 2019b). The number of components to extract is chosen here based on the method described by Gualandi et al. (2016) (see Text S4 in Supporting Information S1).

Position time series are organized as a matrix $X_{obs}$ of size $M \times T$ where $M$ is the number of time series (three times the number of GNSS stations as we use East, North and Up components) and $T$ is the number of epochs. $X_{obs}$ is assumed to be a linear combination of signals produced by each source:

$$X_{obs} = A\Sigma + \varepsilon, \qquad (1)$$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Position time series of stations LEVN (panels a–f) and CNST (panels g–l) (blue diamonds in Figure 1) from data set A. Black dots and gray errorbars in panels to the left show the raw position time series and associated uncertainties. Rightmost panels show the detrended time series of both stations (see Figure 2 and Table 1 for details about the detrending). The red dots represent the reconstructed time series from the ICA decomposition taking only the components related to slow slip events (SSEs). The pink bars are their associated uncertainty. Light blue show the predicted time series based on the SSE kinematic inversions (Section 4.1). The results for all stations are shown in Figures S8 and S9 in Supporting Information S1.</figcaption></figure>

where *A* is called the mixing matrix (giving the weight of each source), Σ the sources matrix and ε the noise assumed here to be Gaussian noise centered on 0. Assuming a number of sources *L*, the size of *A* and Σ are *M* x *T* and *M* x *L*, respectively. Based on this structure, *A* describes the spatial pattern of each source while Σ describes their temporal evolution. In this study, we normalize each source in Σ between 0 and 1 so that their amplitude in *A* may be interpreted as the spatial pattern of maximum deformation.
For details on the technique, refer to Choudrey and Roberts (2003) and Gualandi et al. (2016). Additional details
on how we apply vbICA on our data set are available in Text S4 in Supporting Information S1.

**3.2. Signal Extraction**

We apply the same strategy for both data sets A (before Kaikoura) and B (after Kaikoura). The method of selection of the number of sources to extract from Gualandi et al. (2016) is quite conservative and is controlled by the signals of largest amplitude. Small signals of interest might be rejected by the method due to the dominating signals. To minimize this effect, we perform the ICA twice. First, we extract the sources of largest amplitude, and second, we decompose the time series corrected from the sources of large amplitude to retrieve smaller signals. We define here sources of large amplitude those with a maximum displacement above 25 mm. Following this approach, 8 and 15 components are extracted for data set A, while 9 and 13 components are extracted for data set B. The spatial pattern and temporal evolution of a selection of different sources are shown in Figures 4 and 5 for both data sets, respectively. All components are shown in Figures S1 and S2 in Supporting Information S1.

We then interpret each of those as either being of tectonic origin or not. Any seasonal signal or signal that looks
temporally like noise are put aside. The other sources, in total 10 and 11 components for data sets A and B,
respectively, are considered as of tectonic origin and kept for slip reconstruction (Table S3 in Supporting In-
formation S1; Section 4). Almost all of these identified components show temporal behavior that can be inter-
preted as potentially originating from tectonic sources and surface deformation spatial patterns consistent with
what is expected to result from slip on the subduction interface (i.e., deep and shallow SSEs will imprint an
oblique and perpendicular pattern relative to the trench, respectively). For data set B, one or two components seem
to be related, if not completely at least in part, to post-seismic deformation on crustal faults due to the 2016
Kaikoura earthquake. Those will nevertheless be kept for the inversion scheme. We acknowledge that within the
selected components there might still be some deformation due to volcanoes. We note that contrarily to similar
studies (e.g., Michel et al., 2019b, 2022), we do not extract any clear common mode, which have probably been
reduced thanks to the continental stabilization made with PYACS (Section 2, Klein et al., 2023).

**4. SSEs Kinematic Inversions and Detection**

With the selected components, we aim to detect and constrain the occurrence of SSEs along the subduction zone, reconstruct the slip evolution in space and time and provide a new catalog of SSEs for the region, consistent over the 2009–2023 period. Since we work on detrended time series, any slip distribution we would reconstruct from these components will hence be relative to the inferred trend for each station, which might not be representative of the interseismic period. We, therefore, start by reconstructing the history of slip on the megathrust from the selected independent components. We then correct this slip history by iteratively repeating a two-step approach: first, we detect any deviation from the trend (i.e., an SSE), and second, we remove a linear term estimated based on the period in between each detection. Our final slip history is, therefore, relative to the most stable slip rate for each point on the fault over the whole observation period.

**4.1. Slip History on the Hikurangi Subduction Zone**

We first invert for the slip history on the Hikurangi subduction interface using the components obtained in Section 3. We thus invert the spatial pattern of each selected component (i.e., the mixing matrix *A* in Section 2) for the slip distribution on the fault. This slip distribution is then combined with the temporal evolution enclosed in the corresponding temporal component. The resulting slip evolutions are then summed over all components to produce the final time-dependent slip distribution along the fault. Once again, since the position time series have been detrended, this slip history is relative to an average slip rate which is different in both data sets (see Table 1 and Section 2).

The inversion strategy is similar for both data sets and each component. We consider the section of the sub-
duction which is facing North Island and the northern part of South Island (roughly 715 km long along the fault
trench). The geometry of the subduction is based on Williams et al. (2013). We discretize the subduction
interface into 343 triangular subfaults. The area of the subfaults increases as they are located further away from
the eastern coast where the resolution is expected to be the highest. The average subfault area is ~585 km² and
the deepest subfault node is at ~128 km depth. For components 1 and 2 of ICA 1 of data set B, we add three
rectangular crustal faults, reaching from the surface to 25 km depth, near the northern tip of South Island to
<figure><img src="image_6.png" /><figcaption>Figure 4.</figcaption></figure>
roughly model the post-seismic processes due to the Kaikoura earthquake. The geometry of those faults is coarse as very few stations used in this study are close to the faults and able to constrain the slip inversion. The Green's functions relating slip on the megathrust to surface displacement included in the matrix G are computed assuming a layered elastic half-space (Text S5 and Table S4 in Supporting Information S1) and consider both strike and dip slip.

We consider slip in each component can be either positive or negative since each spatial pattern is the projection of the observed displacements on a given component. Assuming a Gaussian prior for all model parameters (i.e., strike slip and dip slip) and considering a Gaussian likelihood, the general solution to the linear inverse problem aiming at finding the slip distribution for each component is a multimodal Gaussian PDF describing the distribution of all possible models (e.g., Tarantola, 2005). The mean of this multimodal Gaussian distribution writes as

$$ \tilde{m} = m_0 + C_m G^T (GC_m G^T + C_d)^{-1} (d - Gm_0), \quad (2) $$

where d is the geodetic displacement vector, m the modeled slip vector on the fault, m₀ the prior slip model (here set to the null vector), and Cₘ and C₉ the model and data a priori covariance matrices, respectively. Uncertainties can be described by the a posteriori covariance Cpost such as

$$ C_m^{\text{post}} = C_m - C_m G^T (GC_m G^T + C_d)^{-1} GC_m. \quad (3) $$

We consider the prior covariance between slip parameters follows an exponential decay based on the distance between patches (Radiguet et al., 2011) The covariance of the model Cₘ is defined as

$$ C_{m \ i,j} = \left( \sigma_0 \frac{\lambda}{\lambda_0} \right)^2 \exp \left( -\frac{d_{i,j}}{\lambda} \right), \quad (4) $$

where σ₀ is the model parameters a priori uncertainties, λ₀ is a scaling factor set to 24 km (the characteristic subfault size), λ is a correlation length set to 45 km (almost twice λ₀), and dᵢⱼ is the distance between subfaults i and j. We explore the value of σ₀ which defines L-curve plots between the data misfit and the Euclidean norm of the model, ||m||₂, and select a model located at the bend of the L-curve corresponding to a good balance between misfit and the smoothness of the model (Table S5 in Supporting Information S1; Figures S3 and S4 in Supporting Information S1). The slip models of all components together with their misfit are shown in Figures S5 and S6 in Supporting Information S1. As indicators to evaluate the robustness of the inversions, the resolution, restitution and sensitivity maps are shown in Figure S7 in Supporting Information S1. The surface displacements predicted by the final model are shown in Figure S8 in Supporting Information S1 (for Data set A) and Figure S9 in Supporting Information S1 (for Data set B) superimposed on the signal extracted from the ICA and the detrended initial data.

## 4.2. SSEs Detection and Inter-SSE Trend Correction

The time-dependent slip evolution obtained in Section 4.1 is relative to a linear slip trend which is dependent on the time period used to initially detrend the GNSS position time series and on the number of potential SSE within this period (Section 2; Figure 2). The linear trend as well as the SSEs detection threshold thus vary spatially. To minimize the impact of potentially heterogeneous detection threshold, we use an iterative process that combines a SSE detection procedure based on slip rate with an inter-SSE linear trend correction determined from the previous iteration of SSEs' detection (Figure 6).

**Figure 4.** A selection of components extracted by the ICA for data set A (pre-Kaikoura data set). Top panel of each component shows the corresponding temporal evolution while the bottom panel shows the spatial pattern. Horizontal and vertical displacement are indicated by the black arrows and colored dots, respectively. The temporal evolution of each component is normalized and the amplitude of the spatial pattern shows the maximum displacement possible. The length of the black arrow legend is set to the maximum horizontal displacement for each component and allows thus to compare the maximum amplitude between the components. Similarly, the extent of the colorbars are set to the maximum vertical absolute amplitude of each component. Components 1 to 4 are interpreted as related to slow slip events while components 5 and 6 are interpreted as related to non-tectonic sources. Full results of the ICA are shown in Figure S1 in Supporting Information S1.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Same as Figure 4 but for data set B (post-Kaikoura data set). Components 1 to 5 are assumed related to slow slip events and/or post-seismic slip from the Kaikoura earthquake while component 6 is assumed non-tectonic. The full results of the ICA are shown in Figure S2 in Supporting Information S1.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Illustration of the slow slip event (SSE) detection and inter-SSE trend correction procedure (Section 4.2) for a node of the fault located at longitude and latitude (178.40; -38.64). For each iteration, we follow the same series of steps. (a) The slip time series (black dot) at the fault node is temporally smoothed using a zero-phase digital filter with a window of a chosen temporal size *T*<sub>win</sub> here fixed at 30 days (blue curve). (b) The slip rate time series at the node (black curve) is computed based on the smoothed slip time series and a SSE detection is applied based on a slip rate detection threshold, δ̂<sub>thresh</sub>, fixed here at 75 mm/yr (horizontal green dashed line). The detection periods are indicated by the orange rectangles at the bottom of the panel. (c) The inter-SSE periods of the slip time series are isolated based on the SSE detections. We then fit a linear trend with Heaviside functions at the SSEs timing through the inter-SSE periods of slip (red curve). Finally the estimated linear trend is remove from the slip time series before proceeding to the next iteration. The black curve in Panel (a) of the 2nd iteration correspond to the slip time series at the fault node corrected from the slope evaluated at the end of the first iteration.</figcaption></figure>

SSEs detections are based on the slip rate history of the fault. We first apply a temporal filter on the slip history and calculate the time-dependent derivative of slip (Figure 6b). Temporal smoothing is performed using a zero-phase digital filter with a window of a chosen temporal size, $T_{win}$ (matlab filtfilt function, with $b$ = ones (1, $T_{win}$)/$T_{win}$ and $a$ = 1). The filter is a moving average window that runs over the data both forward and backward, smoothing out variations with periods shorter than approximately $2T_{win}$. We detect a SSE when the slip rate over contiguous fault nodes exceeds a fixed detection threshold, $\delta_{thresh}$. Multiple SSEs can hence coexist at a given time step. We then track those SSEs in time. Overlapping spatial patterns of SSEs at time *t* and *t* + 1 are assumed to be part of a same SSE. SSEs that were initially spatially separated but merged at some point in time are assumed to be part of a same SSE. Similarly, a single SSE can divide into multiple SSEs at some point in time but they will actually also be assumed to be part of the same SSE.

Initially detected SSEs within less than $T_{win}$ apart in time are assumed as one SSE and are combined. SSEs shorter than 1 day and that affect only one node of the subduction are removed. Remaining SSEs predicting surface displacements lower than $2 STD_{disp}$ are also removed, where $STD_{disp}$ is the standard deviation of the residual displacement between the ICA time series and their smoothed version.

For each fault node, we select the periods between each SSE and estimate a linear trend valid for this inter-SSE period (Figure 6c). We remove this estimated trend from the slip time series. In doing so, the corrected slip time series do not accumulate any slip over the overall defined inter-SSE periods. Resulting slip rates are then relative
to an inter-SSE slip rate and the detection threshold on the fault is spatially homogeneous. We repeat these steps (SSE detection and inter-SSE secular trend correction) three times and converge toward a stable solution (Figure S10 in Supporting Information S1).

We test multiple values of $T_{win}$ and $\dot{\delta}_{thresh}$ but only show here the results using $T_{win} = 30$ days and $\dot{\delta}_{thresh} = 75$ mm/yr (Figures 7–9). The threshold value is slightly above the maximum slip rate during the 2012–2013 and 2018–2019 periods (~68 mm/yr), periods which seem quiescent. This choice of parameters leads to a catalog of 14 SSEs before (events A1 to A14; Figure 7) and 13 SSEs after the Kaikoura earthquake (events B1 to B13; Figure 8), therefore 27 events in total, with a spatial separation between deep (>20 km depth) and shallow (<20 km depth) events as documented in previous studies (Wallace, 2020, and references therein). SSEs magnitudes range from $M_w$ 5.5 to 7.0 for the shallow events and from $M_w$ 6.0 to 7.5 for the deeper ones. Their duration last from 25 to 216 days for the shallow ones and 35–1005 days for the deeper ones. Note that the durations indicated here are affected by the temporal smoothing applied to retrieve the slip rate of the fault (see Section 7.3).

## 5. Inter-SSE Coupling Map Bayesian Inversion

In order to build a complete budget of slip along the megathrust, we evaluate the inter-SSE coupling map of the Hikurangi subduction zone. We use a fully Bayesian inversion procedure experienced in several earlier studies to derive the distribution of fault coupling along the megathrust and associated uncertainties (e.g., Jolivet et al., 2020). Combined with the occurrence of SSEs, we will hence be in a position to derive a probabilistic budget of slip over the 14 years covered by our data set.

Before any inversion, we need to estimate as accurately as possible the velocity field related to inter-SSE loading. We estimated in Section 2 the velocity fields of both data sets A and B from the displacement time series detrending procedure. These two velocity fields must be corrected to agree with the refined detection of SSEs. For each data set, we first reference the velocity field relative to a fixed Australian block (see Figures 10a and 10e; Text S6; Wallace, Barnes, et al., 2012) and remove horizontal deformation from the block model of Wallace, Barnes, et al. (2012) (Figure 10b). We additionally correct the velocity field from a trend estimated from the reconstructed time series using only components interpreted as noise (Table S6 and Figure S11 in Supporting Information S1). The resulting velocity field can thereafter be assumed related to tectonic sources. We then predict the GNSS velocity field corresponding to the iterative slip corrections made in the SSE detection procedure and subtract it to the Australia fixed, block motion corrected velocity field (see Section 4.2; Figures 10c and 10g). For the velocity field of data set B, we further correct for post-seismic processes due to the Kaikoura earthquake. We estimate the crustal post-seismic velocity field using the predicted surface displacement due to the Kaikoura faults from the inversion in Section 4.1. We then average the velocities from data set A and B (Figure 11b), weighting each data set by their total duration (about 8 and 6 years, respectively) to obtain a final, inter-SSE velocity field representative of secular ground motion over 14 years.

We use all three components of these GNSS-derived velocities, east, north and vertical, to derive the posterior PDF of coupling given the data. This posterior PDF is a mapping of the measurement and associated uncertainties in the space of possible coupling models accounting for our prior knowledge of the coupling distribution (Tarantola, 2005). This posterior PDF writes as

$$p(m|d) \propto p(m) \exp \left[ -\frac{1}{2}(d - Gm)^T C_{\chi}^{-1} (d - Gm) \right], \quad (5)$$

where $p(m)$ is the model prior PDF and $C_{\chi}$ the covariance matrix designed as the sum of the data covariance matrix, $C_d$, and the prediction error matrix, $C_p$, which describes uncertainties on the elastic model (Duputel et al., 2014). The elastic model used here is a 1D layered model simplified from the 3D model of Eberhart-Phillips et al. (2020) (Text S5 and Table S4 in Supporting Information S1). We assume 10% of uncertainty on the elastic properties of each layer. Uncertainties on the velocity field are discussed in Text S7 in Supporting Information S1 and Section 7.3.1 (Figure S14 in Supporting Information S1), and set here to 1.9 mm/yr for the horizontal components and 2.5 mm/yr for the vertical one. We include in the modeling a Mogi source to model an apparent subsidence signal (~17 mm/yr) at the location of the Rotorua Caldera. Since we do not aim to constrain faithfully the source of the subsidence, we simply place the Mogi source at 20 km depth which corresponds roughly to the spatial length scale of the signal. The direction of slip on the subduction interface is aligned with the direction of
<figure><img src="image_9.png" /><figcaption>Figure 7.</figcaption></figure>
long-term slip which varies along the megathrust. The long-term slip rates amplitudes and directions on the fault interface, necessary to compute the coupling, are taken from Wallace, Barnes, et al. (2012). The a priori knowledge of coupling on the fault is a uniform distribution between 0 (i.e., fully creeping) and 1 (i.e., fully locked). We sample the posterior PDF of inter-SSE coupling defined above using AlTar2, a cascading Metropolis algorithm based on the Catmip algorithm (Minson et al., 2013), to obtain an ensemble of 250,000 models.

Figure 11a shows the mean of the 250,000 samples, a faithful approximation of the mean of the posterior PDF of coupling. This model shows that the interface is mostly locked at the location of SSEs but also up-dip of the deep SSEs (Figure 11). As described in previous studies, the down-dip extent of the locking zone drops from ~20 to ~60 km depth south of latitude ~39.5°S (e.g., Maubant et al., 2023; Wallace, Beavan, et al., 2012). We also observe a low coupling separating the shallow and deep SSEs at latitude 40°S (i.e., point 5 in Figure 11a). Similarly, a decrease in coupling is suggested in the middle of a locked area at longitude 175°E and latitude 41°S (i.e., point 6 in Figure 11a). Such intriguing features have although large uncertainties considering the posterior uncertainties (Figures 11j and 11k). The inter-SSE moment deficit rate (MDR) of the whole subduction derived from this model is equal to 4.1 ± 0.2 10¹⁹ N.m/yr (inset in Figure 1).

Figure 11f to m show coupling uncertainties at different points of the fault interface. Location with low coupling in the vicinity of points 1 and 5 are relatively robust. Locked areas at points 2, 7 and 8 can also be considered as robust. Coupling inferred at point 3 shows larger uncertainties but the overall shape of the marginal posterior excludes models with coupling lower than 0.5. To the contrary, coupling at point 4 is very uncertain. Near the trench, mean coupling is close to 0.5, which is expected since this part of the megathrust is far from any GNSS station. Mean coupling at the trench therefore only reflects the average value of a uniform distribution between 0 and 1. The uncertainty of the inter-SSE coupling is illustrated in Figure 12 which shows the maps of 16 and 84 percentile of the coupling (i.e., equivalent to the standard deviation extent of Gaussian uncertainties).

The fit of the mean inter-SSE coupling model to the data is shown in Figure 11b–11e. Horizontal residuals are on average about 2.4 mm/yr (2.0 mm/yr for the vertical component). While the fixed rake constraint may explain part of the observed misfit, surface displacements are clearly not well predicted by the mean model in several places. Large residuals are observed in the vicinity of the faults that ruptured during the Kaikoura earthquake. Potential bias may arise in the local correction of our velocity field from the effect of the post-seismic processes following the Kaikoura earthquake. Any interpretation of the distribution of coupling near this location should be taken with caution. Additionally, we observe large residuals near the Rotorua Caldera, where the single Mogi source used in the inversion cannot capture the full complexity of surface deformation at this location.

The mode of the inter-SSE coupling map, shown in Figure 13e, highlights the locked areas already visible in the mean inter-SSE coupling map.

## 6. SSEs Cumulative Slip Map Bayesian Inversion

In addition to the time dependent evolution of slip during SSEs along the interface and to the distribution of coupled asperities in between SSEs, we here evaluate the total slip released by SSEs between the beginning and the end of our observation period. This slip distribution should equal to the final cumulative slip obtained in our time dependent slip inversion. However, since we could not enforce the non-negativity of slip for the time dependent inversion due to the nature of the independent components, we here aim to faithfully derive the total aseismic slip released over the period with associated uncertainties and compare them with our previous results.

We evaluate the surface displacement field due to the SSEs cumulative slip as the difference between the beginning and the end of each time series corrected from the inter-SSE trend evaluated earlier. The Bayesian approach used for the slip inversion is similar to the one in Section 5 except that we do not include a Mogi source as the apparent subsidence observed at Rotorua Caldera only affects the secular displacement rates. The a priori

**Figure 7.** Slow slip event (SSE) catalog during the pre-Kaikoura time period (2009 to Kaikoura earthquake; data set A). The top panel for each event shows its moment rate function (red curve). The Y-axis have always the same extent, which is not the case for the X-axis. The bottom panel shows the SSEs cumulative slip pattern. The slip distributions are normalized by their maximum slip which is given at the top right of the panels. The events duration, start and end, which estimates are impacted by temporal smoothing, are indicated at the bottom right of the panels. Events are ordered by the timing of their beginning. Note that multiple events can occur simultaneously. The gray contour shows portions of the fault with Restitution values over 0.8 and is used as is an indicator of well recovered inverted slip (Radiguet et al., 2011).
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Same as Figure 7 but for the post-Kaikoura time period (Kaikoura earthquake to 2023; data set B).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Slow slip event (SSE) history. (a), (b), and (c) Show the cumulative slip due to SSEs before (8 yrs; data set A), after (6 yrs; data set B) the Kaikoura earthquake, and both periods combined (14 yrs; data set A + B), respectively. The gray dashed lines in panel (a) indicate the 20 and 60 km iso-depths (C. A. Williams et al., 2013) which are used to separate the shallow (<20 km depth) and deep (between 20 and 60 km depth) sections. (d) and (e) Show the maximum slip rate evolution of the shallow (in blue) and deep (in red) sections, respectively, over the period of both data sets as a function of the latitude as shown on the maps in panel (a) and (b). (f) Shows the cumulative slip evolution at the location of the colored dots in panel (c). The vertical black dashed line in panel (d–f) indicates the date of the M7.8 Kaikoura earthquake. The white stars in panels (a–e) indicate the $M > 6.5$ earthquakes that occurred within the time period associated with each panel (U.S. Geological Survey, 2017). The green lines in panel (b) and (c) represent the surface rupture of the Kaikoura earthquake (Zinke et al., 2019).</figcaption></figure>

PDF of the total slip on the fault is set to a uniform distribution between 0 and 1.5 times the total displacement on the fault if it was fully creeping.

The mean model of SSEs cumulative slip (Figure 13a) obviously shows slip at the same location as the least square inversion (Figure 9c). It shows also some slip at the trench although it only corresponds to the mean of the a priori PDF being so that the model is not constrained at this location. Other areas with low resolution are shown outside the green contour line in Figure 13 (see discussion in Section 7.3.2, Figure S15 and S16a in Supporting Information S1), as for example, at the bottom edge of the modeled fault between North and South Islands. Those areas of slip are zones where the slip marginal PDF is close to a uniform distribution. The mean and mode of SSEs cumulative slip models show both patches of shallow high slip at 39°S and roughly three other patches in the deeper portion, one linked with the Kaikoura post-seismic slip. Slip at latitude 40°S is less visible in the Bayesian inversion than with the least square inversion (Figure 9c).
<figure><img src="image_12.png" /><figcaption>Figure 10.</figcaption></figure>
Maximum and minimum SSEs cumulative slip maps are here illustrated in Figures 13g and 13j by showing the average slip models of the 5% models with the lowest and highest moment released, respectively.

# 7. Discussion

## 7.1. The Behavior of Shallow and Deep SSEs

As shown in previous studies (Wallace, 2020, and references therein), shallow events are of shorter duration but also more frequent (20 events in 14 years, according to detections) than their deeper relatives (7 events). Extended periods of quiescence in the shallow section occur during the years 2012 and 2017 (Figure 9d). We observe that some events propagate along strike (Figure 14) as mentioned in other studies (Wallace et al., 2012). Even if these events do not necessarily fully overlap, they cover a region larger than that of each individual rupture. Hence, it is difficult to isolate distinct asperities in the shallow section. In fact, we do not observe any systematics in the spatial organization of the shallow SSEs with small (e.g., events A4 and A9 with maximum slip of ~10 and 20 mm) and large events (e.g., event A1 and A10 with maximum slip of ~100 mm) that may happen at the same place.

We observe several changes in the behavior of the shallow SSEs following the Kaikoura earthquake. To visualize it, we show in Figure 15 the difference in the frequency content of the slip rate before and after the Kaikoura earthquake, at different location of the fault. To first order, the power spectrum decreases in the shallow section after the Kaikoura earthquake south of latitude 39.5°S for frequencies below 3 months, while it increases north of this latitude. Additionally, shallow events seem able to occur or propagate toward the south (south of 40°S) more often before the Kaikoura earthquake than after. We observe a clear shallow SSE that directly follows the Kaikoura earthquake (event B2; Figure 8), as observed in Wallace et al. (2018). This SSE affects mainly the portion of the fault south of 40°S and the following SSE at this location occurs only 4 years later (event B8; Figure 8). Additionally, the average moment release rate along the shallow section of the megathrust is equal to 7.7 and 7.4 10¹⁸ N.m/yr before and after Kaikoura, respectively, thus they are quite similar. However, if one removes the contribution from the SSEs triggered by the Kaikoura earthquake by starting the evaluation after 2017.5, the average moment release rate drops to 4.9 10¹⁸ N.m/yr, which would suggest a lower value for the period after the Kaikoura earthquake. We could consider that regional changes in the ambient stress following the Kaikoura earthquake led to these changes but this would require further investigation.

Finally, the deep SSEs can be divided roughly into 2 asperities located north and south of latitude 40°S, respectively called the Kapiti and Manawatu asperities (as defined in Wallace, 2020). The deep SSE in 2013–2016 that ruptured both asperities shows a complex moment rate function (e.g., event A7; Figures 7 and 9e). While we observe some continuity in time and space during the event, it is not excluded that this event can be divided into smaller events (Frank et al., 2018; Jolivet & Frank, 2020). In general, these deep events are of longer duration and repeat less often than their shallow counterparts. Further investigations could consider the difference in terms of incoming sediment thickness and roughness to explain this change in spatial distribution and depth of occurrence for these shallow and deep SSEs (Gase et al., 2022).

Whether there is some temporal correlation between the shallow and deep SSEs is unclear (Figure 9). Before the Kaikoura earthquake, shallow SSEs seem more frequent when deep SSEs occur. But it doesn't seem the case after the Kaikoura earthquake. Additionally, the deep SSEs in 2010–2011 seems to roughly connect spatially and

**Figure 10.** (a) and (b) Show the velocity fields estimated from detrending the GNSS time series for data set A and B, respectively. The velocity field in (a) is associated to the detrending applied on the 2009–2013 time period and is used to evaluate the velocity field for the pre-Kaikoura period (data set A; see Table 1 and Figure 2). The velocity field in (b) is associated to the detrending applied on the time period between 2009 and the Kaikoura earthquake, and is used to evaluate the velocity field for the post-Kaikoura period (data set B). Both velocity fields in (a) and (b) are referenced relative to a fixed Australian block. (c) Velocity field due to internal crustal deformation based on the block model from Wallace, Barnes, et al. (2012) and referenced relative to a fixed Australian block. (d) and (e) Show the correction term applied so that the corresponding velocity fields are relative to an inter-SSE trend. (f) Correction term to remove the effect of the Kaikoura post-seismic slip occurring on the crustal faults. It is only applied to Data set B. (g) and (h) Show the final velocity fields for the pre- and post-Kaikoura period (i.e. data set A and B). In short, Final velocity field = Velocity field relative to Australia—Block model correction term—Inter-SSE correction term (– Kaikoura post-seismic). The scale for each panel is the same and is indicated by the arrow legend in panel (a) and the colorbar beside panel (b). The maximum velocity for the horizontal and vertical components of each panel, $V_{Hz}^{Max}$ and $V_{Vert}^{Max}$, respectively, are indicated at their bottom right. Enlarged panels are shown in Figures S12 and S13 in Supporting Information S1.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) Mean inter-slow slip event (SSE) coupling map. (b) Horizontal inter-SSE velocity field used to invert for the inter-SSE coupling map (black arrows) and the horizontal predicted velocity field from the mean inter-SSE coupling map (red arrows). (c) Residuals between the inter-SSE velocity field and the one predicted from the mean inter-SSE coupling map. The blue and red arrows correspond to stations with a velocity smaller and larger than twice data error (i.e., 1.9 mm/yr for the horizontal component), respectively. (d) and (e) Are the same as (b) and (c) but for the vertical component with a data error of 2.5 mm/yr. All panels with velocity fields (i.e., (b) to (e)) have the same scaling, and the length of the arrow legend in each panel is set equal to the maximum velocity present within the panel. (f–m) Show the pdf of coupling at 8 points on the fault as located in panel (a), and the black line depicts the mean value.</figcaption></figure>

temporally to the onset of a shallow SSE. Note that due to the inversions' spatial smoothing, certain SSEs might smudge through the spatial delimitation between shallow and deep SSEs.

## 7.2. Long Term Coupling and Local Slip Deficit Slip Buildup Rate

With the inter-SSE slip deficit rate and the total slip released over the period (Figures 12a and 12b, respectively), we can evaluate the total slip deficit over the period and provide an evaluation of the parts of the megathrust that have accumulated or lost elastic energy over time. We estimate the long-term slip deficit rate as the difference
<figure><img src="image_14.png" /><figcaption>**Figure 12.** (a) 16 percentile inter-slow slip event (SSE) coupling map. (b) 84 percentile inter-SSE coupling map.</figcaption></figure>

between the inter-SSE slip deficit rate (Section 5) and the SSEs long-term slip rate (Section 6). Compared to the local expected convergence rate, we infer hence the long-term coupling model. We thus have an ensemble of 250,000 long-term coupling models (Figures 13c, 13f, 13i, and 13l). The resolved area for long-term coupling (within the green contour in Figure 13; see Section 7.3.2 and Figures S15 and S16 in Supporting Information S1) is assumed as the intersect between the resolved areas from the inter-SSE coupling and SSE cumulative slip models.

Figure 13c shows the mean model. We can see that coupling ranges from -1 (more slip has been released by SSEs than should have accumulated over the period) to 1 (the area is coupled, regardless of any SSE over the period). Negative long-term coupling areas are located at the location of the Kaikoura post-seismic slip, but also where the resolution is poor such as near the trench (outside the green contour in Figure 13c). Such negative long-term coupling can be explained by the combination of priors where data has no influence on the model (Section 6).

Where long-term coupling is close to zero, the moment budget can be assumed balanced over the considered period. Where we have detected SSEs, the long-term coupling is low (≤0.4) suggesting that the slip budget is balanced. We infer areas with high coupling (>0.4), mostly up-dip of the deep SSEs, between latitudes 42° and 40°S. The mean and related standard deviation of the MDR in this area is equal to 4.9 ± 0.5 10¹⁸ N.m/yr, equivalent to a $M_w$ 6.4 per year. Other coupled areas (>0.4) along the shallow portion of the megathrust in the north, between latitude 39° and 38°S, are also observed but less certain (Figure 13i). Our analysis confirms that a large locked section of the megathrust accumulates slip deficit despite the occurrence of SSEs beneath the peninsula hosting Wellington. Such asperity clearly stands out of uncertainties.

## 7.3. Potential Bias and Limitations

Several biases might affect our budget of slip and the time dependent catalog of SSEs we have derived along the Hikurangi subduction interface. We discuss here the influence of our modeling choices and inherent limitations of our data set and methods.

### 7.3.1. Velocity and Displacement Field Uncertainties

The uncertainties on the velocity and displacement fields used to invert the inter-SSE coupling map and SSE cumulative slip, respectively, are difficult to evaluate. First, the uncertainty on the trend calculated from the detrending procedure (Figures 10a and 10b) are very likely underestimated as the trajectory model used here in Section 2 does not model SSEs. Typical uncertainties derived from a trajectory model reach values lower than 1 mm/yr. Second, the uncertainties on the velocities predicted by the block model we use to reference our velocity field reach 6 mm/yr (Wallace, Barnes, et al., 2012; Figure 10b). However, since these uncertainties derive from a
<figure><img src="image_15.png" /><figcaption>**Figure 13.** (a) and (d) Correspond to mean model and mode of the slow slip event cumulative slip. (g) and (j) Correspond to the maps of average cumulative slip of the 5% models with lowest and highest total moment release, respectively. (b) and (e) Correspond to mean model and mode of the inter-SSE coupling. (h) and (k) Correspond to the maps of average inter-SSE coupling of the 5% models with lowest and highest moment deficit rate (MDR), respectively. (c) and (f) Correspond to mean model and mode of the long-term coupling. (i) and (l) Correspond to the maps of average long-term coupling of the 5% models with lowest and highest MDR, respectively. The green curve encompasses the regions for which the posterior PDF is substantially different from its prior (see Section 7.3.2, and Figures S15 and S16 in Supporting Information S1), and thus indicates the regions that are interpreted as resolved.</figcaption></figure>

block model, they do not consider vertical motion, hence only the horizontal component of the velocity field have their uncertainties inflated by this correction. Therefore, complete propagation of uncertainties from the data to the velocity and displacement estimates would lead to underestimated uncertainties on the vertical component (<1 mm/yr) and overestimated ones on the horizontal components (>6 mm/yr).

In order to include meaningful and representative uncertainties, we isolate in both data sets periods that are long enough (>6 months) during which we do not detect any transient deformation due to SSEs. Over those periods, we compute the effect of our inter-SSE slip correction term (as in Section 5) on each component of the position time
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Comparison between the slip rate evolution of three slow slip events and their associated seismicity. (a), (b), and (c) are associated to events A5, A10, and B5 (Figures 7 and 8). Brown dots indicate tremors (Romanet & Ide, 2019; Todd & Schwartz, 2016), brown stars indicate earthquake swarms (Nishikawa et al., 2021; Yarce et al., 2019) and white stars indicate regular earthquakes (U.S. Geological Survey, 2017).</figcaption></figure>

series (i.e., East, North, and Up) for each GNSS site. The apparent velocity of the corrected time series at those periods is expected to be zero but a dispersion is observed and used as uncertainty in our slip and slip rate inversion since it corresponds to the effective disagreement between our models and data. We therefore set the uncertainties to 1.9 and 2.5 mm/yr for the horizontal and vertical components of the velocity field (Text S6 in Supporting Information S1), respectively, and 26.6 and 35 mm for the horizontal and vertical components of the displacement field (i.e., 1.9 and 2.5 mm/yr multiplied by 14 years, the observational period). We also show in the supplement the total slip budget we obtain while applying uncertainties twice larger and half smaller (Figures S17 and S18 in Supporting Information S1), indicating that the main features of our coupling and total slip distributions remain.

Note that other sources of uncertainties might persist, such as local sources of crustal deformation or viscoelastic processes due to previous large events for instance (Li et al., 2015).

### 7.3.2. Uncertainties on the Slip Budget

To illustrate the uncertainty of our slip budget in general, we show the average coupling of the 5% models with the lowest and highest MDR and the total slow slip model of the 5% models with the lowest and highest moment released (Figure 13). Those models are interpreted as maps of minimum and maximum possible coupling or slip given the data. In terms of inter-SSE coupling, we observe that the locked areas at the location of the deep SSEs and some areas up-dip, together with the locked areas of the northern shallow SSEs, are robust features since they appear in both cases (Figure 13h). We observe the same behavior for the resulting long-term coupling which systematically shows a large locked patch beneath the Wellington peninsula.

We provide also a metric which evaluates whether the posterior marginal PDF of inter-SSE coupling at each fault node is substantially different from the uniform prior PDF (Figure S15 in Supporting Information S1). For each node, we calculate the histogram of the posterior and prior coupling PDFs, and compute the difference between the two PDFs. Small differences between prior and posterior marginal PDFs indicate that those fault nodes are poorly constrained by the data (i.e., the posterior resembles the prior and no information came from the data). Green contours in Figure 13 show areas with such difference above 0.15 (Figure S16b in Supporting Information S1). We consider that areas outside of these contours should be taken with caution as local estimates of coupling and/or aseismic transient slip were simply not resolved.

### 7.3.3. Uncertainties on the Slip Temporal Reconstruction

The components of the decomposition we interpret as signal due to slip on the megathrust (Section 3.2) may be affected by some signals we would consider as noise. For instance, although we do not include any Mogi source accounting for the Rotorua Caldera in the SSE time dependent slip inversion, such source of deformation is needed for the secular inter-SSE coupling inversion. Therefore, this volcanic system might locally affect the reconstruction of transient signals by perturbing the temporal baseline of the signal which is not only of tectonic origin.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Difference in the frequency content of the slip rate history before (data set A) and after (data set B) the Kaikoura earthquake. (a), (b), and (c) show the cumulative slip due to slow slip events before (8 yrs; data set A), after (6 yrs; data set B) the Kaikoura earthquake, and both periods combined (14 yrs; data set A + B), respectively. The white stars indicate *M* > 6.5 earthquakes that occurred within the time period associated with each panel (U.S. Geological Survey, 2017). The green lines in panel (b) and (c) represent the surface rupture of the Kaikoura earthquake (Zinke et al., 2019) (d) and (e) correspond to the power spectrum of the nodes indicated in panel (a) (black dots) as function of their latitude corresponding to the maps. Note that the colorbars of the power spectrum are saturated. (f) Difference of power spectrum before and after the Kaikoura earthquake at the location indicated by the black dots in panel (a). Note that the temporal filter applied on the slip time series to extract the slip rate history removes information for periods below 2 months.</figcaption></figure>

The evolution of slip on the megathrust over our observation period also depends on the choice of regularization imposed in the a priori covariance matrix for the inversion of each component (Section 4.1). We have explored various regularization weights (effectively, the value of $\sigma_0$ in Equation 4) and decided to choose the value corresponding to the change in curvature in the relationship between the norm of the misfit and the norm of the Laplacian of slip along the fault. A choice leading to smoother SSEs would degrade too much the residuals while a choice for sharper features would lead to data overfitting. Within the range of acceptable values, we believe the variability in slip values actually lies within what is acceptable considering the posterior model covariance.

The characterization of post-seismic slip following the Kaikoura earthquake on crustal faults is not the purpose of this study. Therefore, we did not select in this study the GNSS stations close to this area and the post-seismic
process is here only coarsely modeled. In addition, we do not consider any potential viscous relaxation in the lower crust or mantle at this location. As a consequence, post-seismic slip on the subduction interface due to this earthquake should not be over-interpreted.

### 7.3.4. SSE Detections

The SSE catalog proposed here depends on the temporal smoothing (parameter $T_{win}$ in Section 4.2) and on the slip rate detection threshold (i.e., $\delta_{thresh}$; Section 4.2). Using different values, the procedure might detect more or less events (e.g., 21 and 18 events before and after the Kaikoura earthquake, respectively, taking $T_{win}$ = 20 days; Figures S19 and S20 in Supporting Information S1). The clearest example is the 2013 deep SSE in the current catalog (event A7 in Figure 7). Changing $T_{win}$ or $\delta_{thresh}$ either merges this event with SSEs in the shallow portion, thus creating a large and long-lasting event, or instead divides this event into multiple small ones. Such choice of smoothing factor impacts more strongly the deeper SSEs than the shallow ones since the former have lower slip rates, hence they stand out less above the threshold. We here chose to apply only one detection threshold on the whole fault but one could potentially set it to vary along depth. In short, the catalog proposed in this study does not match one by one the SSEs as described in other studies (e.g., Wallace, 2020, and references therein) and it is thus one catalog among many potential other ones (Figures S21 and S22 in Supporting Information S1).

For instance, Bartlow et al. (2014) detected 12 events between late 2009 and early 2012 while we extract only 5 events within the same period of time. Our SSE detection procedure leads to merging multiple events that Bartlow et al. (2014) have considered to be independent events (Figure S21 in Supporting Information S1). Our catalog can be potentially extended/modified if taking a smaller detection threshold or changing the amount of temporal smoothing (only used to detect SSEs), but with a risk of including artefacts.

In addition to the detection itself, the SSEs magnitude and duration estimates are also dependent on the chosen $T_{win}$ and $\delta_{thresh}$ (Section 4.2). Longer smoothing window will spread slip events temporally. Additionally, the magnitude of these events might be affected by the fact that we use here only onshore GNSS stations and slip near the trench is mostly unresolved. SSEs magnitude estimates of the shallow events are thus potentially truncated in terms of spatial pattern but also in time. For those reasons, while we show the magnitude-duration scaling law of the extracted SSE catalog in the supplement (Figure S23 and Text S8 in Supporting Information S1), we choose not to interpret it.

Aside from all these considerations on uncertainties, one intriguing feature stands out. We observe that the cumulative slip from SSEs in the shallow section overlaps with an area that is weakly coupled in between SSEs at latitude 39.5°S (Figure 1). This feature is also seen in another study that additionally uses InSAR time series (Maubant et al., 2023). This feature stands out both in the cumulative slip Bayesian inversion and in the temporal slip evolution model although it is less obvious on the Bayesian model as uncertainties are quite large. Since uncertainties are large on this part of the megathrust, this feature should be investigated more.

Finally, we used for the inversions a 1D layered elastic model. It has been shown that a 3D model will produce significantly different amplitude of slip than a 1D model (C. A. Williams & Wallace, 2018). This is likely one of the causes for the difference in amplitude for event B4 between our results (e.g., 85 mm slip) and the one from Woods et al. (2024) (e.g., ~200 mm), together with the fact that we applied different spatial smoothing and have not used offshore instrumentation. Further work on 3D Green's functions is required but beyond the scope of the present study.

## 8. Conclusion

In this study, we evaluate the full slip evolution of the Hikurangi subduction interface between 2009 and 2023 based on the vbICA decomposition technique. This slip evolution is not affected by any temporal smoothing. We use this product to derive a catalog of 27 SSEs that can be divided between deep/long-term and shallow/short-term SSEs. This catalog includes events detected in other studies (e.g., Bartlow et al., 2014; Wallace & Beavan, 2010; Wallace, Beavan, et al., 2012; Wallace et al., 2014, 2018; Yan et al., 2023) and our events characteristics compare well with previous studies. This long-term, continuous slip history represents a homogeneous product—characterized by a single method, consistent parametrization, and unified dataset—which enables to highlight the changes in fault behavior before and after the Kaikoura earthquake.
Based on the SSEs detection, we build maps of cumulative slip due to SSEs, and of the inter-SSE and long-term coupling. Those maps match to first order the results from previous recent and older studies (Maubant et al., 2023; Wallace & Beavan, 2010; Wallace et al., 2012), although the Bayesian approach used here allows to evaluate uncertainties at each point of the fault. Additional details on the inter-SSE coupling are observed such as the uncoupled areas up-dip of the deep-SSEs and would require further investigation. Over 14 years, SSEs mostly balance the moment budget where they occur, but a locked region remains up-dip of the deep-SSEs, beneath the Wellington peninsula (Figure 13c). This portion of the fault last ruptured in the year 1430–1480 (Clark et al., 2019). The MDR estimated in Section 7.2 (i.e., 4.9 ± 0.5 10¹⁸ N.m/yr) combined with this waiting time indicates enough elastic energy accumulated to fuel a seismic event between M<sub>w</sub> 8.1 and 8.3. However, such magnitude might be underestimated since the area and long-term MDR in the northern section (north of 40°S) is not considered. Historical earthquakes might have ruptured this section, such as, for instance the 1080–1135 event. Furthermore, slip deficit may accumulate near the trench where the resolution is poor, even though SSEs have been documented in the north to propagate near the trench (latitude ~39°S; Wallace et al., 2016). In our long-term coupling map (Figure 13c), moderate MDR lingers in the northern section and could potentially participate to such large event, although uncertainties are large (coupling of 0.5 ± 0.2). On the other hand, it has been suggested that the 1855 ~M<sub>w</sub> 8 Wairarapa earthquake, that occurred in the upper crust of the Wellington peninsula, also ruptured a portion of the subduction interface, as the Kaikoura earthquake (Herman et al., 2023; Manighetti et al., 2020; Rodgers & Little, 2006). Such event would then consume some slip deficit on the subduction interface. Considering the subduction is now being loaded since 1855, the estimated MDR would result in a seismic event between M<sub>w</sub> 7.8 and 7.9 if it were to occur presently.

Data Availability Statement

Data and results of this paper are freely available at (Michel et al., 2025). More recent solutions of the position
time series are distributed through (Geodesy Plotter, 2025). We acknowledge the New Zealand GeoNet pro-
gramme and its sponsors NHC, GNS Science, LINZ, NEMA and MBIE for providing data used in this study. The
RINEX files and metadata were acquired from (GNS Science, 1995, 2019).

Acknowledgments

This project has received funding from the
European Research Council (ERC) under
the European Union's Horizon 2020
research and innovation program (Grant
758210 for project Geo4D). Romain
Jolivet acknowledges funding from the
Institut Universitaire de France. A
significative portion of the analyses
reported in this paper were done using
MATLAB (The MathWorks Inc., 2023).
Python tool-boxes were also used:
PYACS + PYEQ (https://github.com/
JMNocquet/pyacs36; version 0.65.66),
CSI (http://www.geologie.ens.fr/~jolivet/
csi). Color scales of Crameri et al. (2020)
were used. We thank the anonymous
reviewers and the Editors for insightful
comments which helped improve the study
substantially. We thank also A. Gualandi,
L. Wallace and C. Williams for fruitful
discussions.

References

Aden-Antoniów, F., Frank, W. B., Chamberlain, C. J., Townsend, J., Wallace, L. M., & Bannister, S. (2024). Low-frequency earthquakes downdip of deep slow slip beneath the North Island of New Zealand. *Journal of Geophysical Research: Solid Earth*, 129(5), e2023JB027971. https://doi.org/10.1029/2023JB027971

Altamimi, Z., Métivier, L., Rebischung, P., Rouby, H., & Collilieux, X. (2017). ITRF2014 plate motion model. *Geophysical Journal International*, 209(3), 1906–1912. https://doi.org/10.1093/gji/ggx136

Bartlow, N. M., Wallace, L. M., Beavan, R. J., Bannister, S., & Segall, P. (2014). Time-dependent modeling of slow slip events and associated seismicity and tremor at the Hikurangi subduction zone, New Zealand. *Journal of Geophysical Research: Solid Earth*, 119(1), 734–753. https://doi.org/10.1002/2013JB010609

Behr, W. M., & Bürgmann, R. (2021). What's down there? The structures, materials and environment of deep-seated slow slip and tremor. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 379(2193), 20200218. https://doi.org/10.1098/rsta.2020.0218

Choudrey, R. A. (2002). *Variational methods for Bayesian independent component analysis*. University of Oxford. Retrieved from http://www.robots.ox.ac.uk/~parg/projects/ica/riz/Pubs/thesis.ps.gz

Choudrey, R. A., & Roberts, S. J. (2003). Variational mixture of Bayesian independent component analyzers. *Neural Computation*, 15(1), 213–252. https://doi.org/10.1162/089976603321043766

Clark, K., Howarth, J., Litchfield, N., Cochran, U., Turnbull, J., Dowling, L., et al. (2019). Geological evidence for past large earthquakes and tsunamis along the Hikurangi subduction margin, New Zealand. *Marine Geology*, 412, 139–172. https://doi.org/10.1016/j.margeo.2019.03.004

Crameri, F., Shephard, G. E., & Heron, P. J. (2020). The misuse of color in science communication. *Nature Communications*, 11(1), 5444. https://doi.org/10.1038/s41467-020-19160-7

Dragert, H., Wang, K., & James, T. S. (2001). A silent slip event on the deeper Cascadia subduction interface. *Science*, 292(5521), 1525–1528.
https://doi.org/10.1126/science.1060152

Duputel, Z., Agram, P. S., Simons, M., Minson, S. E., & Beck, J. L. (2014). Accounting for prediction uncertainty when inferring subsurface fault slip. *Geophysical Journal International*, 197(1), 464–482. https://doi.org/10.1093/gji/ggt517

Eberhart-Phillips, D., Bannister, S., Reyners, M., & Henrys, S. (2020). New Zealand wide model 2.2 seismic velocity and Qs and Qp models for New Zealand [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.3779523

Fasola, S. L., Jackson, N. M., & Williams, C. A. (2023). Deep short-term slow slip and tremor in the Manawatu region, New Zealand. *Geophysical Research Letters*, 40(21), e2023GL105428. https://doi.org/10.1029/2023GL105428

Frank, W. B., Rousset, B., Lasserre, C., & Campillo, M. (2018). Revealing the cluster of slow transients behind a large slow slip event. *Science Advances*, 4(5), 1–7. https://doi.org/10.1126/sciadv.aat0661
Gase, A. C., Bangs, N. L., Van Avendonk, H. J. A., Bassett, D., & Henrys, S. A. (2022). Hikurangi megathrust slip behavior influenced by lateral variability in sediment subduction. *Geology*, **50**(10), 1145–1149. https://doi.org/10.1130/G50261.1

Geodesy Plotter. (2025). Solution NZ_GNSS solENS [Dataset]. GNSS. Retrieved from https://emea01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.poleterresolide.fr%2Fgeodesy-plotter%2F%23%2F%3bounds%3D164.78393554687503%2C-45.83645405018772%2C185.24047851562503%2C-35.362176059146805%2expand%3D1%26solution%3DNZ_GNSS solENS&data=0%7C02%7C%7Cde03626b0bc04dca760b0ddada6e97b%7C84df9e7fe9f640afb435aaaaaaaaaaaa%7C1%7C0%7C638857654562035282%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIIYiOiWlLjAuMDAwMCIsIIAiOiJXaW4zMiIsIKFOiJoiTWFpbCIlsIdUljoyFQ%3D%3D%7C0%7C%7C%7C&sdata=sl0rY51kp9v6tDFHByhD4B87y4ILw4E8G%2BGKNU%2BiAe4%3D&reserved=0

GNS Science. (1995). GeoNet aotearoa New Zealand continuous GNSS network RINEX files [Dataset]. GNS Science. https://doi.org/10.21420/RXKE-AZ44

GNS Science. (2019). GeoNet Aotearoa New Zealand stations metadata repository [Dataset]. GNS Science, GeoNet. https://doi.org/10.21420/0VY2-C144

Gualandi, A., & Liu, Z. (2021). Variational Bayesian independent component analysis for InSAR displacement time-series with application to central California, USA. *Journal of Geophysical Research: Solid Earth*, **126**(4), e2020JB020845. https://doi.org/10.1029/2020JB020845

Gualandi, A., Serpelloni, E., & Belardinelli, M. E. (2016). Blind source separation problem in GPS time series. *Journal of Geodesy*, **90**(4), 323–341. https://doi.org/10.1007/s00190-015-0875-4

Hamling, I. J. (2020). A review of the 2016 Kaikōura earthquake: Insights from the first 3 years. *Journal of the Royal Society of New Zealand*, **50**(2), 226–244. https://doi.org/10.1080/03036758.2019.1701048

Herman, M. W., Furlong, K. P., & Benz, H. M. (2023). Substantial upper plate faulting above a shallow subduction megathrust earthquake: Mechanics and implications of the surface faulting during the 2016 Kaikoura, New Zealand, Earthquake. *Tectonics*, **42**(5). https://doi.org/10.1029/2022TC007645

Herring, T. A., King, R. W., & McClusky, S. C. (2010a). GLOBK: Global Kalman filter VLBI and GPS analysis program release 10.4.

Herring, T. A., King, R. W., & McClusky, S. C. (2010b). GAMIT reference manual, GPS Analysis at MIT, Release 10.4.

Hirose, H., Hirahara, K., Kimata, F., Fujii, N., & Miyazaki, S. (1999). A slow thrust slip event following the two 1996 Hyuganada Earthquakes beneath the Bungo Channel, southwest Japan. *Geophysical Research Letters*, **26**(21), 3237–3240. https://doi.org/10.1029/1999GL010999

Jacobs, K., Savage, M., & Smith, E. (2016). Quantifying seismicity associated with slow slip events in the Hikurangi margin, New Zealand. *New Zealand Journal of Geology and Geophysics*, **59**(1), 58–69. https://doi.org/10.1080/00288306.2015.1127827

Jolivet, R., & Frank, W. B. (2020). The transient and intermittent nature of slow slip. *AGU Advances*, **1**(1), e2019AV000126. https://doi.org/10.1029/2019AV000126

Jolivet, R., Simons, M., Duputel, Z., Olive, J. A., Bhat, H. S., & Bletery, Q. (2020). Interseismic loading of subduction megathrust drives long-term uplift in northern Chile. *Geophysical Research Letters*, **47**(8), e2019GL085377. https://doi.org/10.1029/2019GL085377

Klein, E., Vigny, C., Duputel, Z., Zigone, D., Rivera, L., Ruiz, S., & Potin, B. (2023). Return of the Atacama deep slow slip event: The 5-year recurrence confirmed by continuous GPS. *Physics of the Earth and Planetary Interiors*, **334**, 106970. https://doi.org/10.1016/j.pepi.2022.106970

Larochelle, S., Gualandi, A., Chanard, K., & Avouac, J. (2018). Identification and extraction of seasonal geodetic signals due to surface load variations. *Journal of Geophysical Research: Solid Earth*, **123**(12), 2018JB016607. https://doi.org/10.1029/2018JB016607

Li, S., Moreno, M., Bedford, J., Rosenau, M., & Oncken, O. (2015). Revisiting viscoelastic effects on interseismic deformation and locking degree: A case study of the Peru-north Chile subduction zone. *Journal of Geophysical Research: Solid Earth*, **120**(6), 4522–4538. https://doi.org/10.1002/2015JB011903

Litchfield, N., Van Dissen, R., Sutherland, R., Barnes, P., Cox, S., Norris, R., et al. (2014). A model of active faulting in New Zealand. *New Zealand Journal of Geology and Geophysics*, **57**(1), 32–56. https://doi.org/10.1080/00288306.2013.854256

Manighetti, I., Perrin, C., Gaudemer, Y., Dominguez, S., Stewart, N., Malavieille, J., & Garambois, S. (2020). Repeated giant earthquakes on the Wairarapa fault, New Zealand, revealed by Lidar-based paleoseismology. *Scientific Reports*, **10**(1), 2124. https://doi.org/10.1038/s41598-020-59229-3

Maubant, L., Frank, W. B., Wallace, L. M., Williams, C. A., & Hamling, I. (2023). Imaging the spatiotemporal evolution of plate coupling with interferometric radar (InSAR) in the Hikurangi subduction zone. *Geophysical Research Letters*, **50**(19), e2023GL105388. https://doi.org/10.1029/2023GL105388

Maubant, L., Pathier, E., Daout, S., Radiguet, M., Doin, M.-P., Kazachkina, E., et al. (2020). Independent component analysis and parametric approach for source separation in InSAR time series at regional scale: Application to the 2017–2018 slow slip event in Guerrero (Mexico). *Journal of Geophysical Research: Solid Earth*, **125**(3), e2019JB018187. https://doi.org/10.1029/2019JB018187

Michel, S., Gualandi, A., & Avouac, J. (2019a). Interseismic coupling and slow slip events on the Cascadia megathrust. *Pure and Applied Geophysics*, **176**(9), 3867–3891. https://doi.org/10.1007/s00024-018-1991-x

Michel, S., Gualandi, A., & Avouac, J.-P. (2019b). Interseismic coupling and slow slip events on the Cascadia megathrust. *Pure and Applied Geophysics*, **176**(9), 3867–3891. https://doi.org/10.1007/s00024-018-1991-x

Michel, S., Jolivet, R., Klein, E., & Maubant, L. (2025). 14 years of slip on the Hikurangi subduction zone [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.14014197

Michel, S., Jolivet, R., Lengliné, O., Gualandi, A., Larochelle, S., & Gardonio, B. (2022). Searching for transient slow slips along the san Andreas fault near park field using independent component analysis. *Journal of Geophysical Research: Solid Earth*, **127**(6), e2021JB023201. https://doi.org/10.1029/2021JB023201

Minson, S. E., Simons, M., & Beck, J. L. (2013). Bayesian inversion for finite fault earthquake source models I—Theory and algorithm. *Geophysical Journal International*, **194**(3), 1701–1726. https://doi.org/10.1093/gji/ggt180

Nishikawa, T., Nishimura, T., & Okada, Y. (2021). Earthquake swarm detection along the Hikurangi trench, New Zealand: Insights into the relationship between seismicity and slow slip events. *Journal of Geophysical Research: Solid Earth*, **126**(4), e2020JB020618. https://doi.org/10.1029/2020JB020618

Okuwaki, R., Hicks, S. P., Craig, T. J., Fan, W., Goes, S., Wright, T. J., & Yagi, Y. (2021). Illuminating a contorted slab with a complex intraslab rupture evolution during the 2021 Mw 7.3 East Cape, New Zealand earthquake. *Geophysical Research Letters*, **48**(24), e2021GL095117. https://doi.org/10.1029/2021GL095117

Radiguet, M., Cotton, F., Vergnolle, M., Campillo, M., Valette, B., Kostoglakov, V., & Cotte, N. (2011). Spatial and temporal evolution of a long term slow slip event: The 2006 Guerrero Slow Slip Event. *Geophysical Journal International*, **184**(2), 816–828. https://doi.org/10.1111/j.1365-246X.2010.04866.x
Rodgers, D. W., & Little, T. A. (2006). World's largest coseismic strike-slip offset: The 1855 rupture of the Wairarapa Fault, New Zealand, and implications for displacement/length scaling of continental earthquakes. *Journal of Geophysical Research*, **111**(B12), B12408. https://doi.org/10.1029/2005JB004065

Rogers, G., & Dragert, H. (2003). Episodic tremor and slip on the Cascadia subduction zone: The chatter of silent slip. *Science*, **300**(5627), 1942–1943. https://doi.org/10.1126/science.1084783

Romanet, P., & Ide, S. (2019). Ambient tectonic tremors in Manawatu, Cape Turnagain, Marlborough, and Puysegur, New Zealand. *Earth Planets and Space*, **71**(1), 59. https://doi.org/10.1186/s40623-019-1039-1

Tarantola, A. (2005). Inverse problem theory and methods for model parameter estimation. In *Inverse problem theory and methods for model parameter estimation*. Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898717921.fm

Todd, E. K., & Schwartz, S. Y. (2016). Tectonic tremor along the northern Hikurangi margin, New Zealand, between 2010 and 2015. *Journal of Geophysical Research: Solid Earth*, **121**(12), 8706–8719. https://doi.org/10.1002/2016JB013480

U.S. Geological Survey, E. H. P. (2017). Advanced national seismic system (ANSS) comprehensive catalog of earthquake events and products: Various. https://doi.org/10.1145/annurev-earth-071719-055104

Wallace, L. M. (2020). Slow slip events in New Zealand. *Annual Review of Earth and Planetary Sciences*, **48**(1), 175–203. https://doi.org/10.1145/annurev-earth-071719-055104

Wallace, L. M., Barnes, P., Beavan, J., Van Dissen, R., Litchfield, N., Mountjoy, J., et al. (2012). The kinematics of a transition from subduction to strike-slip: An example from the central New Zealand plate boundary. *Journal of Geophysical Research*, **117**(B2), B02405. https://doi.org/10.1029/2011JB008640

Wallace, L. M., Bartlow, N., Hamling, I., & Fry, B. (2014). Quake clamps down on slow slip. *Geophysical Research Letters*, **41**(24), 8840–8846. https://doi.org/10.1002/2014GL062367

Wallace, L. M., & Beavan, J. (2010). Diverse slow slip behavior at the Hikurangi subduction margin, New Zealand. *Journal of Geophysical Research*, **115**(B12), B12402. https://doi.org/10.1029/2010JB007717

Wallace, L. M., Beavan, J., Bannister, S., & Williams, C. (2012b). Simultaneous long-term and short-term slow slip events at the Hikurangi subduction margin, New Zealand: Implications for processes that control slow slip event occurrence, duration, and migration. *Journal of Geophysical Research*, **117**(B11), B11402. https://doi.org/10.1029/2012JB009489

Wallace, L. M., Beavan, J., McCaffrey, R., & Darby, D. (2004). Subduction zone coupling and tectonic block rotations in the North Island, New Zealand. *Journal of Geophysical Research*, **109**(B12), 1–21. https://doi.org/10.1029/2004JB003241

Wallace, L. M., Hreinsdóttir, S., Ellis, S., Hamling, I., D'Anastasio, E., & Denys, P. (2018). Triggered slow slip and afterslip on the southern Hikurangi subduction zone following the Kaikōura earthquake. *Geophysical Research Letters*, **45**(10), 4710–4718. https://doi.org/10.1002/2018GL077385

Wallace, L. M., Webb, S. C., Ito, Y., Mochizuki, K., Hino, R., Henrys, S., et al. (2016). Slow slip near the trench at the Hikurangi subduction zone, New Zealand. *Science*, **352**(6286), 701–704. https://doi.org/10.1126/science.aaf2349

Williams, C. A., Eberhart-Phillips, D., Bannister, S., Barker, D. H. N., Henrys, S., Reyners, M., & Sutherland, R. (2013). Revised interface geometry for the Hikurangi subduction zone, New Zealand. *Seismological Research Letters*, **84**(6), 1066–1073. https://doi.org/10.1785/0220130035

Williams, C. A., & Wallace, L. M. (2018). The impact of realistic elastic properties on inversions of shallow subduction interface slow slip events using seafloor geodetic data. *Geophysical Research Letters*, **45**(15), 7462–7470. https://doi.org/10.1029/2018GL078042

Woods, K., Wallace, L. M., Williams, C. A., Hamling, I. J., Webb, S. C., Ito, Y., et al. (2024). Spatiotemporal evolution of slow slip events at the offshore Hikurangi subduction zone in 2019 using GNSS, InSAR, and seafloor geodetic data. *Journal of Geophysical Research: Solid Earth*, **129**(8), e2024JB029068. https://doi.org/10.1029/2024JB029068

Yan, L., Sun, Y., Li, M., El-Mowafy, A., & Ma, L. (2023). Slow slip events associated with seismic activity in the Hikurangi subduction zone, New Zealand, from 2019 to 2022. *Remote Sensing*, **15**(19), 4767. https://doi.org/10.3390/rs15194767

Yarce, J., Sheehan, A. F., Nakai, J. S., Schwartz, S. Y., Mochizuki, K., Savage, M. K., et al. (2019). Seismicity at the northern Hikurangi margin, New Zealand, and investigation of the potential spatial and temporal relationships with a shallow slow slip event. *Journal of Geophysical Research: Solid Earth*, **124**(5), 4751–4766. https://doi.org/10.1029/2018JB017211

Zinke, R., Hollingsworth, J., Dolan, J. F., & Van Dissen, R. (2019). Three-dimensional surface deformation in the 2016 Mw 7.8 kaikōura, New Zealand, earthquake from optical image correlation: Implications for strain localization and long-term evolution of the pacific-Australian plate boundary. *Geochemistry, Geophysics, Geosystems*, **20**(3), 1609–1628. https://doi.org/10.1029/2018GC007951

## References From the Supporting Information

Eberhart-Phillips, D., & Bannister, S. (2015). 3-D imaging of the northern Hikurangi subduction zone, New Zealand: Variations in subducted sediment, slab fluids and slow slip. *Geophysical Journal International*, **201**(2), 838–855. https://doi.org/10.1093/gji/ggv057

Hamling, I. J., D'Anastasio, E., Wallace, L. M., Ellis, S., Motagh, M., Samsonov, S., et al. (2014). Crustal deformation and stress transfer during a propagating earthquake sequence: The 2013 Cook Strait sequence, central New Zealand. *Journal of Geophysical Research: Solid Earth*, **119**(7), 6080–6092. https://doi.org/10.1002/2014JB011084

Ide, S., & Beroza, G. C. (2023). Slow earthquake scaling reconsidered as a boundary between distinct modes of rupture propagation. *Proceedings of the National Academy of Sciences*, **120**(32). https://doi.org/10.1073/pnas.2222102120

Ide, S., Beroza, G. C., Shelly, D. R., & Uchide, T. (2007). A scaling law for slow earthquakes. *Nature*, **447**(7140), 76–79. https://doi.org/10.1038/nature05780

Michel, S., Gualandi, A., & Avouac, J. (2019). Interseismic coupling and slow slip events on the Cascadia megathrust. *Pure and Applied Geophysics*, **176**(9), 3867–3891. https://doi.org/10.1007/s00024-018-1991-x

Pearson, K. (1901). LIII. On lines and planes of closest fit to systems of points in space. *The London, Edinburgh and Dublin Philosophical Magazine and Journal of Science*, **2**(11), 559–572. https://doi.org/10.1080/14786440109462720