Low-frequency earthquakes accompany deep slow slip beneath the North Island of New Zealand

,


Introduction
The term "slow earthquake" is commonly used to describe fault-slip events at rupture velocities below standard earthquake rupture velocities and encompasses a range of phenomena such as tectonic tremor (Obara, 2002;Rogers & Dragert, 2003), low-frequency earthquakes (LFEs) (Shelly et al., 2007), very-low-frequency earthquakes (Ito et al., 2007) and slow slip events (Rogers & Dragert, 2003;Radiguet et al., 2012).These events are interpreted to represent shear slip along a fault (Shelly et al., 2007), similar to classic earthquakes, but with longer durations and less radiated seismic energy.Shelly et al. (2007) demonstrated that tectonic tremor (hereafter referred to as tremor) is, at least partially, the composite signal of many LFEs superposed over one another in time, suggesting that tremor and LFEs are different manifestations of the same phenomenon.The spatio-temporal correlation of tremors and LFEs with slow slip events has been extensively reported, especially where dense seismic networks have been installed, namely in Mexico, Cascadia, and Japan (e.g.Kostoglodov et al., 2010;Shelly et al., 2006;Bostock et al., 2012;W. B. Frank et al., 2014).LFEs in particular are now considered a seismic indicator of slow slip and can be used as in-situ monitor of slip (W.B. Frank et al., 2015;W. B. Frank, 2016;W. B. Frank & Brodsky, 2019).Uncovering previously undetected slow slip events using LFEs provides a means to improve the spatio-temporal resolution of images of slow slip along a plate boundary (W.B. Frank et al., 2018).We focus here on the Hikurangi margin of New Zealand, where there have been many reports of slow slip (Wallace, 2020), but there does not yet exist a catalog of low-frequency earthquakes.
New Zealand is located at the plate boundary between the Pacific and Australian Plates (Figure 1).Beneath the North Island, the Pacific Plate is subducting below the Australian Plate along the Hikurangi Subduction Zone, at a convergence rate ranging from 60 mm/year at the northern Hikurangi trough to 20 mm/year in the south, based on elastic block modeling of GNSS data (Wallace et al., 2004).The overall Pacific-Australian Plate convergence rate through the New Zealand region is 39-45 mm/year DeMets et al. (1990) and Beavan et al. (2002).This southward decrease in relative plate motions is due primarily to rapid clockwise rotation of the Hikurangi Subduction Margin's forearc (Wallace et al., 2004), which also leads to back-arc extension in the Taupō Volcanic Zone.
The northern South Island region is the site of the transition from Hikurangi subduction to strike-slip dominated motion along the Marlborough Fault System and onto the Alpine Fault.Slow slip events in the Hikurangi Subduction Zone exhibit diverse durations, magnitudes, and recurrence intervals that vary spatially (Douglas et al., 2005;Wallace & Beavan, 2006, 2010;Wallace & Eberhart-Phillips, 2013;Koulali et al., 2017;Wallace, 2020).
Deep slow slip can be observed at the Hikurangi subduction zone at three main locations along the plate interface: beneath the Kapiti coast region, northwest of Wellington (25-50 km deep), beneath the Manawatu region (15-50 km deep) and beneath the Kaimanawa ranges in the central North Island (30-40 km deep) (Figure 1) (Wallace, 2020).
Kapiti and Manuwatu slow-slip episodes have durations of 1 to 2 years, recurrence times of approximately 5 years and typically involve a large amount of slip (≥ 30 cm) (Wallace & Beavan, 2010).Manawatu slow slip events have been observed geodetically with GNSS positioning in 2004-2005, 2010-2011, and 2014-2015.These slow slip events are observed to often occur soon after Kapiti events, suggesting a northward migration from the Kapiti to the Manawatu region (Wallace et al., 2014).In contrast, slow slip events beneath the Kaimanawa range have shorter duration, lasting for 2 to 3 months and generating approximately 2-5 cm on the plate interface and were most clearly observed in 2006 and 2008 (Wallace & Eberhart-Phillips, 2013).Tectonic tremors have been documented at multiple locations throughout New Zealand.
In the Hikurangi margin many studies have reported observed ambient tectonic tremor at shallow depths offshore Gisborne and along the East Coast of North Island (Kim et al., 2011;Todd & Schwartz, 2016;Todd et al., 2018;Romanet & Ide, 2019).Deeper tremor has also been observed beneath the Manuwatu region (Fry et al., 2011;Ide, 2012;Romanet & Ide, 2019).Recently Romanet and Ide (2019) observed tremor beneath the Marlborough Fault System, but it remains unclear whether this is due to slip on the faults of the Marlbrough Fault system, or the underlying Hikurangi subduction system. A. Wech et al. (2012) observed tremor on the deep extent of the central Alpine Fault, and Romanet and Ide (2019) observed tremor beneath Fiordland in southern South Island, likely associated with slip on the Puysegur/Fiordland subduction system.So far, the only observations of LFEs in New Zealand have been made on the Alpine Fault (Chamberlain et al., 2014;Baratin et al., 2018) where they have been used to infer quasi-continuous slip on the deep extent of the fault.
Matched-filtering or template matching techniques, in which the seismograms of a template event are correlated with continuous data to detect similar waveforms, have been widely used to study tremors and LFEs (e.g.Obara, 2002;Shelly et al., 2007;Brown et al., 2008;Ide, 2010;Bostock et al., 2012;W. B. Frank et al., 2013;Chamberlain et al., 2014;Baratin et al., 2018;Sáez et al., 2019;Romanet & Ide, 2019).Several methods successfully extracted LFEs from tremor waveforms in the past (Brown et al., 2008;W. Frank & Shapiro, 2014;Poiata et al., 2018).We develop here a workflow to construct, pick and locate LFE templates with high-precision and generate the first LFE catalog beneath the North Island of New Zealand making use of the tremor catalog published by Romanet and Ide (2019).

Low-frequency earthquake detection workflow
Low-frequency earthquakes (LFEs) have potential to provide a powerful in-situ monitor of where and when slow slip events occur, but a detailed catalog is necessary to fully exploit this relationship (W.B. Frank & Brodsky, 2019).As LFEs are tiny signals buried in tremor waveforms (Shelly et al., 2006;W. B. Frank et al., 2013), we designed an approach inspired by Shelly and Hardebeck (2010), Chamberlain et al. (2014), Beaucé et al. (2019) and Park et al. (2020) with the key difference that we used tremors directly as template waveforms (Figure 1) to extract LFEs in an iterative approach of template matching, clustering and stacking.
For the matched-filter search, we have used the efficient GPU-based Fast-Matched-Filter routine developed by Beaucé et al. (2018).The workflow presented in this work is summarized in Figure 2 and incorporates three iterations of a matched-filter search (steps 3, 9, 13).In the initial matched-filter search (step 3), we use a composite template catalog containing 257 automatic detections of LFEs made using the BackTrackBB code (Poiata et al., 2016) (step 1) and 323 events from the tremor catalog of Romanet and Ide (2019) beneath the Kaimanawa Ranges (step 2).For the automated detection using BackTrackBB, we used parameters similar to those adopted byPoiata et al. (2018) and a 1D velocity model sampled from the 3D model of Eberhart-Phillips et al. (2010) at the center of the tremor cluster (step 2).For the templates issued from the tremor catalog Romanet and Ide (2019), the S-wave arrival times at each stations are derived from the theoretical arrival time from the location given in the catalog.For the templates extracted with BackTrackBB of Poiata et al. (2016), the S-waves arrival are obtained from the maximum of correlation between each pair of stations.
We scanned the continuous seismograms from 12 GeoNet broadband stations sampled at 100 Hz recorded between 2008 to mid-2020 (Figure S1) (step 3) and bandpass filtered at 2-10 Hz, representing the frequency-band in which the tremors have been observed (Romanet & Ide, 2019).We used 10 s-long templates, starting 2 s before and ending 8 s after the S-wave to account for picking uncertainties and to include the coda of the S-wave in the template waveform.We used an initial detection threshold of 7 times the daily median absolute deviation (MAD) and rely on later steps in the workflow to remove false detections associated with this low threshold.This low initial threshold allows us to detect closely-spaced but not necessarily co-located LFEs as well as repeats of the same family (W.B. Frank & Abercrombie, 2018).This first iteration resulted in 35779 detections.
Following initial matched-filter detection, we applied a signal enhancement step.This involves constructing a waveform matrix for each template, station and component, composed of waveforms starting 90 s before and ending 90 s after each detection to provide an estimate of the overall signal-to-noise ratio of the detection (step 4 and Figure S3a).Below, we refer to all the detections made by a specific template as a family, defined by the template's waveforms and location.Families with fewer than 10 detections were discarded.
We next applied a deblurring filter (step 5) (Lim, 1990), otherwise known as a localmean filter or Wiener filter (Moreau et al., 2017;Beaucé et al., 2019), to each family's waveform matrix.For this filter, we defined a sliding window of size N ×M , where N is the number of detections and M is the number of samples in the time domain.The deblurring filter is implemented as follows: if x is the input signal then the output y is such that where m x and σ 2 x are the local estimate of the mean and variance, respectively, and σ 2 is the noise threshold estimated as the average of all the estimated local variances.Applying this filter enhances the portion of the waveforms in which the local variance, i.e.
the variance of all the waveforms contained in the sliding window, is low compare to the noise threshold.It smooths the waveforms where the local variance is higher than the noise threshold by averaging each waveforms in the sliding window (Figure S3).After several tests of the effects of varying M , the width of the sliding window, after visual inspection, we conclude that 100 samples, which corresponds to 1 s for the 100 Hz-sampled waveforms, seems to be the optimal size for improving the signal-to-noise ratio.
Previous iterative stacking and matched-filtering routines have employed a simple linear stack of events within a family to enhance the signal-to-noise ratio of subsequent generations of templates (e.g.Chamberlain et al., 2014).However, we found that when we attempted this that our second-generation templates were degraded due to high noise levels.Instead of stacking all waveforms, we therefore identified the most similar waveforms within a family and stacked only these subsets of events to obtain higher stacking gain.With this approach, we constructed correlation/dissimilarity matrices for all families and employed an unsupervised learning approach to identify the most similar waveforms.
To construct the similarity matrix for each station and component, we correlated the detections with each other in a 10 s window centered on each detection (step 6).We allowed a ±0.5 s shift in the alignment of waveforms for correlation in order to find the maximum correlation coefficient between any two waveforms.To obtain a unique similarity matrix CC for each family, the correlation coefficient between event u and v, CC u,v , is computed as the mean of the correlation between u and v at all stations and components available.This unique similarity matrix is then converted into a dissimilarity ma- This dissimilarity matrix provides an indication of the distance between events such that low values in D indicate events close in correlation-space to each other and high values indicate disparate events in correlation-space (Chamberlain et al., 2018).Note that here, we set negative correlations to zero, and hence the maximum of D is 1.
To then construct clusters of similar events, we used the Hierarchical Agglomerative Clustering (HAC) (Müllner, 2011) unsupervised learning algorithm (Figure 2).The HAC algorithm begins with a forest of clusters as each detection is a single cluster.At each iteration the two closest clusters are merged, forming a branch.The algorithm continues until all the detections are gathered into a unique cluster known as the root.This is a "bottom-up" approach starting with no events grouped, and is in contrast to the Hierarchical Divisive Clustering (Kaufman & Rousseeuw, 2009) which involves a "top-down" approach in which all events are initially linked in a single cluster.
To compute the distance between clusters we used the average linkage method (Sokal, 1958), which estimates the distance between clusters as the average of the pairwise distances between potential cluster members.As described by Park et al. (2020), we judged the average linkage method more suitable to our approach in comparison to other methods such as the single or the complete linkage methods that consider the closest and the farthest detections in each cluster, respectively, to merge them together.The average linkage between clusters A and B can be written as: where, D u,v is the distance between event u and v.
We use a dendogram to illustrate the arrangement of the observations (Figure S2).
To delineate clusters within each family, we must select a dissimilarity threshold above which events are treated as unclustered.We assume that families contain both low-quality (noisy) detections and high-quality (clear) detections, and that the high-quality detections are more similar to each other than to the poor-quality detections whose correlations may be dominated by uncorrelated noise.We can thus define a dissimilarity threshold to select the high quality waveforms for stacking to create a high signal-to-noise ratio waveform representative of the entire family of detections.
To define the dissimilarity threshold for each family we chose the highest dissimilarity that allowed us to regroup 80% of the family's detections into one cluster.By doing so, we are effectively excluding from the stack the 20% of the detections that are the least similar to the rest of the family (Figure 2 and S1).Removing approximately 20% of the least-similar detections is a compromise between the completeness of the catalogue and the quality of the extracted waveforms constituting the main cluster.

Location and relocation of the low-frequency earthquake candidates
Once extracted, we linearly stacked the deblurred waveforms of the main cluster for each station and component to create a new template of higher signal-to-noise ratio (step 7).We visually confirmed the higher quality of the new templates (step 8) before proceeding to a second iteration of the matched-filter search (step 9).With the second iteration we obtained 85856 detection and after repeating steps 4-7 with this new set of detections (step 10), we were able to manually pick S-waves arrivals on all families and P-waves, which were previously below the noise level, emerged for a majority of the LFE candidates (step 11).
After two iterations we were able to manually pick 445 P-and 963 S-phases on the stacked waveforms for the 111 remaining families at this stage (Figure 3a).We then computed absolute hypocenter locations for the picked events using the NonLinLoc algorithm (Lomax et al., 2000) and the same 1D velocity model adapted from Eberhart-Phillips et al. ( 2010) used for preliminary locations in step 1 with BackTrackBB (Poiata et al., 2016).Our selection of templates for the final matched-filter iteration relies on the quality of the location of these events, evaluated based on the maximum length of the 68% (or 1σ) error ellipsoid's three semi-axes.We retained the events that had all the 68% error ellipsoid semi-axes less than or equal to 20 km (Figure 3b and S2).This threshold represents a good compromise between the number of events that we would use for the next sections and not too large location uncertainties (Figure S4).A total of 108 over 111 families could be located with picks at 4 or more stations and 77 families were located with a 68% error ellipsoid with semi-axis smaller than 20 km.
The remaining families are concluded to contain true LFE detections for three main reasons: (1) the waveforms of the templates and detections are dominated by energy in the 2-8 Hz band as is characteristic for low-frequency earthquakes (Shelly et al., 2006); (2) the families' hypocenters are found to be located near the plate interface as expected (Brown et al., 2009); (3) their detections are dominated by bursts of activity or event swarms, the defining feature of LFE activity (W.B. Frank et al., 2014) (Figure 3c and Figure 5a).Swarms of regular earthquakes are also commonly associated with shallow SSEs at the Hikurangi subduction margin as well (Delahaye et al., 2009), although the depth and frequency content mark our detections as distinct from regular earthquakes.
After locating the LFE families with NonLinLoc (Lomax et al., 2000) and after selecting a subset of families for further analysis based on the size of their 68% error ellipsoid, we find that most events locate around the location of the plate interface from the interface model of Williams et al. (2013) (Figure 4).We expect LFEs to occur on the plate interface (Brown et al., 2009) and we thus compute relative locations of the LFE families to discern whether their locations collapse to the plate interface.To do so, we use the relocation algorithm GrowClust (Trugman & Shearer, 2017), a relocation algorithm that employs hierarchical clustering to find neighbouring events, and relocation within clusters based on the minimization of differential travel-time residuals within clusters.GrowClust makes use of differential travel time observations, cross-correlation values, and reference starting locations to group and relocate events.
For a given event pair u and v, the GrowClust algorithm computes a similarity coefficient Z ij of each distinct event pair.This similarity coefficient is the sum over the crosscorrelation values r ij;k for the k common stations within a maximum station distance ∆ max and that exceed a minimum value r min : To relocate the LFE candidates, we set the maximum inter-station distance ∆ max to 102 km, to be higher than the maximum of the inter-station distance of the our network which is approximately 101 km.We tested different values for r min ranging from 0.1 to 0.6, implying a low to high acceptance threshold of the similarity between event pairs.A higher threshold would lead GrowClust to create more localized clusters with higher accuracy.We noticed that increasing the value of r min led to a smaller number of relocated events but ultimately produced similar results for those events relocated.
We found that for all values of r m in, the LFE families were always relocated closer to the interface (see Figures 4, S3 and S4 for r min = 0.6, 0.4 and 0.2 respectively).With r min = 0.6, 44 LFE stacks were relocated.

Source Mechanism Estimation
In addition to LFEs having locations consistent with the plate interface, we expect LFEs to have source mechanisms consistent with shear failure on planes with similar geometry to the regional subduction interface.Determining the source mechanism of lowfrequency earthquakes has been challenging in the past due to the events' characteristically low signal-to-noise ratios, but past observations of shear mechanisms have provided strong evidence of LFE relationships to plate boundary slip (Shelly et al., 2007;Ide et al., 2007;W. B. Frank et al., 2013;Baratin et al., 2018).Because our events represent the first LFE detections on the deep extent of the Hikurangi Subduction Zone, we attempted to determine their geometric consistency with slip on the plate interface, and rule-out their possible association with deep volcanic events (e.g.Reyners, 2010;Hurst et al., 2016).
We first attempted to fit the waveforms of the LFE stacks with synthetic ones for a range of different assumed geometries, following a similar approach to that of W. B. Frank et al. (2013).We generated synthetic waveforms using the Axitra code (Coutant, 1989) and Green's functions derived from the 1D velocity model sampled from the 3D model of Eberhart-Phillips et al. (2010).We aligned the P and S arrivals of the LFE stacks and synthetic waveforms using cross-correlation, and we normalized the stacked LFE waveforms by the amplitude at the station where it is maximum for the S-wave and normalized the synthetic waveforms so the amplitude at the same station is 1.Finally, we computed the average root-mean-square amplitude difference between real and synthetic waveforms over all stations and components.We were unable to obtain a compelling result using this method and could not find a common mechanism for the LFE candidates.It is likely that our workflow, combining the deblurring filter and stacking (Section 2) may have altered the waveforms in a non-linear sense relative to the true LFE source mechanism.
As an alternative approach, we opted to investigate the amplitude ratio between P and S waves at the different stations of the network used in this study.Measured S/P manuscript submitted to JGR: Solid Earth ratios reflect the radiation patterns of P-and S-waves and are thus indicative of an earthquake's focal mechanism (Hardebeck & Shearer, 2003).
To simplify the determination of LFE source mechanisms, we assume a double couple LFE source and search over the strike, dip and rake of focal planes.We generated a data set of synthetic waveforms with a sampling rate of 20 Hz for a source located at the barycenter of the relocated LFE candidates (approximately at a latitude of -39°, longitude of 176°E and depth of 50 km) with a strike ranging from 0 to 360°sampled every 20°, a dip between 0 and 90°, sampled every 10°and we used a rake fixed at 120°.This assumption is based on observations during deep slow slip events beneath Manawatu (Wallace & Beavan, 2010) and the Kaimanawa ranges (Wallace & Eberhart-Phillips, 2013) where the direction of slip on the interface appears to be oblique (component of right lateral and reverse) and parallel to the Pacific-Australia plates motion (Wallace et al., 2004;Wallace & Beavan, 2010, and references therein).
We then compared these synthetic waveforms to the LFE stacks and searched for the best strike and dip angles that describe the distribution of the measured amplitude ratios.Before computing the amplitude ratios of both the real and synthetic waveforms, we down-sampled the real LFE candidates waveforms to 20 Hz and bandpass filtered both the real and synthetic data to 2-4 Hz to simplify the waveforms.For each candidate and at each station, we measured the maximum amplitude of the P wave on the vertical component and divided this by the mean of the S-wave amplitude maxima measured on both horizontal components for both real and synthetic data.For individual stations we compared the ratio measured on the synthetic waveforms for each source mechanism to the distribution of the amplitude ratios for the real data, represented by the violin-plots in Figure 6a.
To then determine the optimal mechanism, we computed a score representing the proportion of stations for which the synthetic amplitude ratio falls into the 10-90% interquartile range of the distribution of the observed amplitude ratio.This allowed us to map the scores obtained for different strikes and dips (Figure 6b).We found that the mechanisms that have the highest score (0.83) for all LFE families have a strike of 240°and a dip of 30°(Figure 6a and b) for a rake of 120°for a source located at the barycenter of the LFE stacks.This source mechanism is consistent with the expected geometry of the Hikurangi subduction (e.g.Wallace et al., 2004;Williams et al., 2013).

Discussion and Conclusions
In the work presented here, we developed an original methodology to extract lowfrequency impulsive signals buried in tremor that we interpret as low-frequency earthquakes (LFEs).Our workflow combines matched-filtering, clustering and stacking in an iterative approach to increase the signal-to-noise ratio sufficiently to manually pick P and S waves arrivals.After only two iterations, we were able to build a catalog containing more than 300 times the initial number of events in the tremor catalog of Romanet and Ide (2019).

Evidence for low-frequency earthquakes
To demonstrate that these events are indeed low-frequency earthquakes, we investigated different spatio-temporal characteristics of the families and estimated their source mechanism, assuming that they involve reverse faulting.As seen in Figures 3 and 5, the detected events tend to cluster in time in burst-like episodes; this is also true if we consider the activity of single families.Such behavior has been observed for low-frequency earthquakes in other regions and concluded to be associated with slow slip activity (e.g.Shelly et al., 2006;W. B. Frank et al., 2014;Chamberlain et al., 2014;W. B. Frank, 2016).
Further evidence supporting the tectonic nature of these low-frequency earthquakes is the location of the families.We first located our events using NonLinLoc (Lomax et al., 2000) using manually picked P-and S-wave arrivals for 71 families.Our subsequent relocation of these events using GrowClust (Trugman & Shearer, 2017) provides higherresolution locations close to the plate interface, at approximately 50 km depth (at -39°N latitude and 176°E longitude).However, our location procedure only made use of a 1D velocity model sampled from the 3D model of (Eberhart-Phillips et al., 2010).In the context of a subduction zone, and the nearby Taupō Volcanic Zone, this approximation may have a pronounced effect on the travel times computed by NonLinLoc, potentially introducing errors in the relocation.However, given the methodology presented here which stacks filtered waveforms which can lead to hardly quantifiable changes in the observed arrival times, we judged that the use of 3-D model would not bring significant information.
The final evidence we present in this work is the likely source mechanism of the detected events.By comparing the amplitude ratio of P-and S-waves between the LFE candidates and synthetic waveforms, we were able to identify the strike and the dip of the most representative double-couple source mechanism of the candidates.We additionally tested two different depths, 45 and 55 km (see Figure S9) which resulted in similar results with strike angles ranging from 220 to 240°and dip angles ranging from 30 to 50°knowing the direction of the dipping slab and assuming the convention that the interface is dipping to the right.The resulting source parameters are in agreement with the parameters expected for a rupture on the plate interface at these depths (Wallace et al., 2004).
The elastic-block model of Wallace et al. (2004) shows that all the strike-slip component has been accommodated by the rotation of the East part of the North Island and crustal faults.This implies that at depth, rupture on the interface should be oriented West-East parallel to the direction of convergence between the Australian and Pacific plates (DeMets et al., 1990).Manawatu's SSEs also exhibit a generally east-west direction of slip along the interface (e.g.Wallace & Beavan, 2010).
Taking these three key observations into account, we are confident that these events are LFEs associated with deep slip along the the Hikurangi subduction plate boundary.
These LFEs thus represent a unique opportunity to investigate the slip history of the deeper portion of the seismogenic zone beneath the North Island.

Workflow Enhancements
Although our methodology has successfully detected LFEs for the first time on the deep extent of the Hikurangi margin, several areas can be further examined and potentially improved.The use of the hierarchical agglomerative clustering algorithm on the dissimilarity matrix (built from the correlations between all the detections in a family) is less conservative than using the correlation between the detections and the parent template.The comparison of all detections with all other detections potentially provides a means to extract a new generation of templates that represent a somewhat different source to the original template, but is more representative of the entire family of detections.We also observed that for a few families the hierarchical agglomerative clustering was generating secondary smaller, but significant, clusters that could be used to construct templates for "sub-families".These secondary clusters could reflect events at other locations and/or with a different source mechanism to the original templates.
The deblurring of waveforms before the linear stack significantly improves the signalto-noise ratio.The number of the samples chosen for the deblurring moving-window is important as a short window length would be unlikely to improve the signal-to-noise ratio (e.g. with 25 samples in Figure S10) while longer windows will actually distort and smooth the waveforms (e.g.1000 samples in Figure S11).Moreau et al. (2017) suggested a rule of thumb where the number of samples in the deblurring moving window should be set to M ≥ F s /F m , where F m is the highest frequency in the signal and F s the sampling rate.By choosing M = 100 (≥ F s /F m = 100Hz/10Hz), we smooth the signal by choosing a longer time window.
After the deblurring filter and the clustering steps we chose to use a linear rather than non-linear stacking technique, such as a phase-weighted stack (Schimmel & Paulssen, 1997;Thurber et al., 2014).We did not use non-linear stacking because, although these methods can greatly improve the signal-to-noise ratio, they also distort the waveforms, as pointed out by Baratin et al. (2018) and Beaucé et al. (2019).This kind of stacking could nevertheless be useful for phase-picking (Thurber et al., 2014;Baratin et al., 2018).
The deblurring filter we employed is not the only possible signal-enhancement method available, and we also tested the Singular Value Decomposition based Wiener filter (SVDWF) proposed by Moreau et al. (2017) and later used by Beaucé et al. (2019) in a deblurring and stacking routine as well.The SVDWF includes spectral filtering, which keeps a certain number of singular vectors extracted from the singular value decomposition, and a deblurring filter also called a Wiener filter.We found that applying the deblurring filter on a limited number of singular vectors did not improve significantly the signal to noise ratio in comparison of using only the deblurring filter on the largest cluster of waveforms.

LFE Occurrence
Regarding the weekly count of detection shown in Figure 5a, we observe a strong burst of activity in 2008, starting at the beginning of the catalog.We also observe that the inter-event time (Figure 5b) is reduced and this period could be interpreted as a very large and continuous burst of LFEs.However, we note that the stations we chose to use for this study were not all fully operational until 2015 (Figure S1), especially for this particular period in 2008, when only 3 stations were available.This implies that with fewer stations, we have more detections but probably with a larger false detection rate.Nevertheless, the peak in detections seen after 2015, similar to the one in 2018, does not correspond to a significant change in the number of stations.In addition, when we used a matched-filter threshold of 10 x MAD, resulting in fewer, but higher-quality detections, this initial burst of activity remained (Figure S6).We suggest that the bursts of LFEs observed after 2008 are real and represent episodes of increased slip-rate on the deep extent of the Hikurangi subduction margin.
While Hardebeck and Shearer (2003) highlights the usefulness of the S/P amplitude ratio in the estimation of source mechanism parameters, they sometimes observed a scattering of the S/P ratio at the same station for groups of events.Location and velocity model errors can partly explain this scatter as well as potential differences in source mechanism.However, they found that the noise level had a large effect on the scattering of the S/P ratio.Here, we computed the amplitude ratio at a single location i.e. the barycenter of the LFE families for a first order analysis.This obviously led to a scattering of the ratio and potentially large distributions (Figure 6a) for the LFE stack measurements at each station, thus justifying a rather large interquartile range (10%-90%) when comparing the observed to synthetic amplitude ratio to estimate the fitting score.
We speculate that after more iterations of our approach, the signal-to-noise ratio would increase, potentially reducing the scatter of the amplitude ratio.This approach could potentially be used with several barycenters, accounting for smaller clusters of events rather than considering all LFE families at a single location.
A detailed spatio-temporal analysis of the LFE catalog with respect to the continuous GNSS positioning is an important next step.As noted by Romanet and Ide (2019) and as we observe in Figure 5, only one known slow slip event coincides with increased LFE activity in 2010.We can see that this increase of activity into two bursts that correspond to the beginning and the end of the SSE in the GNSS time series.This SSE was located beneath the Kaimanawa ranges and Manawatu region (Wallace & Beavan, 2010) and is located updip of the cluster of LFEs but does not seem to overlap with it (Fig- gions are analogous to the frequent tremor episodes observed below the geodetically detectable (and less frequent) episodic tremor and slip events (e.g.Obara et al., 2010;A. G. Wech & Creager, 2011).That said, we have selected the families based on their location accuracy obtained from NonLinLoc (Lomax et al., 2000), and these well-located families may not represent the complete spatial extent of LFE occurrence here.To identify welllocated events, we used the three axes of the 68% error ellipsoid derived from the 3D location PDF given by NonLinLoc and in particular half of their total length (Figure S5).
Setting a threshold on one or several of the semi-axis lengths allowed us to filter out events that present a poor location accuracy, i.e. large semi-axis.Visual inspections of the distribution of the dimensions of these ellipsoids for each LFE stack (Figure S12) showed that imposing a threshold on only the first and second axes would not filter many of the events; several ellipsoids would still exhibit large third semi-axis (≥ 40 km).However, applying a threshold of 20 km only to the third axis allows us to filter out LFE stacks that are poorly located.We noted that the families presenting both a characteristic burstlike behavior and a location close to the plate interface have a third semi-axis length below 20 km.We then relocated the LFE families with GrowClust (Trugman & Shearer, 2017).In comparison, horizontal and vertical uncertainties obtained from GrowClust are on the order of 3-4km (Figure S13).We noticed that increasing the correlation threshold r min from 0.1 to 0.6 did not significantly affect the location uncertainties while the new locations can change drastically (see Figures 4, S7 and S8), as well as an expected decrease in the number of relocated events.This means that the quality of the relocation is independent of r m in, hence the differences in locations obtained with different threshold may come from the 1D velocity model sampled from the 3D model of Eberhart-Phillips et al. (2010) and the initial locations given to GrowClust.It is difficult to discuss location accuracy in more detail as we were unable to propagate the uncertainties from NonLinLoc into GrowClust with the aim to consider the events location as a probability distribution and not as a single point in space.Although the LFEs appear to be located largely down-dip of the geodetically detectable SSEs, given the uncertainties presented above about the LFE stacks and the uncertainties and the location of the slowslip events, we can't be certain that there is no overlap between the LFEs and known SSEs.Other known SSEs in 2008 and 2014/2015 (Wallace & Beavan, 2010;Wallace, 2020) do not correspond to significant bursts of LFEs activity.If LFEs detected here are driven by deep slow slip, then geodetic observations would potentially have difficulty capturing the surface signature of such deep slow slip, particularly if such events are small and relatively frequent, as the LFE bursts suggests.The lack of significant LFE bursts during geodetically-detect SSEs also suggests that drawing a direct correlation between slow slip and LFEs (and using the LFEs as a way to monitor slow slip) may be less straightforward in the central Hikurangi margin, compared to other slow slip regions ((W. B. Frank et al., 2015;W. B. Frank, 2016;W. B. Frank & Brodsky, 2019)).This may be due in part to the spatial separation between the LFE region and the geodetically-detectable SSE source region.
We present in this work an original approach to extract low frequency earthquakes from the noisy signal of tremors.We applied this approach to the tremor activity occurring beneath the Kaimanawa range of the North Island, New Zealand to build the first catalog of low-frequency earthquakes in the Hikurangi margin.Future work investigating precisely where and when LFEs are occurring with respect to slow slip will help to improve our understanding of the potential interplay between the aseismic and seismic component of the earthquake cycle at depth in the Hikurangi subduction zone.
manuscript submitted to JGR: Solid Earth  -19-manuscript submitted to JGR: Solid Earth see Figure 4) and the recurrence time between consecutive events against the cumulative number of events (for all families) along time.The shaded areas correspond to known slow slip events and their color correspond to their location (see Figure 1).
-21-manuscript submitted to JGR: Solid Earth Figure S4) and the recurrence time between consecutive events against the cumulative number of events (for all families) along time.The shaded area correspond to known slow slip Events and their color correspond to their location (see Figure S4).

Figure 1 .
Figure 1.Tectonic setting of the study.The tremors (Romanet & Ide, 2019) used as templates in this study are represented by red circles.The low-frequency earthquakes (LFEs) detected in this study are represented by the purple circles.The GeoNet stations used are the inverted yellow triangles and the contour lines mark the accumulated slip which occurred during the most recent Slow Slip Events (Wallace, 2020, and reference therein).The black line marks the plate boundary between the Pacific and the Australian plates while the dashed lines shown the position of the subducting slab at depth following Williams et al. (2013) plate interface model.

Figure 2 .Figure 3 .
Figure 2. Schematic view of the worklow constructed to extract and refine LFE waveforms based on the Hikurangi tremor catalog developed by (Romanet & Ide, 2019).

Figure 4 .
Figure 4. Relocation of the LFE candidates using GrowClust with a rmin = 0.6.The original location obtain with NonLinLoc is represented by the empty orange circles.The black lines show the distance between the initial and final locations obtained with GrowClust and represented by purple circles.The bottom plot shows a East-West projection at depth. the dashed lines shown the position of the subducting slab at depth following Williams et al. (2013) plate interface model.The two squares represent the GPS stations THAP (green) and VGMO (blue) used for comparison in this study (Figure 5)

Figure 6 .: 6 :Figure S3 .:Figure S4 .:
Figure 6.Source mechanism strike and dip estimation based on S/P amplitude ratios.(a) Comparison of the S/P amplitude ratios at each station between the final low-frequency earthquakes templates (violins) and the synthetic waveforms assuming a reference source mechanism with a strike of 240°, a dip of 30°, and a rake of 120°.The score represents the percentage of synthetic amplitude ratios that fall into the 10%-90% interquartile range (IQR) of the observed distribution.(b) Representation of the score for different strikes and dips for a given rake of 90°and a depth of 50 km.The green-outlined beachballs represent the mechanisms with the highest score (here 0.83).The background levels of gray show the distribution of the score with light colored areas representing higher scores.