Slow Slip Events in New Zealand: Irregular, yet Predictable?

Current earthquake forecasting approaches are mainly based on probabilistic assumptions, as earthquakes seem to occur randomly. Such apparent randomness can however be caused by deterministic chaos, rendering deterministic short‐term forecasts possible. Due to the short historical and instrumental record of earthquakes, chaos detection has proven challenging, but more frequently occurring slow slip events (SSE) are promising candidates to probe for determinism. Here, we characterize the SSE signatures obtained from GNSS position time series in the Hikurangi Subduction Zone (New Zealand) to investigate whether the seemingly random SSE occurrence is governed by chaotic determinism. We find evidence for deterministic chaos for stations recording shallow SSEs, suggesting that short‐term deterministic forecasting of SSEs, similar to weather forecasts, might indeed be possible over timescales of a few weeks. We anticipate that our findings could open the door for next‐generation SSE forecasting, adding new tools to existing probabilistic approaches.


Introduction
Tectonic plate movements can lead to various modes of deformation, ranging from slow slip to fast earthquakes.Accurately forecasting such slip events remains a significant challenge, with current methods primarily relying on probabilistic considerations such as frequency-magnitude distributions (Jordan et al., 2011).These types of forecasts are based on the recordings of past events, and likelihood of future earthquakes occurrence is expressed as a probability that an event with a certain magnitude will occur within a certain time and area, with the assumption of a time-independent, linear stochastic (or random) process (e.g., Gerstenberger et al., 2020), with recent models also incorporating time-dependent characteristics (e.g., Field et al., 2015;Gerstenberger et al., 2022).
A fundamentally different approach is to consider deterministic forecasts which yield a single possible outcome only.The determinism could stem from the understanding of underlying physical processes, such as assuming that earthquakes originate from frictional sliding, which is a nonlinear process (Urbakh et al., 2004).In nonlinear systems, small changes in the input can lead to disproportionate, and often non-intuitive outputs.Such deterministic nonlinearity can result in chaotic systems, which are characterized by their irregularity and apparent randomness, yet are still deterministic and predictable at short timescales (Kantz & Schreiber, 2003).If analyzed with conventional linear tools (e.g., Fourier transform), chaotic time series appear random, but are in fact full of patterns only revealed with nonlinear analysis techniques (Abarbanel, 1996).Previous studies have suggested that both slow and fast fault slip might be governed by deterministic chaos (Barbot, 2019;Gualandi et al., 2020;Huang & Turcotte, 1990a, 1990b;Kato, 2016;Nie & Barbot, 2021;Poulet et al., 2014;Shelly, 2010;Veveakis et al., 2017).
Due to the long recurrence times of major earthquakes, data to identify deterministic behavior for seismic slip is currently insufficient.However, the frequently occurring slow slip events (SSE) with long rise times and slow rupture velocities do not suffer from this limitation and provide an opportunity for the detection of determinism.Land surface displacement during SSEs has been detected over the last two decades with continuously operating GNSS at many plate boundaries.In some cases, such as Cascadia and New Zealand, these networks have been operating for 20 years or more, and have recorded more than a dozen slow slip cycles (Schmidt & Gao, 2010;Szeliga et al., 2008;Wallace, 2020;Wallace & Beavan, 2010).
In the North Island of New Zealand, SSEs occur frequently but irregularly on the Hikurangi subduction zone (SZ) interface (Figure 1a).A large variety of events have been detected by the continuously operating GNSS monitoring network of GeoNet (www.geonet.org.nz;Wallace, 2020;Wallace & Beavan, 2010).Shallow SSEs (<15 km depth) that are short in duration (<1 month) and occur frequently (every 1-2 years) are observed offshore the northern and central East Coast.At the southern Hikurangi subduction zone, where most of the interface appears to be interseismically locked to ∼25-30 km depth, SSEs occur at greater depths (>30 km), last longer (>1 year), and seem less frequent (ca.every 5 years).
In this study, we aim to explore the presence of deterministic chaos in the continuous GNSS time series recording SSEs in the Hikurangi SZ using nonlinear analysis techniques.If SSEs indeed exhibit characteristics of deterministic chaos, it would enable to go beyond probabilistic approaches and predict both timing and magnitude of SSEs at the short term.

Phase Space Reconstruction
We use the daily position timeseries data of various GNSS stations as the surface manifestation of the spatiotemporal evolution of slip events at depth, to investigate the dynamical properties of the SSEs (Figure 1) (GNS  (GNS Science, 2000).SSEs are identified as distinctive jumps in the GNSS position timeseries, followed by slow westward drifting Inter-SSE periods.Black lines in the bottom of each subfigure show the derived daily velocities.Science, 2000).We therefore employ low-pass Finite Impulse Response filtering using a Hamming window with a cut-off frequency f c of 1/50 and a window length L of 60 days to reduce measurement noise (see Text S1 in Supporting Information S1).We then derive the position differences, yielding daily velocities v(t) (Figure 2a), and use singular system analysis (SSA) with an SSA window width M of 60, to decompose the signal into different singular vectors (SV) (Figure 2b).The values t i of these vectors at a specific time can then be used as new coordinates to reconstruct the multidimensional phase space of the SSE system (Broomhead & King, 1986;King et al., 1987) (see Text S2 in Supporting Information S1) (Figure 2c).The phase space describes all possible states of a dynamical system, with each state t i defined as a point, and each degree of freedom represented by an axis.The temporal sequence of all states forms a trajectory, along which the system evolves in time.If the trajectories are structured in phase space, they form an attractor that is an invariant of the system and can reveal otherwise undetected patterns (Abarbanel, 1996).Since the reconstructed attractor represents a proxy of the original system, it provides the basis for the detection of deterministic chaos and forecasting with nonlinear analysis tools (Kantz & Schreiber, 2003).
Figure 2b displays the three first left singular vectors (SV) obtained through the SSA of v(t) from GNSS station GISB.The phase space, reconstructed from these SVs, reveals a structured attractor (Figure 2c).This attractor comprises well-defined, loop-shaped trajectories with upward-folded sides, where each loop corresponds to an SSE.The size of each trajectory is directly proportional to the SSE's magnitude.During Inter-SSE periods, the system occupies a small subspace near the attractor node only, while significant changes in the SVs occur during SSE activity.Analyzing data from multiple GNSS stations, we found that stations clearly recording shallow SSEs have well-defined attractors with similar structures, such as GISB and PAWA (Figures 2c and 2d).The attractors of deep SSEs, exemplified by the station TAKP, exhibit a similar dynamical fingerprint as well, but with significantly fewer loops (Figure 2e).Due to their longer recurrence intervals of 4-5 years and fewer SSE cycles observed, the attractors of deep SSEs are more sparsely mapped.Therefore, we do not attempt to further quantify properties of the deep SSEs.At stations which were not significantly affected by SSEs, we were not able to obtain well-defined attractors (see RIPA, Figure 2f).
The phase space's dimensionality approximates the degrees of freedom of the underlying physical problem (Theiler, 1990), and low-dimensional spaces are characteristic for deterministic chaos (see Section 2.3).The dimensionality of the phase space, called embedding dimension d, can be estimated directly from the spectrum of singular values (Broomhead & King, 1986;Ghil et al., 2002;King et al., 1987) (see Text S2 in Supporting Information S1).We find dimensionalities d of approximately six for shallow and 35 for deep SSEs (Figure S2 in Supporting Information S1).The significantly larger dimensionality for the deep SSEs is likely due to noise, since larger values of M decrease the significance of the individual SVs and smooth the signal (Figures S3-S5 in Supporting Information S1) (Vautard & Ghil, 1989;Walwer et al., 2016).Since d approximates a theoretical upper bound of the number of degrees of freedom (Vautard & Ghil, 1989), it provides a first-order indication of the presence of a low-dimensional system, at least for the shallow SSEs.It is important to note that d is not an invariant of the system, but depends on data quality, sampling rate, and window width (Broomhead & King, 1986;Vautard & Ghil, 1989).Varying both low-pass filtering and M does however not significantly change the values of d (Figures S6-S8 in Supporting Information S1).

Nonlinearity
The observation of irregular recurrence times and the unraveling of structure in phase space with geometrically complex attractors, consisting of many different-sized trajectories, arguably point to a chaotic system.Nonlinearity is a necessary but not sufficient condition for deterministic chaos (Kantz & Schreiber, 2003), and fortunately, the detection of nonlinearity is relatively insensitive to noise (Theiler, 1992).A popular method to test for nonlinearity in a time series from an unknown dynamical system is surrogate data hypothesis testing (Kantz & Schreiber, 2003;Lancaster et al., 2018;Theiler, 1992) (see Text S3 in Supporting Information S1).We use wavelet iterative amplitude adjusted Fourier transform (WIAAFT) surrogates to test the null hypothesis that the signals can be explained by a non-stationary linear stochastic process, potentially with a nonlinear measurement transformation (Keylock, 2006).If the values of a discriminatory statistic, here Takens' best estimate of the correlation dimension D 2 (T) , for 100 realization of surrogates differ significantly from the statistics computed from the measured time series, we reject the null hypothesis (Lancaster et al., 2018;Theiler, 1992).Since a deterministic nonlinear system supposedly leads to a lower dimension of the attractor, we reject the null hypothesis if D 2 (T) of the measured time series is smaller than the fifth percentile of the surrogates (Lancaster et al., 2018;Theiler, 1992).
Using this approach, we can reject the null hypothesis for 11 different GNSS time series, all recording shallow East Coast SSEs (Figure 3).These measurements thus cannot be explained by a linear stochastic process, and we can assume nonlinearity.Interestingly, the detected nonlinear signals appear in two spatial clusters, both located in close proximity to the regions with maximum cumulative slow slip (Figure 3a).This suggests that the shallow SSEs in the nearby source regions are nonlinear systems that might even be chaotic (Theiler, 1992).For all other stations recording only a few SSEs, including stations recording deep SSEs, we could not reject the null hypothesis, which does however not explicitly exclude nonlinearity.

Indications of Deterministic Chaos
In contrast to finding evidence for nonlinearity, detecting deterministic chaos in short and noisy time series has proven more challenging (Casdagli, 1992;Kantz & Schreiber, 2003;Theiler, 1992).The revealed structure in phase space (Figure 2), as well as the detected nonlinearity (Figure 3), gives a first hint for the presence of chaos (Abarbanel, 1996).Further invariants that are indicative of chaotic systems are finite-dimensional fractal (or strange) attractors and exponential divergence of nearby states.However, due to the noisy and limited data available, quantifying these invariants has proven difficult (Eckmann & Ruelle, 1992;Theiler, 1992).Therefore, we further analyze the data from the station GISB only, which records shallow SSEs exceptionally well and is, in our view, the most suitable site in the Hikurangi SZ for this analysis.
The qualitatively self-similar attractor shown in Figure 2c suggests a geometrically complex, potentially fractal attractor (Casdagli, 1992).Estimating the correlation dimension D 2 with an entropy-based approach (see Text S4), we were unable to identify a scaling region or derive a value of D 2 , regardless of the filters used (Figures S9-S12 in Supporting Information S1).This is not surprising, given that the measurement noise obscures the smallscale features of the attractor (Datseris & Parlitz, 2022;Kantz & Schreiber, 2003).A stochastic process should however result in a constant increase of D 2 with increasing d (Casdagli, 1992;Kantz & Schreiber, 2003).We thus argue that the saturation of D 2 with increasing d (Figures S9-S12 in Supporting Information S1) can be viewed as a first-order indication of a low-dimensional system.
Chaotic systems are governed by exponential divergence of nearby states, expressed by a positive Lyapunov exponent λ 1 (see Text S5 in Supporting Information S1) (Kantz & Schreiber, 2003).Directly deriving λ 1 from the attractors (Skokos et al., 2016), we found positive λ 1 values in the range of 8-12, dependent on the filtering parameters (Figures S13-S16 in Supporting Information S1).Due to the exponential divergence of nearby states, small variations in the initial conditions lead to vast differences in the outcome, rendering chaotic systems unpredictable at longer time-scales despite their deterministic nature.The Lyapunov time t λ , defined as the inverse of λ 1 , serves as a measure of the characteristic timescale on which a system is chaotic and provides a theoretical limit of the predictability.Values of λ 1 between 8 and 12 suggest that the shallow SSEs recorded at GISB should be predictable to within 1/8 to 1/12 of the approximate recurrence period of the shallow SSEs of about one year, implying that the onset of shallow SSEs might be forecast about 30-45 days in advance.

Discussion
Compared to the Cascadia SZ, SSEs in New Zealand occur quite irregularly.This is not surprising, given the more heterogeneous structure (Barnes et al., 2020;Gase et al., 2022) and slip behavior of the Hikurangi SZ, with substantial spatial variation in the distribution of SSEs and interseismic coupling (Wallace, 2020).Our findings imply that, despite appearing random, the SSEs in New Zealand might still be deterministic and their onset could be forecast a few weeks in advance, remarkably similar to observations from the Cascadia SZ (Gualandi et al., 2020).We argue that the irregular recurrence of SSEs in New Zealand might be governed by stress interactions between different SSE source regions, since deterministic chaos can arise from the coupling of different systems (Huang & Turcotte, 1990b).
As the amount of currently available data is limited, only a subset of the SSE attractors is mapped, and the attractors likely consist of an infinite number of different-sized trajectories, explaining the system's variability in magnitudes.This self-similarity observed at larger scales suggests that SSEs with very small magnitudes and short durations of a few days might occur during Inter-SSE periods.Such hidden SSEs are however obscured by measurement noise, and stronger filtering for noise reduction would remove the signals of such tiny SSEs as well, making their identification difficult (Frank, 2016;Jolivet & Frank, 2020;Kato & Nakagawa, 2020).We suggest that SSEs in New Zealand occur more frequently than observed with the currently available data.
The existence of deterministic chaos opens the door for next-generation SSE forecasting, since the system's shortterm behavior can directly be derived from the phase space, for example, using local modeling (Abarbanel, 1996;Engster & Parlitz, 2006) (see Text S6 in Supporting Information S1).To illustrate this conceptually, we use the GNSS data from 2002 to September 2020 of station GISB as training data to reconstruct the phase space (Figure 4a).The trajectories closest to the most recent measurement then approximate the future evolution of the system in phase space (Figure 4a), which can be converted back to the measured GNSS position time series (Figure 4b).In the case of the example in Figure 4, the short-term forecast of the daily positions appears relatively accurate.Beyond the theoretical forecast limit, the chaotic nature of the system renders long-term forecasts impossible, which is illustrated by the mismatch of the SSE forecast in late 2021.However, the theoretical limit derived from λ 1 is a rather conservative estimate of the system's predictability (Datseris & Parlitz, 2022).For accurate short-term forecasts, a good description of the phase space is required, which is not currently fulfilled with the limited amount of data available, which maps the phase space sparsely (Figures 2c and 4a, and Figure S17 in Supporting Information S1).Potential ways to overcome these limitations could be enhanced data quality (i.e., higher signal-to-noise ratio), longer time series, or the incorporation of physics-based models.In our view, the well-recorded shallow SSEs represent the most promising systems for effective short-term SSE forecasting in New Zealand.

Conclusion
In summary, indications for nonlinearity, sensitivity to initial conditions, and self-similarity suggest that the irregularly occurring SSEs in New Zealand may exhibit deterministic chaos.We demonstrate that, despite the apparent randomness, with long and precise enough datasets, there appears to be some degree of predictability in SSEs.Similar to forecasts of chaotic weather systems, both timing and magnitude of SSEs have potential to be predicted a few weeks in advance.As SSEs lead to stress changes on adjacent parts of the subduction interface which can trigger large earthquakes (e.g., Kato et al., 2012;Koulali et al., 2017;Radiguet et al., 2016;Segou & Parsons, 2020), short-term SSE forecasting might also allow to identify phases of enhanced earthquake hazard (Obara & Kato, 2016).Evidence that also seismic slip might be chaotic (e.g., Barbot, 2019;Kato, 2016;Shelly, 2010) suggests that similar deterministic earthquake forecasts could one day be possible when long enough time series of the seismic cycle become available.

Figure 1 .
Figure 1.(a) Map of the GNSS stations of the GeoNet monitoring network (www.geonet.org.nz;blue dots) and cumulative slow slip (2002-2014; red contours in mm) recorded in the Hikurangi subduction zone (after Wallace, 2020).Gray dashed lines represent depth contours of the subduction interface (in km below sea level; data from Williams et al., 2013).(b) Time series of the daily displacement signals (blue: data; red: low-pass filtered data) measured at four different GNSS stations (all E component) (GNS Science, 2000).SSEs are identified as distinctive jumps in the GNSS position timeseries, followed by slow westward drifting Inter-SSE periods.Black lines in the bottom of each subfigure show the derived daily velocities.

Figure 2 .
Figure 2. (a) Time series of measured displacements (blue; filtered signal in red) and daily velocities (black).(b) The three first singular vectors (SV) derived with the singular system analysis (SSA) of the daily velocities for the station GISB.(c) Reconstruction of the phase space for data from the station GISB, using the first three SVs as coordinates.(d)-(f) Examples of the phase space for shallow (station PAWA), deep (station TAKP), and no obvious SSE signals (station RIPA).

Figure 3 .
Figure 3. (a) Map of the GNSS stations, colored after the results of the surrogate data hypothesis tests.Stations with detected nonlinearity (red squares) cluster in close proximity to locations of largest onshore displacement during shallow SSEs offshore the east coast.We note that the stations recording deep SSEs only incorporate a few SSE cycles (due to their longer recurrence intervals), so they do not record a sufficient number of SSEs to assess nonlinearity (though this cannot be ruled-out).(b) Boxwhisker plots of the surrogate data test for a selection of stations along the east coast.The whiskers indicate the fifth, and 95th percentiles of the surrogate distributions.Blue diamonds show Takens' best estimate of the correlation dimension D 2 (T) of the measured time series.

Figure 4 .
Figure 4. Conceptual example of local modeling-based short-term SSE forecast for the station GISB, using the GNSS position time series from 2002 to September 2020 as training data.(a) Training data (green to yellow) and forecast (orange) in phase space.(b) Forecast converted to the GNSS position time series (orange line), compared to the measured displacement (blue; filtered signal in red) and the derived daily velocities (black).The vertical solid line marks the end of observations and start of the forecast (current state), and the dashed white line indicates the theoretical forecast limit defined by the Lyapunov time t λ .