The 2019 East Coast Slow Slip Event, New Zealand: Spatiotemporal Evolution and Associated Seismicity

Abstract Slow slip events (SSEs) are interpreted as the transient quasi-static fault deformation in the deep transition zone from locked to freely slipping in many subduction zones. Using continuous Global Positioning System (cGPS) data collected in New Zealand, we estimate the spatiotemporal evolution model during the 2019 SSE and analyze the influence of subduction interface heterogeneity on seismicity during SSEs at the Hikurangi margin. The results reveal that the 2019 SSE extends from the northern (Gisborne) to the central (Hawke’s Bay) Hikurangi subduction interface and decays rapidly within approximately 3-4 weeks. It releases a total seismic moment of about 4.83 × 1019 N·m (Mw 6.8), with a significant slip in Gisborne and a secondary slip in Hawke’s Bay. The slip depths are similar, but peaks, durations, and rates differ slightly. By combining previous SSEs (2011-2019), diverse characteristics are concluded, i.e., shorter duration and more frequency in Gisborne and relatively longer duration and less frequency in Hawke’s Bay. The seismicity offshore and onshore indicates along-strike variations, which appear to be spatially correlated with the variations in topography, such as subduction seamounts. The heterogeneities on the subduction interface are related to the spatiotemporal distribution of SSEs and seismicity along the Hikurangi margin.


Introduction
Numerous slow slip events (SSEs) have been recorded by continuous Global Positioning System (cGPS) time series at subduction margins, including Cascadia, Nankai, and Nicoya (Hirose et al. 1999;Obara and Hirose 2006;Shelly et al. 2006;Kano and Kato 2020). They occur in areas representing a transition zone between a locked seismogenic zone in shallow regions where large earthquakes nucleate and an aseismic creep zone in deep regions along the plate interface (B€ urgmann 2018). A widely accepted view is that the effective normal stress is reduced by the low frictional resistance of fault rocks and high pore pressure along the subduction interface and contributes to the occurrence of SSEs (Liu and Rice 2007;Skarbek, Rempel, and Schmidt 2012). Numerous SSEs are considered to be a type of slow earthquake, which periodically releases accumulated strain that can trigger microseismic sequences (Beroza and Ide 2011;Yarce et al. 2019) and/or interplate earthquakes (Obara and Kato, 2004;Ito et al. 2013;Radiguet, Perfettini, and Cotte 2016;Husker et al. 2019). In addition, SSEs are not inherently hazardous but can change the stress field, which could further induce potentially damaging earthquakes along the plate interface or elsewhere (Bekaert, Hooper, and Wright 2015;Obara and Kato 2016). Several ideas have been presented to explain the physical characteristics of the transient slow slip behavior, however, there is still no consensus (Shibazaki and Iio 2003). Therefore, it is essential to investigate the generation and evolution of SSEs to assess potential earthquake-tsunami hazards.
In the past few decades, more than 20 SSEs have been observed along the northern Hikurangi subduction interface, as illustrated in Figure 1 Wallace et al. 2012;Wallace et al. 2014;Wallace et al. 2018;Wallace et al. 2016). These SSEs are characterized by shallow slip (<15 km), short duration (lasting a few weeks), and relatively frequent recurrence intervals (1-2 years) Wallace et al. 2012;Todd et al. 2018;Wallace 2020). Some studies suggested that shallow SSEs can be related to the seismicity rate at Hikurangi (Delahaye et al. 2009;Wallace et al. 2012; Bartlow et al. 2014). Numerous SSEs affect or trigger large intra-slab earthquakes . Some theoretical and observational evidence also demonstrated that local and/or regional transient seismic events influence SSEs and associated seismic sequences or suppress ongoing SSEs, such as the 2014 Mw 6.3 Eketahuna and 2016 Mw 7.8 Kaikoura events Wallace et al. 2017). In addition, seismicity along the Hikurangi margin is concurrent with migrating SSEs, with most seismic events occurring within the subducting Pacific plate and fewer at the plate interface Yarce et al. 2019). At present, the relationship between SSEs and seismicity at the Hikurangi subduction margin is poorly understood. Here, we retrieve the spatiotemporal evolution of the east coast SSEs along the Hikurangi subduction interface using geodetic data, with an objective to further investigate the relationship between SSEs and seismicity.
This study uses the principal component analysis-based inversion method (PCAIM) (Kositsky and Avouac 2010) to investigate the detailed characteristics of the spatiotemporal evolution of the 2019 SSE in the northern Hikurangi subduction zone. The obtained principal components are recombined to obtain the best-fitting time-dependent slip model, and the relationship between the time-dependent slip model and seismicity is analyzed. The effects of subduction interface heterogeneity are discussed for the SSE distribution and seismicity by integrating previous seismic events. Gisborne earthquake (https://github.com/GeoNet/data/blob/main/moment-tensor/GeoNet_CMT_ solutions.csv), and its aftershocks are shown as blue dots; the dashed gray lines indicate subduction interface contours from Williams et al. (2013). The corner inset shows the tectonic setting of the 2019 SSE along Australia (AUS) -Pacific (PAC) plate boundary zone. (b and c) Northsouth and east-west components of the cGPS detrended time series from MAHI station, showing the east coast SSE in 2019 with a gray color ramp.

Extraction of slow slip signals
The daily time series of 45 stations are adopted from a cGPS network operated by the New Zealand's GeoNet (http://www.geonet.orz.na) and PositioNZ (http://apps.linz.govt.nz/positionz) agencies (Figure 1, green dots). The time series obtained ranges from 2010.00 to 2019.91. The coordinates and their formal uncertainties were extracted from the combined daily solution and converted to three-directional (east, north, up) displacements with respect to a priori position and epoch in the ITRF2008 realization. Generally, the resulting position time series have no adjustments or slow slip signals may be offsets-contaminated due to earthquakes, equipment changes at individual sites, seasonal variations (annual and semiannual), and measurement noise. The position time series model is typically written as follows (Nikolaidis 2002;Langbein 2004): where yðt i Þ is the observed coordinate in time with t i starting from 1 January 2010; y 0 is the initial position; v is the long-term secular velocity; A, B, C, and D are amplitudes of seasonal deformation (annual and semiannual); H denotes the Heaviside step function; O j is the jth nontectonic offset associated with GPS instrument issues (update, replacement, and failure); k shows the number of nontectonic offsets; c l is the lth seismic offset; n indicates the number of seismic offsets; b p is the final displacement caused by the pth SSE with p indicates the number of SSEs; s exp is the slow slip time constant approximated by an exponential function; tq is the onset time of the SSE; e is the observational error. In this study, using an exponential function ) to model the 2019 SSE and extract the SSE signals, four primary steps are performed. (1) The outliers were eliminated, and the offset terms were estimated from the position time series. Some observations were removed when the root mean square error values exceeded 20 mm and 40 mm for the horizontal and vertical components, respectively. Besides, observations were removed by applying an iterative outlier rejection criterion when the residuals were more than 3-sigma using the normalized root mean square (Herring 2003). Offset terms associated with GPS instrument issues and large earthquakes, such as the 2016 Mw 7.8 Kaikoura earthquake, were also removed.
(2) The amplitudes of sinusoidal signals with annual and semiannual variations were estimated for each component at every GPS station.
(3) The observational epochs visibly affected by past SSEs were manually eliminated, except for the 2019 SSE ( Figure S1). Additionally, the post-seismic signals of the 2016 Kaikoura earthquake captured by the GPS sites of North Island were defined as mainly slow slip signals Jiang et al. 2018). In order to improve the signal-to-noise ratio of the 2019 SSE, we manually removed the relevant data from the position time series. Moreover, the secular velocities of GPS stations were robustly estimated by the least-square algorithm. Generally, earth surface is moving at a constant rate with respect to ITRF. In this study, the linear trend was removed from the observations when we cleaned up the time series. This step effectively made the cleaned time series comparable with the model. (4) The slow slip signals are extracted from the cleaned position time series from Figure 2, by correcting the long-term rate terms, seasonal (annual and semiannual) signals, nontectonic and seismic offsets, and outliers in the position time series. For the low signal-to-noise SSEs, it is necessary to appropriately filter the correlated common mode nontectonic signals and noise in the GPS times series. This can be achieved by filtering tools built into the spatiotemporal filtering algorithm during slow slip inversion with PCAIM in the next section.

Slow slip inversion with PCAIM
To estimate the spatiotemporal evolution of the 2019 east coast SSE, we utilized the PCAIM to perform the slip distribution inversion using the GPS position time series at the Hikurangi margin. Position time series (north, east, up) was decomposed into a linear combination of several principal components, for which each component is composed of a spatial vector, a singular value, and a time function (see details in S1). Notably, the PCAIM assumes that the SSE signal is synchronous everywhere, which is generally not the case for SSEs characterized by slip propagation. For example, in Figure 3c, the southern GPS sites show the landward motion and uplift, which is not expected from thrust faulting due to the nearby SSE patch. This normal faulting patch could cancel out the earlier onset of SSE than the actual timing due to the First PC. It is obvious because of the different onset timing of transient motion at CKID and MAHI (Figures 2  and 4). Considering the first-order features of the SSE displacement, this inversion technique is thought to be still able to effectively restore the major character of migration.
Similar to Wallace et al. (2018), a three-dimensional (3 D) fault geometry was used, in which the parameter dip was set to 12 for shallow and 23 for deep segments with a strike of 221 and rake of 120 . The "make_fault_model.m" routine of the PCAIM package was adopted to mesh the fault subduction interface with approximately regular triangular subpatches and then use the Laplacian operator to smooth the fault slip during the inversion (Perfettini and Avouac 2014). The optimal spatial smoothing factor and the number of principal components were determined using the reduced chi-square test for all observations. As shown in Figure 3d, the reduced chi-squares were tested for four components with v 2 n ¼ 0:525, five components with v 2 n ¼ 0:535, and three components with v 2 n ¼ 0:522; hence the three components are slightly smaller than the others. Therefore, the three components were applied to extract the most useful information from the GPS position time series. In addition, the optimal Laplacian smoothing factor was set to 10e3, which regularizes the ill-posed problem and imposes significant penalties on SSE patches. Indeed, the proposed inversion model is efficient for the spatiotemporal evolution of the 2019 SSE at the Hikurangi subduction interface.

Cumulative slip distribution
As shown in Figure 4, the positions of the east components are changing with an exponential function. The horizontal displacements are nearly perpendicular to the Hikurangi trench, and the movement pattern is consistent with the significant thrust movement at the plate interface with a maximum value of about 4 cm (MAHI, Figure 4). These three components are obtained from PCAIM, which satisfactorily explains the surface deformation characteristics. The temporal functions and spatial weights of the three principal components are presented in Figure 3. The first principal component dominates ( Figure 3a) and explains 77.0% of the variance in the filtered data, illustrating the primary spatiotemporal feature of the 2019 SSE following an exponential decay mechanism. In particular, the time vectors indicate a complex pattern of variation with multiple accelerations and decelerations. Spatial vectors show that slow slips exhibit thrust-dominated movement (Figure 3a), concentrated in the Gisborne and Hawke's Bay areas, similar to previous east coast SSEs Koulali et al. 2017;Warren-Smith et al. 2019). The second and third principal components (Figure 3b and c) present only 11.7% and 11.3% of the variance in the reconstructed data, respectively. Regarding their spatial functions, the second and third components still show a few slow slip motions in Gisborne and Hawke's Bay, respectively (Figure 3b and c). Although some seasonal deformation and observational noises, primarily in the GPS position time series, are not properly corrected in the second and third principal components; however, these components are still necessary for the inversion. Figure 5a shows the cumulative slow slip distribution over about 3 months, which is dominated by thrust faulting. It is located at a depth of <15 km, with a peak slip of about 10.5 cm near the Gisborne region, similar to previous northern Hikurangi SSEs Wallace et al. 2018;Wallace 2020). The shallow slow slip indicates that the major slip is concentrated in offshore Gisborne, while the moderate fault motions are relatively concentrated offshores in central Hawke's Bay (Figure 5a). This is opposite to the 2016 east coast SSE, with a significant slip near Hawke's Bay and the second slip near the Gisborne region Jiang et al. 2018). The released seismic moment is calculated based on a uniform shear modulus of 30 GPa and a Poisson's ratio of 0.25, a total of about 4.83ⅹ10 19 NÁm, which is equivalent to an Mw $6.8 event.
Figures 4 and 5b show a satisfactory fit between observations and predictions in the preferred slow slip model. However, the MAHI and PAWA stations demonstrate high-frequency fluctuations, which may be caused by noise. In order to verify the model's robustness, a resolution test on the Hikurangi subduction interface is performed by simulating the displacements that would result from a test slip with a prescribed slip. Two slow slip patches near Gisborne and Hawke's Bay and along the shallow (about 6-20 km) Hikurangi subduction interface are only partially recovered in the model according to the resolution test ( Figure S2). The result indicates that the proposed inversion model can recover major slip regions of Hawke's Bay and Gisborne.

SSE evolution characteristics
In order to investigate the spatiotemporal evolution characteristics of the 2019 SSE, a period of three months is divided into eight durations, 10 days each, during the inversion. As shown in Figure 6, the spatiotemporal evolution model reveals a four-stage process. A brief acceleration of the fault slip is characterized as stage 1 (see Figure 6a and b) in the Gisborne region. The maximum fault slip and corresponding slip rate reach 0.011 m and 1.1 mm/day from 21 March to 1 April, respectively (Figure 6a). The fault slip is then continued along the subduction margin with the enlarged area and magnitude. Its slip rate significantly increases to 5.5 mm/day with a maximum fault slip of 0.055 m from 2 April to 11 April (Figure 6b). Afterward, the slip slows down and continues for approximately two weeks and begins to spread to Hawke's Bay in late April (Stage 2, Figure 6c and d). For the two periods (12 April-21 April and 22 April-1 May), the peak slip magnitude is reduced to 0.016 m and 0.014 m on the subduction interface, and the corresponding slip rate is reduced to 1.6 mm/day and 1.4 mm/ day (Figure 6c and d). In early May (2 May-11 May), a slow slip patch appeared in Hawke's Bay region (Stage 3, Figure 6e). The fault slip rate significantly increases to 2.4 mm/day with a maximum fault slip of 0.024 m. Afterward, there are no obvious variations (i.e., slip peaks and rates) in the last three periods from 12 May to 10 June (Stage 4, Figure 6f-h).
In addition, the 2019 SSE migrated southwestward from Gisborne to Hawke's Bay along the Hikurangi margin in the opposite direction, as compared to the 2011 and 2015 SSEs (Table S1; Wallace et al. 2012;Koulali et al. 2017). For the 2012 and 2014 SSEs, the Tolaga Bay patch starts to slip first, followed by a larger SSE at Gisborne, which are different from the 2019 SSE (Todd et al., 2016;Wallace et al. 2016). The above results indicate that the 2019 SSE has two prominent evolution characteristics: one is that the 2019 SSE stars in Gisborne and is accompanied by a small slip patch in Hawke's Bay; the other is that the 2019 SSE shows an accelerateddecelerated-accelerated-decelerated motion during interface slip when compared to previous SSEs along the northern Hikurangi margin. Those evolution characteristics can be verified by the displacements at some GPS stations CKID, CNST, KOKO, and MAHI (Figures 2 and 4).

Linkage between SSE and seismicity
The 2019 east coast SSE was accompanied by abundant seismicity (Figure 6). These contemporaneous earthquakes were recorded by the national GeoNet seismicity catalog (https://www.geonet.org.nz/data/types/ eq_catalogue). It indicates that seismicity during the 2019 SSE is concentrated in a NE-SW trending band (parallel to the slab strike). The depth range of these seismic events is from 12 km to 33 km, with the highest seismicity at a depth of 25 km to 30 km.
Similar to previous SSEs, such as 2004, 2010, and 2014 east coast SSEs (Delahaye et al. 2009;Jacobs, Savage, and Smith 2016;Todd and Schwartz 2016;Todd et al. 2018), seismic events during the 2019 SSE mostly occurred within the Hikurangi plate ( Figure 6). Compared to the seismicity before the 2019 SSE, the seismicity increased substantially and lasted for several days at the early stage of the slow slip (Figures 6b and 7). Some seismic events are close to the major slip regions (Figure 6b). Moreover, several seismic events occurred in some secondary slip regions (southern Hawke's Bay and northern Gisborne) over the period between 2 May to 11 May during the 2019 SSE (Figure 6e, f, and h). Compared to the Gisborne region, the seismicity near Hawke's Bay has generally located downdip the SSE slip patches (Figures 6e-h, and 7), which appeared to be a different seismic pattern. The Mw 5.1 Gisborne earthquake is the largest event recorded at the end of the SSE on 14 May 2019 (Figure 6f), occurring updip of the SSE. The Coulomb stress change (DCFS) caused by slow slip ( Figure S3) is much lower than the triggering threshold of 10 kPa (King, Stein, and Lin 1994;Harris 1998;Ma, Chan, and Stein 2005;Guo et al. 2019;Zhang et al. 2021) with the values up to À0.06 kPa and À0.03 kPa on the two nodal planes, respectively. Thus, the 2019 SSE did not induce the Mw 5.1 Gisborne event based on the Coulomb failure stress model (see details in S2).
To explore the relationship between the seismic moment during every stage and seismicity, we added a stage between the April 2 and April 11, Figure 7. The histogram represents the corresponding daily counts of seismic events during the 2019 SSE. Red histograms are the corresponding daily counts of seismic events prior to the 2019 SSE. The green curve is the maximum daily slip rate in different slow slip stages, and the red curve is the seismic moment accumulated every slip stage (time intervals are 10, 4, 6, 10, 10, 10, 10, 10, and 10 days). Notably, we added a stage between April 2 and April 11 from Figure 6, and the time intervals are 4 and 6 days, respectively. The black arrows on the horizontal axis indicate the duration from the beginning to the end of the slip patches. The black arrows show the timing of the Mw 5.1 Gisborne earthquake. and the time intervals are 4 and 6 days, respectively (Figure 7). Notably, the daily number of seismic events during the east coast SSE is temporally correlated with the seismic moment of SSE and the maximum daily slip rate at different slip stages (Figure 7). Accordingly, the seismic events occur concurrently with the onset of moment release during slow slip. Some previous SSEs, such as 2004, 2010, and 2014 Gisborne SSEs, produced an apparent increase in seismicity rate changes during the northern Hikurangi SSEs (Delahaye et al. 2009;Jacobs, Savage, and Smith 2016;Todd et al., 2016;Yarce et al. 2021). In contrast, some earthquake swarms preceded the SSE-induced displacements several days, and the stress loading caused by the SSEs was not a plausible triggering mechanism for these pre-SSE earthquake swarms (Nishikawa, Nishimura, and Okada 2021). However, other shallow SSEs in Costa Rica (Outerbridge et al. 2010) and Nankai (Suzuki et al. 2016) showed different seismicity patterns from the one described in this study. This shows that the relationships between the seismicity and the SSEs may vary with specific subduction zone conditions. Additionally, the seismic moment released only about 0.2% of the total SSE moment (about 6.19ⅹ10 16 NÁm) from the GeoNet catalog during this period; thus, the contamination of seismicity-induced surface deformation in the GPS position time series is negligible.  Figure 8a, b, and f) are mainly spread in the Gisborne region, with about 9 slow slip episodes during 2011-2019 (Table 1). Historical onshore seismicity is primarily confined to the subduction plate near Hawke's Bay ( Figure  9b). Moreover, the depth range of earthquakes in the original GeoNet catalog ( Figure 9a) differs significantly from that of the relocated earthquakes (Figure 9c), especially for offshore seismic events Yarce et al. 2019). The GeoNet catalog shows that a substantial uncertainty can be related to the simplified 1 D velocity model and the lack of seismometers in offshore seismicity (Yarce et al. 2021). Although the offshore seismic events from the GeoNet catalog are not well constrained, the onshore seismic events between Gisborne and Hawke's Bay can illustrate the primarily seismological features. Cross-sections reveal that most historical seismic events are distributed in the seismicity zone (SISZ) near the Gisborne region, while offshore events have migrated near the subducted seamounts (Figure 9a). On the other hand, long-duration (about 3-11 weeks) and less frequent SSEs (such as the 2015a: January to February SSEs) tend to extend over Hawke's Bay, with about 5 slow slip episodes during 2011-2019 (Table 1; Figure 8c and d).

Impacts of subduction interface heterogeneity
Significant differences in the spatial distribution of SSEs and seismicity between Gisborne and Hawke's Bay can be attributed to the highly heterogeneous stress distribution at the northern Hikurangi subduction interface, possibly a result of thick underthrust sediments and shallow interface topography induced by subduction of seamounts (Barker et al. 2009;Bell et al. 2010;Pedley et al. 2010;Sun, Saffer, and Ellis 2020;Wallace et al. 2012;Todd et al. 2018;Wallace et al., 2020;Warren-Smith et al. 2019). Generally, the northern section (Gisborne) of the Hikurangi margin is characterized by a thin sediment cover and seamounts above the sediments, which leads to a steep and narrow accretionary wedge (Wallace et al., 2020). In contrast, the central (Hawke's Bay) and southern segments are covered with thicker sediment with a wider accretionary wedge ).  Koulali et al. 2017). Light blue and light green shades indicate high-amplitude reflectivity zones (HRZs) and low-amplitude reflectivity zones (LRZs), respectively. The three brown shades with numbers S1, S2, and S3 show the subduction seamounts (Bell et al. 2010).
Similar to previous SSEs in 2011, 2014, and 2016, most of the slow slip occurred near the Gisborne region, downdip of the high-amplitude reflectivity zones (HRZs: HRZ_2 and HRZ_3 with light blue), and subduction seamounts (Figure 8f; S2 and S3 with brown shades). It can be associated with the fault instability caused by the enhanced megathrust normal stress on the seamount's leading edge (Sun, Saffer, and Ellis 2020). The fault instability further facilitates fluid migration to the subduction interface and increases pore fluid pressure in the Gisborne region (Barker et al. 2018;Ellis et al. 2015;Chesley et al. 2021). The results demonstrate that the rough subducting plate can reduce the effective stresses by increasing the pore fluid pressure in the Gisborne region, and the rough subducting plate plays a crucial role in stimulating shallow SSEs and stick-slip behaviors near the slab (Bell et al. 2010;. This characteristic does not exist in Hawke's Bay (Wallace et al., 2020). Our results reveal that SSEs and seismicity are strongly related to the heterogeneous stress distribution at the plate boundary. A robust analysis of seismic behavior along a subduction boundary will require more relocated seismic data.

Conclusions
This study assesses the spatiotemporal evolution of the 2019 SSE at the Hikurangi subduction interface. The findings demonstrated that slow slip started offshore Gisborne and then migrated into Hawke's Bay irregularly along the subduction interface. Maximum slow slip for a depth range of <15 km is up to about 10.5 cm. It releases a total moment of about 4.83ⅹ10 19 NÁm, corresponding to an Mw $6.8 event. The depths of the primary slip region beneath the Gisborne and the secondary slip region beneath Hawke's Bay are similar, but their peaks, durations, and rates are slightly different. Combining the previous SSEs of 2011-2016, short-duration and frequent SSEs are mainly distributed in the Gisborne region, whereas long-duration and less frequent SSEs tend to spread over Hawke's Bay. In addition, concurrent seismic events occur with the onset of moment release during slow slip. Much of the historical offshore and onshore seismicity shows marked along-strike variations, which coincide with the variations in the subducting topography, such as subduction seamounts. Hence, results further imply that the heterogeneity of the subduction interface impacts the spatiotemporal distribution of SSEs and seismicity. GeoForschungsZentrum GFZ, for providing the PSGRN/PSCMP software Generic Mapping Tools were used to generate graphical illustrations.

Disclosure statement
No potential conflict of interest was reported by the authors.

Data and resources
All data in this work are available. Continuous GPS position time series are downloaded from New Zealand's GeoNet (http://www.geonet.orz.na) and PositioNZ (http://apps.linz. govt.nz/positionz) in our work (Figure 1). These contemporaneous earthquakes are from the national GeoNet seismicity catalog (https://www.geonet.org.nz/data/types/eq_catalogue). The supplementary material provides two texts, four figures, and previous SSEs data for this article to support those results and discussions.

Funding
This work is funded by the Sichuan Science and Technology Program (No. 2020GZYZF0010) and the National Natural Science Foundation of China (Nos. 42171429 and 41931076).