The current period of volcanic activity has been ongoing for the last 2 ka., described as periodic low-volume ( < 0.05 km³) phreatomagmatic eruptions which occur every 25–30 years (Kilgour et al., 2013; Conway et al., 2016). Since the first historically recorded eruption in 1830, there have been 35 eruptions with a recorded Volcanic Explosivity Index (VEI) > 1 (Hurst et al., 2018). Major eruptions have a repose period of around 50 years with the last three occurring in 1895, 1945 and finally 1995/1996. Notably, the 1945 major eruption resulted in the creation of a new lava dome and debris barrier confining the crater lake (Johnston et al., 2000). Crater lake breakouts at Mt. Ruapehu are relatively common and have occurred repeatedly throughout the 20th and 21st century (Houghton et al., 1987; Schaefer et al., 2018). In 1953, the sudden collapse of the crater debris wall created a crater lake breakout and the formation of a lahar which destroyed the Tangiwai railway bridge causing the death of 151 people (Houghton et al., 1987; Johnston et al., 2000). The 1995/1996 eruption sequence was the largest in the 50 years prior (Sherburn et al., 1999). The sequence began with a localised phreatomagmatic eruption on 17 September 1995 (Hurst et al., 2018) and developed into a major eruption undergoing 10 eruptive phases, before ending in December 1996 (Bryan and Sherburn, 1999). The most recent recorded eruptions were two steam-driven eruptions occurring on 4 October 2006 and 25 September 2007 (Jolly et al., 2010; Keats et al., 2011; Carniel et al., 2013). The sub-aqueous 4 October 2006 eruption may not have produced an eruption column, but did significantly raise the crater lake water level (Jolly et al., 2010), which has been linked to the 18 March 2007 tephra dam collapse and subsequent lahar on 17 March 2007 (Carrivick et al., 2008). Another, larger phreatic eruption occurred on 25 September 2007 (Jolly et al., 2010). Although there were insufficient precursory data suggesting an eruption was imminent at the time, minor volcano-tectonic earthquakes and tremor were found in the seismic data to have preceded the eruption by around 10 min (Jolly et al., 2010). The main eruption lasted no longer than 1 minute and produced an eruption column which ejected ballistics 2.5 km to the north, with lahars entering two local catchments and disrupting the Whakapapa Ski field (Kilgour et al., 2010).

Mt. Ruapehu hosts its own tectonic setting, with faults circumscribing almost the entire base of the volcano some 20 km from the summit (Leonard et al., 2021). Shallow seismicity ( < 20 km depth) over the last 30 years occurs in distinctive high-density clusters at depths along the Raurimu Fault, which runs north to south along the west flank of Mt. Ruapehu (Figure 2). Analysis of earthquake epicentre distributions from the GeoNet earthquake catalogue, suggests there are two distinctive regions, to the northwest and southeast of the Mt. Ruapehu summit, where earthquakes tend to occur more frequently than others in characteristic clusters (Figure 2). Seismicity originating from deeper sources reflects the slip geometry of the Pacific Plate subducting beneath the Australian Plate (e.g., Yarce et al., 2019) from the southwest at depths 50–100 km, to the northeast at depths > 150 km.

To the southeast, approximately 20 km from the summit, is a band of near-constant seismicity along the Rangipo (Desert Road) Fault, near the town of Waiouru (Hayes, 2004) (Figure 2). Average b-values have been calculated at 1.06 at mid-crust depths ranging from 11 to 21 km implying a tectonic origin (Hayes, 2004). The b-value describes the relationship between large and small earthquakes in a catalogue, commonly used to describe the nature of the seismic activity of a region (Gudmundsson, 2020). Increased seismic activity was detected prior to the 1995 major eruption and once again in 2000/2001, with b-values sustaining high levels through the 1995/1996 eruption sequence reaching 1.71 in 2000 (Hayes, 2004). Although the earthquake cluster seems connected to volcanic episodes at Mt. Ruapehu, it is unlikely to be of volcanic origin itself and instead may be associated with regional stress changes following the 1995 major eruption (Hayes, 2004).

The seismicity to the west of Mt. Ruapehu is mainly concentrated northwest of the summit at the north extent of the Raurimu fault near National Park. There is also another distinct band of seismicity, which can be observed near Erua (Figure 2). Earthquake locations then become more uniformly distributed to the west and southwest towards Horopito and the Tohunga Junction. The cluster of earthquakes near Erua has been defined as precursory signal to the 1995/1996 major eruption sequence (Hurst and McGinty, 1999; Hurst et al., 2018). Studies by White and McCausland (2016) and Meyer et al. (2021) have developed methods, originating from McGarr (1976), for estimating the intrusive volume of magma intrusions from tectonic earthquake sequences at distal faults to suggest they might be suitable mid-term eruption forecasts. Hurst et al. (2018) argues the petrological timescale of erupted materials coincides with the seismicity near Erua, before the 1995 eruptions, suggesting that the seismic signal was an indicator of magma movement. Other studies have also suggested that the 1995 Erua swarm was not of tectonic origin, but a signal of stress changes associated with fluid movement within the volcanic system (e.g., Hayes, 2004). This may not be limited temporally to the 1995/1996 eruption, and may instead be an ongoing signal typical of pore-fluid pressure triggered by fluid movement in the mid-crust (e.g., Keats et al., 2011). However, a lack of observed seismic rate increase prior to the 2006 and 2007 eruptions makes eruption forecasts using this seismic signal uncertain. For the swarm near Erua to be considered a reliable forecasting signal (e.g., Hurst et al., 2018), it must be unique to the months preceding volcanic eruptions at Mt. Ruapehu, and should be absent during periods of volcanic inactivity (e.g., Ardid et al., 2022).

The earthquake catalogue for the Mt. Ruapehu region was downloaded for a period from 1990 to 2023, containing 28,522 earthquakes at shallow depths ( < 20 km) containing information on earthquake event type, time, location, depth and magnitude (Figure 2). Information including hypocentre location error is included in the Supplementary Material. Event types not relating to naturally occurring seismic events (e.g., quarry blasts and explosions) were removed, as were earthquakes located at the Mt. Ngauruhoe and Tongariro summits, which neighbours Mt. Ruapehu to the northeast, along with the near-persistent earthquake hot spot the southeast near Waiouru (e.g., Hayes, 2004) (Figure 2). The final study area was refined to a rectangular area to the west of Ruapehu [−39.09, 175.27, −39.43, 175.58] spanning approximately 1,800 km² to the west of the summit, to focus on activity exclusively to the west of Mt. Ruapehu, where the precursory signals at distal faults were identified by Hurst et al. (2018) (Figure 2).

Approximately half of the earthquakes in the catalogue (14,587) were assigned depths of 0, 5, or 12 km by an internal operator. These approximations for the hypocentre of the earthquake may be necessary when there are fewer seismic stations available to make a measurement, or when there is uncertainty in the velocity model used. Keeping earthquakes with locations biased to the assigned depth would be detrimental to the unsupervised machine learning approach we adopt and were therefore removed from the catalogue. Although this is a large proportion of the earthquakes in the catalogue, we argue that these earthquakes generally fit the underlying geographical trend of the seismicity of the region, and do not add any significant information to the overall catalogue (Figure 2). Therefore, the operator assigned events were removed from the catalogue for the purpose of this research.

The final data processing step was to assess the magnitude of completeness (M _c ) of the catalogue with operator assigned earthquakes removed. This was done by estimating the point of maximum curvature for the cumulative and non-cumulative frequency magnitude distribution (FMD) (Pavlenko and Zavyalov, 2022). Using the Maximum Curvature (MAXC) method (Wiemer, 2000), bin sizes of ΔM _L 0.1 were evaluated to find the best fit M _c value, which was calculated as 1.6 with a b-value of 1.04 ± 0.02. Finally, the earthquake catalogue was filtered to contain only earthquake events ≥M _c , totalling 2,795 individual events.

Density-based clustering algorithms are effective at analysing noisy data sets and can handle irregular cluster shapes, making them a good method for analysing the clustering and zoning of earthquakes (Scitovski, 2018). Clustering algorithms, such as DBSCAN (density-based clustering algorithms with noise) (Ester et al., 1996) have been effectively used to identify and classify seismic swarms in earthquake catalogues (e.g., Petersen and Pankow, 2023). Where the methodology described in Petersen and Pankow (2023) begins with identifying temporal clusters in the catalogue before spatial analysis, this study aims to identify spatial clusters prior to the temporal cluster analysis using density-based clustering algorithms. The seismic activity in the months prior to the 1995 eruption was confined to a distinctive cluster Hurst and McGinty (1999), therefore we attempt to observe whether this kind of seismic signal is unique in specific regions (e.g., Hurst et al., 2018) prior to volcanic eruptions or whether it repeats during periods of volcanic quiescence.

Hierarchical density-based clustering algorithms with noise (HDBSCAN) (Campello et al., 2013) is an unsupervised machine learning clustering algorithm that is particularly useful at identifying clusters of varying shape and densities within a noisy data set. HDBSCAN begins by computing the density of data points based on their proximity to neighbouring points using a similar methodology as the popular density-based algorithm DBSCAN (Ester et al., 1996). DBSCAN begins by selecting a random point in the data set (x) and, using Euclidean distance, calculates whether there are enough points in its neighbourhood to begin building a cluster, this is labelled a core point. Once a core point has been satisfied, the cluster is expanded to detect all core points which are directly-reachable in the latitude, longitude and depth dimensions. Once all of the core points for a cluster have been assigned, the algorithm adds the non-core border points to the cluster. The border points cannot be used to connect any other points to extend the cluster, even if another non-core point is within the distance threshold of the border point. Data points which fall outside the threshold defined by the input parameters are labelled as outliers and given the value −1. These steps are processed sequentially, meaning a non-core point which is within the distance threshold (ɛ) of two or more separate clusters will be assigned to the cluster which was calculated first (Figure 3).

HDBSCAN builds upon the DBSCAN algorithm by allowing the computation of clusters with varying ɛ distance, which is effective for grouping clusters with a range of densities. This is achieved by using the concept of mutual reachability distance (MRD), where the influence of multiple points are considered in relation to their local neighbourhood density, defined as: d m r e a c h − k a , b = max c o r e k a , c o r e k b , d a , b Where d(a, b) is the original distance metric between points a and b. HDBSCAN then constructs a minimum spanning tree (MST) from the MRD values, which is designed to highlight the strongest connections between points. The hierarchical part of the algorithm begins by using a single-linkage approach to repeatedly merge clusters close in density, ensuring that the merged cluster is associated with the cluster having the highest minimum density of the two. As more clusters are merged, the algorithm constructs a condensed tree that captures the hierarchy of cluster merges based on density. Each level of the tree represents a different level of minimum density required for a cluster to exist. This allows for a range of cluster sizes and shapes to be detected. The stability of a cluster is measured by how often it appears across different density levels in the condensed tree. Clusters that are more stable (i.e., they persist across multiple levels) are considered more significant. The final clustering result is obtained by selecting the clusters with the highest stability, while noise points and outliers are also accounted for. This allows HDBSCAN to automatically determine the total optimal number of clusters and handle clusters of varying densities.

Where HDBSCAN inherently selects the optimal number of clusters and corresponding distance between points, the algorithm requires a domain knowledge for the minimum number of samples per cluster, i.e., the minimum number of earthquake events for a cluster, (min_cluster_size). The min_cluster_size parameter determines the size of clusters calculated in the data set and impacts the total number of noise values (Starczewski et al., 2020) and has a default value of 5 data points per cluster. Although there are no constraints as to the number of earthquake events that make up a swarm, we can assume the earthquake signals we are attempting to identify contain greater than 5 events per cluster. In order to select the optimal number of points per cluster, one can run a range of min_cluster_size values through the clustering algorithm and test each result against a fitness function in a method known as Tournament Selection (TS) (e.g., Karami and Johansson, 2014). The best-fit min_cluster_size candidate based on the score of the fitness function is then selected for the clustering algorithm. Although there is no silver bullet fitness function for evaluating clustering algorithms, the Silhouette Coefficient θ (Rousseeuw, 1987) is a simple analytical method that can be used to find the optimal number of clusters for a variety of different clustering algorithms (e.g., Wang et al., 2017; Angmo et al., 2021) and does not require any prior training to produce results (Shutaywi and Kachouie, 2021). The Silhouette Coefficient θ is available in Scikit-learn and begins by comparing the distances of each sample to its own cluster and the neighbouring cluster, with scores ranging from −1 to 1: θ i = b − a m a x a , b Where a is the average distance of sample point i within a to other point within cluster a and b is the distance of sample point i to its nearest cluster. The function then calculates the mean Silhouette Coefficient over all cluster to return a single score for all clusters: θ = 1 n ∑ i = 1 n θ i where n is the number of samples and θ(i) is the Silhouette Coefficient for the i-th sample. A score of 1 would be achieved for two perfectly distinguished clusters at infinite distance from each other. A score of zero means that clusters are insignificant and the distance between them is not significant. The Silhouette Coefficient evaluated the goodness of fit for min_cluster_size 5–40 where the optimal number of samples in a cluster was 31 with a score of 0.19. We found that the cluster size did not change much for min_cluster_size 30–40. However, the signal near Erua was most prominent for min_cluster_size 30–32 and the Silhouette Coefficient reduces as sample size > 33. Therefore, we selected a min_cluster_size of 31 for the HDBSCAN algorithm. We noticed an increase in volatility for the total number of clusters, given a sample size < 30 samples per cluster, where clusters tend be more unstable and form smaller groups. A summary of HDBSCAN results with min_cluster_size 5–40 can be seen in the Supplementary Material.

DBSCAN was used to identify temporal clusters in the time series for each of the regions classified using HDBSCAN. Temporal anomalies represent periods when the number of events within the HDBSCAN cluster region is unusually high, which is a proxy for earthquake sequences such as mainshock-aftershock (MS-AS), foreshock-mainshock aftershock (FS-MS-AS) or swarm sequences. The date and time of each event was converted to numeric using matplotlib’s date2num() with the epoch set to “1990-01-01T00:00:00,” where “1990-01-02T00:00:00” = 1. The DBSCAN method requires two input parameters, the minimum distance between points (ɛ), and the minimum number of points (MinPts) accepted to assign a cluster. The distance metric ɛ selected is then based on the number of days, converted to numeric time. To identify temporal anomalies, we employ a rolling average window for a 30 days time period and calculate the number of events in each window, which is stepped forward in daily increments. The z-score is then calculated for each rolling average and the inverse cumulative distribution function (PPF) of the standard normal distribution is calculated to define the anomaly threshold. This threshold corresponds to a 99% confidence interval and can be interpreted as the number of events above which rolling averages are considered anomalous, or the earthquake rate above usual background levels, which is considered anomalous. We found that the number of events considered anomalous within a 30-day window was 25 events. Using a 30-day window allows for a sufficient number of data points within each window, whilst facilitating the identification of longer-term patterns, which represents the nature of earthquake swarms. 25 events within a 30-day window also generally corresponds to the earthquake counts to the west of Mt. Ruapehu preceding the 1995/1996 eruption sequence (Hurst and McGinty, 1999; Sherburn et al., 1999; Hurst et al., 2018). The selected DBSCAN parameters are then applied to each of the geospatial clusters identified using HDBSCAN.

Using an unsupervised machine learning clustering approach, such as DBSCAN, is advantageous as it has the ability to identify moments in the time series where there may be anomalously high numbers of earthquakes over a period of time, without prior domain knowledge of earthquake sequences, and it also removes the need to set absolute boundaries on the length of time for each sequence. Meaning one can automatically capture earthquake sequences over a range of time-frames without bias, providing the minimum threshold set by the z-score is achieved.

Evaluating the seismic energy associated with geospatial and temporal anomalies is an effective method for gaining a deeper understanding of the underlying characteristics and mechanisms at play. It is noteworthy that the seismic energy release is intrinsically linked to the magnitude of the earthquake event (McGarr, 1976). Methods for determining the earthquake local magnitude M _L can vary depending on the geology of the region and seismic network administration (Bormann and Di Giacomo, 2011). When performing seismic hazard assessments, it is convenient to convert local magnitude estimates to be consistent with moment magnitude M _w calculations (Rhoades et al., 2021), which is widely considered the best estimation tool for calculating earthquake magnitude (Ristau, 2009). In New Zealand, the relationship between M _L and M _w for shallow earthquakes (i.e., < 33 km depth) (Ristau, 2009), can be described as: M L = 0.88 ± 0.03 M w + 0.73 ± 0.20

Seismic moment M ₀ gives information regarding the total energy release of an earthquake and is proportional to the moment magnitude M _w defined as: M w = 3 2 log 10 M 0 − 9.0

Where M₀ is in N m (Nm) and 9.0 is an empirical calculation of the shear stress of the crust. The cumulative seismic moment of earthquake sequences has been identified as a tool for volcanic eruption forecasting (e.g., Thelen et al., 2010; White and McCausland, 2016; Meyer et al., 2021). These kinds of approaches represent a significant step towards statistically quantifying key volcanic processes, including swarm-like earthquake sequences, departing from conventional descriptive measures and providing a more data-driven foundation for volcanic eruption forecasting.

The HDBSCAN algorithm took 0.05 s to run for the earthquake catalogue containing 2,795 events, identifying six distinct earthquake regions (clusters), for a cluster size minimum of 32 samples per cluster, amounting 63% of earthquakes in the catalogue. The clusters [C0, C1, C2, C3, C4, C5] exhibited various characteristics, including geographical distribution, density, and shape (Figure 4). The algorithm performed well at grouping a variety of cluster densities and shapes, and was also effective at discerning between boundaries of high levels of seismic activity, particularly between C5 and C1 (Figure 4). The location of earthquakes is spread non-uniformly across the study area west of Mt. Ruapehu, with occasional regions of geospatial density, particularly > 5 km northwest of the summit. Earthquake activity decreases north of the summit towards Tongariro National Park. There is no clear correlation between earthquake activity and the active fault lines, with many earthquakes occurring away from documented fault lines, particularly in the region on the west flank of Mt. Ruapehu within 10 km of the summit, although, it is possible that not all active fault lines are documented. There are also few earthquakes located in close proximity to the summit.

Cluster C0 (blue) in Figure 4 is a comparatively low-density cluster, located 5–10 km south of the summit on the eastern limit of the Ohakune Fault, which runs west towards the Tohunga Junction. It contained the fewest number of earthquakes (113 events), with typical magnitudes ranging from 1.75 to 2.16, at depths typically ranging from 14 to 17 km. C1 (green) and C5 (pink) sit in a trajectory orientated approximately 40° northeast, with C1 earthquakes occurring away from the Raurimu Fault at a greater depth, on average between 14.7 and 18.1 km, compared to C5, which typically occur within the range of 11.2–14.6 km below the surface. There is a distinct difference in the magnitude of earthquakes, with significantly more M _L ≥ 3 occurring in C5 (Table 1). There is a subtle difference in the distribution density of events between C1 and C5, with C5 having a higher density of events. The cluster density tends to decrease as earthquakes are detected away from the Raurimu Fault to the west. C2 (purple), situated > 20 km from the Mt. Ruapehu summit, is in proximity to the Mangamaire Fault. Given the cluster has relatively few events, event magnitudes are characterised by some relatively large, M _L ≥ 3 , events, coinciding with a range of smaller 1.73–2.28 M _L events. Although some of the larger events occur along the fault line, there is evidence of a group of events spreading northeast of the fault. C3 (orange) is mainly situated on the junction between the north extent of the Raurimu Fault and the southwest limit of the National Park Fault, but also extents northwest, away from the Raurimu Fault. C3 has the most number of events of all clusters, with 666 earthquakes at an average depth of 12.8 km and 16 events with a M _L > 3.0. C4 (cyan) is located 12 km northwest of Mt. Ruapehu at the eroded extinct volcano Hauhungatahi (Cameron et al., 2010) (Figure 4).

DBSCAN took < 10 s to iterate over the seven clusters identified by HDBSCAN (Figure 5), including the noise cluster containing background seismicity (Figure 6). Using the z-score threshold outlined in the methodology, and a 99% confidence interval, the number of events per cluster above the threshold was calculated as 25 events. This means, for all core points within a 30 day period, DBSCAN checks for any other events within 30 days. These are border events and are added to the cluster, but not considered density reachable for any other events. If there are at least 25 core and border events in total, a cluster is formed representing a period of earthquakes above the background level. For an ɛ distance of 30 days and a minimum number of 25 events per cluster, the DBSCAN algorithm identified eight temporal clusters within the time series derived from the HDBSCAN earthquake regions. The clusters identified by HDBSCAN are labelled by adding the year of when the sequence began to the original cluster name. The occurrence of earthquake sequences identified by DBSCAN are unique to each HDBSCAN cluster. Each cluster “lights up” with activity at different times throughout the 33 years and do not seem to have any correlation with any other earthquake sequence. Some clusters have near constant background activity with occasional bursts in earthquake rate (e.g., C3-2003 and C3-2009), whereas others have a single moment of activity before returning to quiescence (e.g., C0-2006). Earthquake sequences detected above the anomaly threshold vary in duration, but tend to last no longer than 100 days. The temporal clusters shown in the time series graphs (Figure 5), including the background seismicity (Figure 6), were plotted to show the geographical extent of each earthquake sequence (Figure 7).

C0-2006 is situated in a typically inactive region and had a significant increase in earthquake rate in 2006, where 71 events occurred. Activity peaked between 23 August–2 September 2006, with 40 events. On 28 August, a total of 16 events occurred in a single day. The earthquake rate in this period is correlated with a sudden increase in cumulative seismic moment. DBSCAN automatically highlighted a period between 23 August–8 October 2006 as an anomalous earthquake sequence, containing a total of 45 events, with an average magnitude of 2.04 M _L and two events M _L > 3. Overall, DBSCAN attributed 39.8% of earthquakes within C0 to a single earthquake sequence over the 33 year period, with an average earthquake rate of 1.1 events per day.

The time series for the C1 region shows near constant earthquake activity, with an increase of activity between 2004 and 2012, also corresponding to an increase in cumulative seismic moment during the same time frame, before flattening out again. Within the period 2004–2012, when the cumulative seismic rate increases, DBSCAN identified one period (C1-2008) where the earthquake rate surpassed the background seismicity threshold. The sequence consisted of 27 small magnitude earthquakes ranging from 1.59 to 2.53 M _L , lasting 48 days from 6 April 2008–24 May 2008 at a rate of 0.56 earthquakes per day.

There was no seismicity recorded in the C2 region before 2004, perhaps due to a limitation of seismic apparatus available. Since 2014, the seismicity has increased in the region significantly, with 92.6% of total events in the cluster region occurring after July 2014 with a relatively deep source and getting shallower through time. C2-2015 had an average earthquake rate of 0.44 per day, which was deemed higher than the threshold by DBSCAN from 2 July 2015–8 October 2015. During this period, three earthquakes M _L > 3 were recorded including a 3.54 M _L , which is reflected in the sudden increase in cumulative seismic moment.

The C3 time series is characterised as near-constant earthquake activity at relatively shallow depths of < 16 km, with occasional temporal bursts in activity. There are two occasions where DBSCAN detected an increase in seismic rate greater than the background threshold. The first sequence (C3-2003) consists of 34 events at a rate of 0.4 earthquakes per day, located laterally across the north extent of the Raurimu fault, with the largest event occurring long the fault line and smaller events moving east and west. The second sequence (C3-2009) was made up of 60 events, mainly located in the region west of the fault with an approximate 20° orientation. This sequence had a higher earthquake rate than the previous sequence at approximately 0.7 events per day.

The small, dense cluster located in proximity to Hauhungatahi (C4), had three time periods where there was a sequence of earthquakes and small increases in cumulative seismic moment. The first sequence occurred in 2002 at shallow depths from 8 to 12 km and was not identified by DBSCAN. The second sequence (C4-2005), which was identified by DBSCAN due to an increase in earthquake rate of 0.45 events per day over 76 days, began in January 2005 and consisted of 34 relatively small events at depths of 13 km on average. There was a third sequence, similar to the events in 2002, originating from deeper sources (15–19 km), which began with a 3.42 M _L earthquake before following a typical mainshock-aftershock pattern.

C5 is contained to the Raurimu fault in close proximity to the town of Erua. Although earthquake activity can be observed throughout the time series, the majority of events have occurred before 2004 within C5-1995 and C5-2001, with epicentres located in a narrow band between 10 and 15 km depth. The two earthquake sequences highlighted by DBSCAN had anomalous earthquake rates, consistent with sudden increases in cumulative seismic moment. C5-1995 and C5-2001 are similar in geographical orientation, average depths, earthquake rate, and maximum magnitude. However, there were more events M _L > 3 in C5-1995 and a slightly lower average rate of 0.4 earthquakes per day. C5-1995 sequence occurred in close temporal proximity to the onset of the 1995 major eruption. Whereas the C5-2001 earthquake sequence began on 30 December 2000, lasting 88 days with an average daily rate of 0.45 events, and has no temporal proximity to any volcanic activity at Mt. Ruapehu. The number of earthquakes during the peak of seismic activity was different for C5-1995 and C5-2001, with the peak of the C5-1995 sequence totalling 16 events over 15 days, and 26 events over 25 days in total for C5-2001. However, although the largest event in each sequence had a similar magnitude (3.4 M _L ), the mean magnitude for C5-1995 was higher (2.6 M _L ) than C5-2001 (2.2 M _L ).

The slow increase of cumulative seismic moment across the entire region is coupled with occasional rapid increases (Figure 6), particularly leading up to the 1995 eruption. We have labelled this earthquake activity as background seismicity, due to the fact that earthquakes occurred in regions with no geospatial relationship. The largest increase in cumulative seismic moment began in 1995, driven by the first earthquake sequence C5-1995, before the onset of the 1995 major eruption and waned during the main eruption sequence before increasing again before and during the small eruption in 1997. There were also smaller increases in 2003/2004, which occur due to C3-2003, and slightly larger increase in 2015/2016 which coincides C2-2015.

The earthquake events occurring within C0-2006, situated along the eastern boundary of the Ohakune Fault on the southernmost region of the research domain, coincided closely in time with the phreatic eruption of October 2006. This synchronicity may hold significance due to the surprising nature of the eruption, which seemingly erupted without any warning signs of imminent volcanic activity (Jolly et al., 2010). A cursory examination of the earthquake catalog through visual analysis might not yield any apparent indicators of abnormal seismic activity in this area. Hence, the identification of the earthquake sequence from HDBSCAN and the temporal sequence from DBSCAN, was as an unforeseen revelation. C0-2006 was indicative of a complex sequence, exhibiting both MS-AS and swarm behaviour in a short period of time (Figure 8). The largest event (3.1 M _L ) on 26 August 2006 was followed by two successively smaller events following a typical mainshock-aftershock sequence. On 28 August 2006, 16 individual events were recorded with the largest and second largest magnitude of 3.03 M _L and 2.9 M _L . The characteristics of the differences in magnitude and the cumulative seismic moment over a short period of time coupled with a b-value of 1.13, suggests this earthquake sequence is of swarm origin and may be the result of an intrusion of fluids into a dike at shallow depths, which interacted with the hydrothermal system and the crater lake. The sudden intrusion of magma into a crack near the surface may have triggered a rapid expansion and steam driven eruption without necessarily raising the crater temperature.

This cluster of earthquakes is the only set which occurs entirely away from a mapped fault source. C1-2008 is situated > 1 km west of the Raurimu fault and is clearly a separate signal to the neighbouring C5-1995 and C5-2001. A near-constant seismic sequence began around the year 2000 at depths > 12 km and has continued until 2023, with the most active phase from 2005 to 2012. The earthquake depth is most likely linked to the angle of the subducting plate, which is probably why the earthquakes in C5 occur at a shallower depth on average (e.g., Keats et al., 2011). The b-value for C1-2008 was the highest in the region, but also had the largest uncertainty (Table 2), which is probably driven by the small sample size.

C2-2015, on the northwest extent of the study area, may be an artifact of the timing of the station deployment. The earthquake rate increases to the threshold defined in DBSCAN, to a peak activity from 24 August 2015, where there was an increase in magnitude to the 3.31 M _L event on 26 August. The C2-2015 distribution of events is located away from the fault with a b-value of 0.99 (Figure 7). The temporal evolution of C2-2015 develops from a series of earthquakes with a narrow magnitude band, to a gradual increase in magnitude and rate, with the largest event occurring in the middle of the sequence (Figure 9). Extrapolating from the C2-2015 sequence, there is clear evidence of average earthquake source depths from deep to shallow (Figure 5C), which could be evidence of fluid migration through the mid-crust (Yoshida and Hasegawa, 2018; Yoshida et al., 2023).

The C3 cluster is the most seismically active region in the study area. The near-constant rates with occasional bursts is similar to the seismicity described by Hayes (2004) to the southeast at Waiouru, which is outside of our study area. The first and second temporal bursts in C3-2003 and C3-2009 occur on different trajectories across the National Park fault, yet are similar in geographical extent (Figure 7). C3-2003 and C3-2009 exhibit similar swarm-like behavior (Figure 10), with high b-values of 1.51 and 2.05 respectively. We observe multiple temporal bursts in seismicity, which did not reach the threshold set for DBSCAN, and a lower threshold would certainly yield more swarm events. The temporal evolution of the C3 cluster has similar characteristics to swarm bursts discussed by Petersen and Pankow (2023). Here the variable inter-event bursts are related to variability in fluid migration in a fault-valve model (Sibson, 1992). For C3, changes in pressure from fluids below the surface, coupled with precipitation from the hydrothermal system at Mt. Ruapehu (Christenson and Wood, 1993), can cause the fault to open and close, leading to unpredictable variation in inter-event times (Petersen and Pankow, 2023). The fluid-driven source is supported by relatively high b-values for the region. Further sensitivity analysis into the time series for C3 might offer more insight into the processes we see in this region; however that is beyond the scope of this paper.

C4-2005 was detected in close proximity to the extinct volcano, Hauhungatahi (Cameron et al., 2010), which is made up of > 170 ka old eruptive material. The small earthquake sequence detected at the end of 2004 is most dissimilar from the other signals, exhibiting earthquakes with a narrow band of magnitudes typically ≤ Δ0.5 M _L with no significant energy release (Figure 11). C4-2005 has the high estimated b-value of 2.25, which is perhaps driven by a small sample size (e.g., Geffers et al., 2022) coupled with a narrow magnitude band and the majority of earthquakes < 3 M _L .

C5, situated along the northern extent of the Raurimu fault, close to the small town of Erua, is a region with earthquake activity which is often referred to in the literature as a precursor to the major 1995 eruption (Hurst and Vandemeulebrouck, 1996; Hurst and McGinty, 1999; Hurst et al., 2018). Although, the cluster is in the orientation, and in close proximity to the Raurimu fault, seismicity is thought to have not originated from the fault (Keats et al., 2011). The first earthquake sequence, C5-1995, began in January 1995 and ended at the end of July 1995. There were two bursts of activity, a small MS-AS sequence occurred mid-March, before the main burst of activity which began early April. The April burst had 16 individual earthquake events with 2 M _L > 3 and coincided with a warming of the crater lake along with an increase in Mg2+ ions (Nakagawa et al., 1999), which may be an indication of fresh magma interacting with the hydrothermal system (Hurst and Vandemeulebrouck, 1996). The earthquake rate then waned until the beginning of the major eruption sequence on 17 September 1995 (Figure 12A). Conversely, the b-value of the DBSCAN sequence is 0.88, which is comparatively low for earthquake swarms in the region and is consistent with Hurst and McGinty (1999), and the tectonic mechanisms occurring to the southeast during the same time-frame (Hayes, 2004).

The second burst of activity (C5-2001), which began in January 2001, had a similar orientation, geographical extent and total number of events as the 1995 sequence, however did not precede any kind of volcanic activity. C5-2001 had a different temporal evolution to the 1995 sequence, with a buildup of seismic activity preceding the largest event on 3 March 2001, and smaller events following (Figure 12B). The largest event occurring around the middle phase of the temporal evolution is a key characteristic of swarm-like behaviour associated with volcano-tectonic events (Jones, 2005; Pesicek et al., 2018). The difference between the largest two events, 3.45 ± 0.23 and 2.85 ± 0.23, is 0.60 ± 0.46, and perhaps is on the boundary of swarm-like magnitude differences. The b-value of the swarm is calculated as 1.67 which is typical for swarm behaviour driven by volcano-tectonic processes (Gudmundsson, 2020).

The cumulative seismic moment of an earthquake sequence is a proxy for the total amount of energy released and provides insights into the mechanisms driving the seismicity in the region. If there is a lack of pressure or insufficient volume, magma may not have enough energy to overcome stress barriers and can become arrested (Caricchi et al., 2021), which may provide a reason for why there was not an eruption following the earthquake swarm in 2001.

The results of this study have significant implications for understanding the seismic behavior the west of Mt. Ruapehu. By employing density-based clustering algorithms, we have gained a more comprehensive perspective on the spatial and temporal patterns of earthquake activity. Spatial and temporal clustering has shown promise as a quick and effective way of detecting earthquake swarms, particularly within noisy data sets. Retrospective analysis of earthquake catalogues containing many thousands of events can be filtered into high-density regions using a couple of lines of code to call the HDBSCAN algorithm without the necessity of prior domain knowledge. These regions can then be analysed for temporal bursts in seismic activity by using DBSCAN, using the ɛ distance function to calculate distance as time. This is useful for detecting swarms; however it requires the statistical calculation of minimum events for a given time window, which can vary for different systems. Furthermore, the parameter selection for the HDBSCAN can be very sensitive, producing very different results for different minimum cluster size parameters. However, the ability to automatically identify spikes in earthquake rate over different periods of time revealed interesting results, which were less sensitive to the initial parameters, and correlated well with the cumulative seismic rate. Further research could explore including magnitude to the temporal analysis to distinguish between swarm-like and MS-AS sequences. This may help focus on further refining our understanding of the relationships between seismic clusters, volcanic processes, and eruption precursors. Additionally, the application of advanced machine learning may provide deeper insights into the complex interplay between seismicity and volcanic activity in the Mt. Ruapehu region.