# Seasonal Snow Cycles and Their Possible Influence on Seismic Velocity Changes and Eruptive Activity at Ruapehu Volcano, New Zealand

Alexander S Yates, Corentin Caudron, Aurélien Mordret, Philippe Lesage, Virginie Pinel, Thomas Lecocq, Craig A Miller, Oliver D Lamb, Nicolas Fournier

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

To cite this version:

Alexander S Yates, Corentin Caudron, Aurélien Mordret, Philippe Lesage, Virginie Pinel, et al. "Seasonal Snow Cycles and Their Possible Influence on Seismic Velocity Changes and Eruptive Activity at Ruapehu Volcano, New Zealand." Journal of Geophysical Research: Solid Earth, 2024, 129 (12), 10.1029/2024jb029568. ird-04995397

HAL Id: ird-04995397

https://ird.hal.science/ird-04995397v1

Submitted on 18 Mar 2025

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

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

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

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

## RESEARCH ARTICLE

10.1029/2024JB029568

**Key Points:**

* Evidence of preferential timing in eruption chronology at Ruapehu, with higher than expected proportion of eruptions in spring

* Seismic velocities computed via coda-wave interferometry, using cross-wavelet transform approach to better handle variable tremor

* Groundwater changes and snow-loading influence seasonal velocities, with the latter possibly a mechanism to trigger gas-driven eruptions

**Supporting Information:**

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

**Correspondence to:**

A. S. Yates,
alecyates1991@gmail.com

**Citation:**

Yates, A. S., Caudron, C., Mordret, A., Lesage, P., Pinel, V., Lecocq, T., et al. (2024). Seasonal snow cycles and their possible influence on seismic velocity changes and eruptive activity at Ruapehu volcano, New Zealand. *Journal of Geophysical Research: Solid Earth*, 129, e2024JB029568. https://doi.org/10.1029/2024JB029568

Received 23 MAY 2024

Accepted 28 NOV 2024

**Author Contributions:**

**Conceptualization:** Alexander S. Yates, Corentin Caudron, Aurélien Mordret, Philippe Lesage

**Data curation:** Alexander S. Yates

**Formal analysis:** Alexander S. Yates

**Funding acquisition:** Corentin Caudron, Philippe Lesage

**Investigation:** Alexander S. Yates

**Methodology:** Alexander S. Yates, Corentin Caudron, Aurélien Mordret, Philippe Lesage, Virginie Pinel, Thomas Lecocq

**Project administration:**

Corentin Caudron, Aurélien Mordret, Philippe Lesage

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

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

## Seasonal Snow Cycles and Their Possible Influence on Seismic Velocity Changes and Eruptive Activity at Ruapehu Volcano, New Zealand

Alexander S. Yates<fnref n="1" />,<fnref n="2" />,³ Ⓓ, Corentin Caudron<fnref n="1" />,<fnref n="2" />,³ Ⓓ, Aurélien Mordret<fnref n="1" />,⁴ Ⓓ, Philippe Lesage<fnref n="1" /> Ⓓ, Virginie Pinel<fnref n="1" /> Ⓓ, Thomas Lecocq⁵ Ⓓ, Craig A. Miller⁶ Ⓓ, Oliver D. Lamb⁶ Ⓓ, and Nicolas Fournier⁶ Ⓓ

**Abstract** Understanding volcanic eruption triggers is critical toward anticipating future activity. While internal magma dynamics typically receive more attention, the influence of external processes remains less understood. In this context, we explore the relationship between seasonal snow cycles and eruptive activity at Ruapehu, New Zealand. This is motivated by apparent seasonality in the eruptive record, where a higher than expected proportion of eruptions (post-1960) occur in spring (including the two previous eruptions of 2006 and 2007). Employing recent advancements in passive seismic interferometry, we compute sub-surface seismic velocity changes between 2005 and 2009 using the cross-wavelet transform approach. Stations on the volcano record a higher velocity in winter, closely correlated with the presence of snow. Inverting for depth suggests these changes occur within the upper 300 m. Notably, we observe that the timing of the previous two eruptions coincides with a period associated with an earlier velocity decrease at approximately 200–300 m depth relative to the surface. Reduced water infiltration (as precipitation falls as snow) is considered a likely control of seasonal velocities, while modeling also points to a contribution from snow-loading. We hypothesize that this latter process may play a role toward explaining seasonality in the eruptive record. Our findings shed light on the complex interactions between volcanoes and external environmental processes, highlighting the need for more focused research in this area. Pursuing this line of inquiry has significant implications toward improved risk and hazard assessments at not just Ruapehu, but also other volcanoes globally that experience seasonal snow cover.

**Plain Language Summary** Understanding what triggers volcanic eruptions is crucial for anticipating future activity. Previous studies often focus on processes occurring within the volcano; for example, movements of magma. However, there is an increasing need to understand how external, environmental, processes interact with volcanoes. This is apparent at Ruapehu volcano in New Zealand, where we find evidence that more eruptions occur in spring relative to other seasons. To further explore this, we measure small changes in the speed of seismic waves traveling in the vicinity of the volcano. This provides information on processes occurring beneath the surface, with the wave speed (velocity) sensitive to subtle changes in the rock they pass through. Our findings indicate that the presence of seasonal snow is the main driver of near-surface velocity changes. This occurs as rain falls as snow during winter, and also through compaction under the weight of a snow layer. We hypothesize that small variations in subsurface pressures—through the addition and subsequent removal of a seasonal snow layer—may be responsible for increased spring-time activity. This line of inquiry requires further exploration, with significant implications for risk and hazard assessment at volcanoes globally that experience seasonal snow cover.

## 1. Introduction

It has long been recognized that volcanoes show sensitivity to external forcings. For example, tectonic activity (Nishimura, 2017; Yokoyama, 1971), tidal forces (Bredemeyer & Hansteen, 2014; Dinger et al., 2019; Jaggar et al., 1924; Neuberg, 2000), long-period sea-level changes (Mason et al., 2004; McNutt & Bevan, 1987), and precipitation (Farquharson & Amelung, 2020; Gaete et al., 2020; Hort et al., 2003; Matthews

<footnote n="1">¹University Grenoble Alpes, University of Savoie Mont Blanc, CNRS, IRD, University Gustave Eiffel, Grenoble, France,</footnote>

<footnote n="2">²Department of Geosciences, Environment and Society, Laboratoire G-Time, Université libre de Bruxelles, Bruxelles, Belgium, ³WEL Research Institute, Wavre, Belgium, ⁴Department of Geophysics and Sedimentary Basins, GEUS, Geological Survey of Denmark and Greenland, Copenhagen, Denmark, ⁵Seismology and Gravimetry Department, Royal Observatory of Belgium, Uccle, Belgium, ⁶GNS Science, Wairakei Research Centre, Taupo, New Zealand</footnote>
**Software:** Alexander S. Yates, Corentin Caudron, Aurélien Mordret, Virginie Pinel, Thomas Lecocq

**Supervision:** Corentin Caudron, Aurélien Mordret, Philippe Lesage

**Validation:** Alexander S. Yates, Corentin Caudron, Aurélien Mordret, Philippe Lesage, Virginie Pinel, Craig A. Miller, Oliver D. Lamb, Nicolas Fournier

**Visualization:** Alexander S. Yates

**Writing – original draft:** Alexander S. Yates

**Writing – review & editing:** Alexander S. Yates, Corentin Caudron, Aurélien Mordret, Philippe Lesage, Virginie Pinel, Thomas Lecocq, Craig A. Miller, Oliver D. Lamb, Nicolas Fournier

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Topographic map of Ruapehu, New Zealand. The snow cover as seen from MODIS satellite data is shown for the months July 2006 (purple), October 2006 (orange), and May 2007 (green). Different ski-stations distinguished as T (Turoa) and W (Whakapapa). Inset (bottom-right) shows location of map region (red square) within New Zealand.</figcaption></figure>

et al., 2002; Uchida, 2024) have all been linked to volcanic activity. A summary of external (and internal) eruption triggers can be found in the recent review of Caricchi et al. (2021).

Another external process that has been proposed to influence volcanic activity is the unloading of snow and ice. This can occur on longer time scales, through de-glaciation (Conway et al., 2023; Jellinek et al., 2004; Sigmundsson et al., 2010; Sigvaldason et al., 1992), and also shorter time scales due to the unloading of seasonal snow/ice (Albino et al., 2010). Studies of this phenomena have largely relied on statistical analysis of the eruptive history—sometimes combined with modeling of pressure changes within the magma chamber (e.g., Albino et al., 2010; Jellinek et al., 2004; Sigmundsson et al., 2010)—to assess the correlation between unloading and volcanic activity. With this in mind, techniques that allow us to probe the sub-surface directly would enable us to better understand the underlying mechanisms behind supposed external triggering of eruptive activity.

The last two decades have seen great advances in our ability to monitor volcanic systems using coda wave interferometry. By cross-correlating seismic energy recorded by different sensors, we can monitor subtle changes in seismic velocities that can be related to a variety of geophysical processes. For example, velocity changes recorded at volcanoes have been attributed to increased pressurization within the volcano system (Brenguier et al., 2008; Budi-Santoso & Lesage, 2016; Caudron et al., 2021; Mordret et al., 2010; Yates et al., 2019) and also changes in fluid content (Caudron et al., 2015, 2022). External factors are also regularly present in velocity changes recorded at volcanoes. These include meteorological processes, such as precipitation (e.g., Caudron et al., 2022; Rivet et al., 2015; Sens-Schönfelder & Wegler, 2006) and snow-loading/unloading (e.g., Cannata et al., 2017; Caudron et al., 2022; Donaldson et al., 2019; Hotovec-Ellis et al., 2014), and velocity decreases following large earthquakes (e.g., Battaglia et al., 2012; Lesage et al., 2014; Yates et al., 2019). For the purposes of monitoring volcanoes, velocity changes associated with environmental processes are often undesirable as they can mask changes associated with volcanic processes. Thus, several authors attempt to model the velocity changes produced by environmental effects, for example, precipitation and snow fall, and correct for this (Budi-Santoso & Lesage, 2016; Rivet et al., 2015; Wang et al., 2017). However, it is also important to examine these external processes in greater detail toward exploring their potential influence on volcanic systems.

We explore the relationship between snow cycles and seismic velocity changes at Ruapehu volcano in New Zealand (Figure 1). The stratovolcano (2,797 m summit elevation) is located at the southern end of the Taupo Volcanic Zone (Leonard et al., 2021) and accumulates snow above approximately 1,600 m between the months of
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) All events of size 0–5 at Mt Ruapehu in catalog of Scott (2013), excluding events with no listed start month. (b) De-clustered catalog (size 2+, post-1960, 1995/96 sequence removed, 14-day interval between events). (c) Circular representation of events in (b), with corresponding p value following Rayleigh test against a uniform distribution. (d) As before, but for size 2 and 3 events only.</figcaption></figure>

April and November (Manville et al., 2000). Ruapehu has experienced many eruptions over the past century, with the most recent eruptions occurring in 2006 and 2007 (both phreatic). Interestingly, both of these eruptions occurred approximately 1 year apart (on 4 October 2006 and 27 September 2007) during the southern hemisphere spring. Furthermore, evidence of preferential timing of eruptions is identified when analyzing the eruption catalog of Scott (2013) (extending back to the 1830s), with over 50% of eruptions occurring during the months of spring following de-clustering of the catalog (Figure 2; detailed in Section 2.1).

We take advantage of recent developments in passive seismic interferometry by applying the cross-wavelet transform method (e.g., Mao et al., 2020; Mordret et al., 2020) to compute single-station seismic velocity changes between 2005 and 2009. This approach provides greater time-frequency resolution compared to traditional approaches (Mao et al., 2020). Mordret et al. (2010) had previously computed velocity changes associated with the eruptions of 2006 and 2007, identifying a decrease in the days prior to the 2006 eruption. In this study, we target a broader frequency range than the 0.2–0.7 Hz range used by Mordret et al. (2010), including higher frequencies (>1 Hz) influenced by volcanic tremor. In doing so, we can explore velocity changes at a shallower depth than previously investigated. Furthermore, processing a longer and more continuous data set relative to Mordret et al. (2010) provides an opportunity to better understand the sub-surface influence of seasonal snow cycles and whether this could play a role in preferential timing of eruptive activity at the volcano.

## 2. Data and Methods

### 2.1. De-Clustered Eruption Catalog

Evidence of seasonality in the eruptive record at Ruapehu is identified through the preparation of a de-clustered eruption catalog. We begin with the revised catalog of Scott (2013) (Figure 2a), extending back to the 1830s. First, eruptive activity assigned a size of 1 (or less)—which includes small phreatic eruptions confined to the crater lake—are excluded toward ensuring a good level of completeness (see Table S1 in Supporting Information S1 for size definitions at Ruapehu). Furthermore, only eruptions post-1960 are included. This differs slightly from the 1940 cut-off considered by Scott (2013), yet appears preferable toward ensuring improved completeness for size 2
events (Figure S1 in Supporting Information S1) (i.e., phreatic or phreatomagmatic eruptions confined to the crater basin but depositing material outside the crater lake). Then, as done by Scott (2013) and Christophersen et al. (2022), we exclude activity associated with the 1995/96 eruptive episode as this period accounts for a significant proportion of eruptions in the catalog that are related to each other. Finally, since the target is eruptions that may be triggered by external processes, all eruptions occurring within 14-day of a prior eruption are removed (as done by Christophersen et al. (2022), Figure S2 in Supporting Information S1). This excludes events causatively linked to a previous event. The final de-clustered catalog consists of 50 eruptive events, shown in Figure 2b. Further details of the de-clustering process, with justifications, can be found in Text S1 in Supporting Information S1.

From the de-clustered catalog, we observe that over 50% of the eruptions occur within the spring months of September, October, and November (Figure 2b). The majority of these are small eruptions (size 2), with only one event (of 25) classified as phreatomagmatic (Figure S3 in Supporting Information S1). The remainder are described as non-magmatic (14 events), as phreatic/phreatomagmatic (6 events), or undefined (4 events). Such a high proportion of eruptions in spring months—particularly in October, where 30% of all recorded size 2 events occur—appears at odds with a uniformly distributed pattern of eruptions in time. This is supported by relatively low *p* values (*p* < 0.05) using the Rayleigh test (Mardia, 1975), which assesses whether the data cluster in a certain direction when represented as vectors on a unit circle (Figures 2c and 2d). For smaller events (sizes 2 and 3), a *p* value of 0.013 is recovered (Figure 2d) while this falls to *p* = 0.002 (Table S3 in Supporting Information S1) when considering only size 2 events. We also perform Fisher's exact test (Fisher, 1935), which calculates the probability of obtaining the observed distribution of events under the null hypothesis of independence, and achieve similarly low *p* values when comparing the number of eruptions in spring to the rest of the year (*p* ≤ 0.01 for size 2 and 3 events; Table S3 in Supporting Information S1). A full summary of these *p* values and their corresponding tests can be found in Text S2 in Supporting Information S1.

The results of both statistical tests (Rayleigh and Fisher's) indicate a possible relationship between the timing of eruptions and the months of spring. This is particularly the case for small events (sizes 2 and 3), which we might anticipate would be more susceptible to external forcing relative to larger, magmatic, events. However, we stop short of concluding the findings are statistically significant given concerns regarding the use of an arbitrary significance threshold (Di Leo & Sardanelli, 2020; Wasserstein & Lazar, 2016). Nonetheless, the low *p* values recovered suggest further investigation into preferential timing of eruptive activity at the volcano is warranted.

## 2.2. Seismic Data and Cross-Correlation Function Preparation

Three-component seismic data (100 Hz) are downloaded from a network of nine seismic stations in the proximity of the volcano (Figure 1a). These consist of seven short period sensors and two broadband sensors, with data recorded and publicly available via the GeoNet data service (GNS Science, 2021). Since a single-station approach is used to compute seismic velocity changes, cross-correlating between different components of the same station, no instrument response correction is applied.

Cross-correlation functions are processed using the python package MSNoise (Lecocq et al., 2014) following relatively standard processing briefly outlined below. We first group data into continuous chunks, before demeaning and tapering, and merging into 1-day traces. Any gaps smaller than 10 s are filled via interpolation. The data are then downsampled to 25 Hz, after applying an initial bandpass filter between 0.01 and 10 Hz to avoid aliasing. Data for each station component are then spectrally whitened between 0.1 and 8.0 Hz and clipped in the time-domain at three times the RMS. Cross-correlation functions are computed between the different components of a single-station—East-North (EN), East-Vertical (EZ), and North-Vertical (NZ)—in 30-min windows. If any windows still contain gaps (i.e., were larger than 10 s), the data are excluded. The remaining functions are then stacked to give 10-day cross-correlation functions, where the function recorded for a given day reflects the cross-correlation function of that day and the nine previous cross-correlation functions. The auto-correlations (EE, NN, ZZ) are not considered given an inability to apply 1-bit spectral whitening when using the classical cross-correlation approach. This typically results in reduced stability in the final waveforms (e.g., De Plaen et al., 2016; Hobiger et al., 2014; Yates et al., 2019).

## 2.3. Tremor as a Seismic Source

One of the challenges in monitoring seismic velocity changes at volcanoes is the presence of volcanic tremor, which can significantly impact the character of cross-correlation functions (Ballmer et al., 2013; Gómez-García
et al., 2018). Thus it is often the case that studies focus on a frequency band below the band of volcanic tremor (e.g., Cannata et al., 2017; Mordret et al., 2010)—usually dominant between 0.5 and 5.0 Hz (Chouet & Matoza, 2013)—unless it can be verified that the tremor is stable as a seismic source (e.g., Donaldson et al., 2017). At Ruapehu, tremor is common at frequencies above 1 Hz, usually with a single peak between 1.8 and 2.3 Hz (Hurst, 1998; Hurst & Steven, 1993; Latter & Balsillie, 1984; Sherburn et al., 1999). Such energy is clearly seen in Figure 3a, where we compute the spectral width of the network covariance matrix as a function of time and frequency (e.g., Seydoux, Shapiro, de Rosny, & Landès, 2016; Seydoux, Shapiro, De Rosny, et al., 2016) using data recorded by four short-period stations closest to the volcano (DRZ, TRVZ, TUVZ, and WNVZ). This follows the processing of Yates et al. (2023) using the python package CovSeisNet (Covseisnet, 2020), with 1-bit temporal and spectral normalization applied during pre-processing and the covariance matrix estimated in 5,000 s windows. The spectral width reflects the eigenvalue distribution of the covariance matrix and tends to zero when the seismic wavefield at the corresponding frequency is dominated by a single source. Similarly, it approaches a maximum when non-coherent noise is considered (Seydoux, Shapiro, De Rosny, et al., 2016). Thus, reduced spectral width values indicate a more coherent seismic wavefield dominated by fewer seismic sources. Further details can be found in previous studies applying this methodology (e.g., Journeau et al., 2020; Seydoux, Shapiro, de Rosny, & Landès, 2016; Seydoux, Shapiro, De Rosny, et al., 2016; Soubestre et al., 2018; Yates et al., 2023).

It has previously been suggested that volcanic tremor at Ruapehu may be useful as a repeating seismic source toward passive seismic interferometry (Yates et al., 2023). Using hierarchical clustering, it was shown that many cross-correlation functions computed between stations FWVZ and WPVZ were relatively similar in the time-domain and their spectra within the frequency band of tremor (1.8–4.0 Hz) (Yates et al., 2023). Thus we target this energy toward computing velocity changes between 2005 and 2009, where tremor was more regular relative to later years (Figure S5 in Supporting Information S1). The reliability of such an approach, in relation to our results, is discussed further in Section 4.1.

## 2.4. Measuring Velocity Changes Using Wavelet Transform

In this study we apply the cross-wavelet transform approach, as applied by Mao et al. (2020). The method is based on wavelet cross-spectrum analysis (Farge, 1992; Torrence & Compo, 1998), and has the advantage that it allows us to obtain delay time measurements between current and reference stacks of cross-correlation functions with a finer spectral resolution than can be obtained using traditional methods (Mao et al., 2020; Yuan et al., 2021); for example, the doublet (Clarke et al., 2011; Poupinet et al., 1984) or stretching (Sens-Schönfelder & Wegler, 2006) method. Furthermore, the ability to compute phase differences between two waveforms in the time-frequency space allows more precise quality control criteria. Thus, where it might ordinarily be necessary to discard entire cross-correlation functions in the presence of a changing seismic source (e.g., tremor), the wavelet transform allows for the exclusion of specific frequencies in time while retaining reliable data. Such flexibility makes it an ideal choice for processing seismic velocity changes at Ruapehu, where subtle changes in the seismic source (e.g., tremor) may influence the corresponding velocity change measurements. Note, additional wavelet-domain methods for measuring velocity changes exist—in addition to that applied by Mao et al. (2020)—but these are not evaluated here.

Velocity changes are recovered by measuring the travel time shifts between the 10-day (current) stacks and a reference stack of all available cross-correlation functions. Following Mao et al. (2020), the Continuous Wavelet Transform (CWT) of both the reference and current stacks are obtained in the time-frequency domain by performing convolutions between the cross-correlation functions and a series of scaled and translated versions of a mother wavelet. Here, we use the Morlet wavelet (Morlet et al., 1982a, 1982b) for the mother wavelet and compute the wavelet transform between 0.1 and 8.0 Hz. The CWT is complex and thus contains both amplitude and phase information. Thus it is possible to recover the phase difference between the wavelet transform representations of both the reference ($W_X(s, n)$) and current ($W_Y(s, n)$) stacks (where $s$ is wavelet scaling and $n$ the translation). To do this, the cross-wavelet transform ($W_{XY}(s, n)$) is computed between the wavelet representation of both waveforms. The phase difference is then extracted as $\Delta\phi = \arg(W_{XY}(s, n)) = \arg(W_X(s, n)) - \arg(W_Y(s, n))$. This can be divided by $2\pi f$ at each frequency $f$ to give a delay time measurement (Figure S6a in Supporting Information S1). Note, to enhance the signal, we apply an additional smoothing to both $W_X$ and $W_Y$ via a Gaussian filter in the frequency domain (temporal smoothing with standard deviation 1.5 s) followed by a box-car filter across wavelet scales (set to five scales). Finally, apparent velocity changes are recovered by performing a weighted linear regression of delay time against lag time within a chosen coda window between 0.1 and 4.0 Hz, where
<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a) Normalized spectral width measurement. Lower values indicate a more coherent seismic wavefield dominated by fewer seismic sources. Computed using stations DRZ, TRVZ, TUVZ, WNVZ. Eruption timing indicated by red-dashed lines. (b) Velocity changes (average of EN, EZ, and NZ component-pairs) measured at DRZ station. (c) Velocity changes at FWVZ station. (d) Velocity changes for MTVZ station. (e) Wavelet coherence (average of EN, EZ, and NZ component-pairs) between 10-day and reference stacks for DRZ station. (f) Wavelet coherence for FWVZ station. (g) Wavelet coherence for MTVZ station. (h) Snow radius and thickness.</figcaption></figure>
$dt/t = -dv/v$. The error in the slope, and thus $dv/v$, is computed as the square root of the diagonal element of the covariance matrix resulting from the regression (based on Scipy function scipy.optimize.curve_fit; Virtanen et al., 2020). This gives the standard deviation of the slope estimate. We define the minimum lag time for the coda as 5 s and define the length of the coda window to be 20 cycles at a given frequency (up to a maximum of 120 s for the very lowest frequencies). For example, at 1 Hz, the coda window would be from 5 to 25 s. The weighting function (Figure S6d in Supporting Information S1) used in the regression is based on the amplitude of the wavelet cross-spectrum $|WXY(s,n)|$ (shown in Figure S6c in Supporting Information S1), which reflects the local power of both the current and reference stack (Mao et al., 2020). Specifically, we normalize the amplitude of the wavelet cross-spectrum by dividing it by the maximum value along the time-axis at each frequency, resulting in weights that vary between 0 and 1 (following Mao et al., 2020).

In addition to weighting the linear regression based on the amplitude of the wavelet cross-spectrum, we apply additional control quality criteria to ensure that only reliable measurements of velocity change are recorded. Delay time measures with a wavelet coherence less than 0.7, or maximum delay time measurement of 0.3 s, are assigned a weight of zero. Furthermore, we require that at least 25% of the delay time measurements have a non-zero weight, to avoid spurious velocity change estimates based on only a few points. The influence of wavelet coherence threshold is discussed further in Section 4.1.

As mentioned, an advantage of applying the wavelet transform is the ability to retain reliable measurements at certain frequencies while discarding unreliable measurements that relate to varying seismic sources. For example, during time periods where there is minimal tremor signal, the wavelet coherence and weighting function have very low values above 1 Hz (Figure S7 in Supporting Information S1). In this case, the measurements below 1 Hz would be retained and those above 1 Hz discarded (Figure S7e in Supporting Information S1). Similarly, we can observe that sometimes the coherence is reduced at only certain frequencies (Figure S8 in Supporting Information S1), likely corresponding to changes in the tremor source. In this case, the weighting assigned to the majority of measurements above approximately 1.5 Hz is approximately zero (Figure S8d in Supporting Information S1). Thus, using the wavelet approach allows coherent measurements at specific frequencies to be preserved while disregarding measurements where changes in the source characteristics are occurring.

## 2.5. Complementary Data Sets

Snow coverage, as shown in Figure 1a, is estimated from MODIS (Moderate Resolution Imaging Spectroradiometer) satellite data (Hall et al., 2006). After picking the boundary of snow, the area is converted to a radius assuming a circular load (Figure 1b). This represents a simplified approximation, but serves to quantify the coverage. The radius of snow cover for each year, and their average, is shown in Figure S9a in Supporting Information S1.

Snow thicknesses are recorded at Whakapapa and Turoa ski stations—both approximately 2,000 m elevation—during times the stations are operational (data provided by Ruapehu Alpine Lifts Ltd). The Whakapapa measurement is made approximately 1 km north-east of FWVZ station (Figure 1). The Turoa measurement is made in close proximity to TRVZ station. The average yearly trends recorded at both sites between 2005 and 2009 are similar (Figure S9b in Supporting Information S1).

Daily precipitation records obtained from New Zealand's National Climate Database (NIWA, 2020), where we download data from the Chateau weather station located approximately one km north of WPVZ station (approximately 10 km from the volcano summit). These data are then used to compute pore pressure changes associated with fluid diffusion after rainfall events, following Rivet et al. (2015) and Wang et al. (2017). This assumes that groundwater level changes can be related to the deviation of rainfall from the average precipitation rate. While it is anticipated that precipitation will vary over the study region, thus limiting a more precise comparison at individual stations, we anticipate that the longer-term trends will be similar and thus sufficient to evaluate a possible correlation between rainfall and seismic velocities.

Based on the one-dimensional fully coupled diffusion equation of Talwani et al. (2007), and applied by Rivet et al. (2015), pore pressure changes ($P(r,t)$) are computed at a distance $r$ and time $t$, following:
$$P(r, t) = \sum_{i=1}^{n} \delta p_i \operatorname{erfc} \left[ r / (4c(n-i)\delta t)^{1/2} \right] \quad (1)$$

where r refers to the distance (treated as the depth), n the total number of time increments δt between the start of the time series (i = 1) and time t, δpᵢ is the precipitation load change computed as ρgδhᵢ, where δhᵢ is the change in precipitation recorded at time tᵢ from the previous day (tᵢ-1), and c is the diffusion rate in m²/s. We use c = 1m²/s, following Wang et al. (2017) and Caudron et al. (2022), and compute the mean pore pressure in the upmost 4 km (averaging the results from r = 0 m to r = 4000 m, in 500 m steps) using precipitation data starting from the year 2000. The 4 km limit represents an approximate sensitivity of the lowest frequency used (0.1 Hz). The final calculated mean pore pressure changes between 2005 and 2009 are shown in Figure S10 in Supporting Information S1.

# 3. Results

## 3.1. Measured Velocity Changes

Velocity changes recorded at three different stations are shown in Figure 3, where we average the results of the three component-pairs (EN, EZ, and NZ) for each station. Two of these stations, DRZ and FWVZ (Figures 3b and 3c respectively), correspond to the stations closest to the summit of the volcano and thus are expected to have significant snow cover during winter months. The third station, MTVZ (Figure 3d), is located at a greater distance from the volcano, where snow cover throughout the year is minimal (Figure 1a). The wavelet coherence computed between 10-day stacks and the reference function is also shown for each station (Figures 3e-3g). This confirms the dominance of tremor above 1 Hz across the network, with higher coherence measured at the furthest station (MTVZ, Figure 3g) during stronger tremor episodes (Figure 3a). Similarly, reduced tremor is associated with low coherence.

At all three stations, alternating patterns of low and high velocities dominate, with an apparent yearly cycle. The timing of these changes differs, however, between the stations on the volcano and the station at distance from the volcano. At DRZ and FWVZ stations, velocities are generally highest at between the months of July and October (Figures 3b and 3c respectively). At MTVZ, however, high velocities occur from approximately February to July (Figure 3d). We also see that velocities are relatively consistent across the entire frequency range at MTVZ station, especially visible during 2008. In contrast, alternating low and high velocities at DRZ and FWVZ station are mostly seen at frequencies above 1 Hz. Furthermore, a significant proportion of data above 2.5 Hz is excluded at stations DRZ and FWVZ (Figures 3b and 3c) due to low coherence (Figures 3e and 3f). At DRZ, in particular, the timing of low coherence appears to match the periods of snow cover (Figure 3h).

Computing the mean velocity change between 1 and 2.5 Hz (targeting the annual fluctuations) at DRZ and FWVZ station shows that velocities in this frequency band appear correlated with the presence of snow on the volcano. For example, the timing of seasonal velocity increases, and subsequent decreases, closely follows the trend of the approximate radius of snow cover (dashed-line in Figure 4a). Furthermore, the highest velocities are generally recorded during periods of greatest snow thickness each year (shaded gray in Figure 4). At MTVZ, however, velocities tend to be slightly lower during periods of more snow (Figure 4b). Furthermore, a good fit is observed when comparing velocity changes recorded at MTVZ with pore pressure changes computed using daily precipitation data (Figure 4b).

Velocities at each station can be stacked (averaged) over different calendar years to better visualize annual changes (Figure 5), with an average computed after dividing by the total number of years (four). This is similarly performed for the wavelet coherence also (Figure S11 in Supporting Information S1). We focus on the frequency band between 0.25 and 2.5 Hz, given the lack of reliable measurements above 2.5 Hz following quality control (e.g., Figures 3b and 3c). Again, we observe differences in the timing of the seasonal trend at different stations. This is reflected when plotting the approximate timing of the seasonal velocity high (at ~2 Hz) against station elevation (Figure S12 in Supporting Information S1). Note that stations NGZ and OTVZ lie in close proximity to Ngauruhoe volcano (Figure 1), which is also snow-covered in winter. Finally, we once again observe that stations closer to the summit show larger differences between different frequency bands (Figures 5a-5c).
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Averaged velocity changes between 1 and 2.5 Hz for seismic stations in Figure 3. Gray-shaded time periods represent times of greater snow load (from snow thickness data). (a) Velocity changes at DRZ and FWVZ stations. Snow radius shown as black-dashed line. (b) Velocity changes at MTVZ station, compared with 10-day smoothed pore pressure changes (dashed-black line) from Figure S10 in Supporting Information S1 (axis inverted). This smoothing is applied in order to match stack size of velocity changes.</figcaption></figure>

Following the increase in velocities close to the summit during winter, there appears to be a frequency-dependent difference in the time that velocities start to decrease. Of interest is that this coincides with the relative timing of the last two eruptions (in 2006 and 2007), represented as vertical dashed-lines in Figure 5. This is most apparent at DRZ (Figure 5a) and FWVZ (Figure 5c) stations, where velocities decrease first at the lower frequencies (1.0–1.5 Hz) from approximately late-September/early October while velocities remaining elevated at higher frequencies (1.5–3.0 Hz) through the rest of the year. At TRVZ station—similarly located on the volcano—this feature is less apparent (Figure 5b); however velocities (> 1 Hz) do show decrease during this same time period relative to early months (July through September). This similarly coincides with the relative timing of the last two eruptions. Thus, understanding this transition period, and the earlier decrease in velocity at some frequencies (i.e., 1–1.5 Hz), may provide clues toward understanding preferential timing of eruptions at the volcano.

## 3.2. Frequency to Depth Inversion

If Rayleigh waves are assumed to dominate the coda, the relationship between apparent velocity changes as a function of frequency ($\delta C/C(f)$), for any given mode, and the shear wave velocity perturbation as a function of depth ($\delta\beta/\beta(z)$) can be described by:

$$ \frac{\delta C}{C}(f) = \int_{0}^{\infty} K(f, z) \frac{\delta \beta}{\beta}(z) dz \qquad (2) $$

where $K(f, z)$ is a depth sensitivity kernel. This relationship assumes either a Equation 1 fixed Poissons ratio and density or Equation 2 fixed P-wave velocity and density (Haney & Tsai, 2017). It is then possible solve the inverse problem:

$$ d = G \cdot m \qquad (3) $$

where $d$ is the measured velocity changes at each frequency for a given day ($\delta C/C$), $G$ the matrix of sensitivity kernels ($K$) and $m$ the velocity changes with depth ($\delta\beta/\beta$). Thus we can estimate where the velocity changes may
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Velocity changes for all nine stations plotted against calendar day, stacking (averaging) the velocity of each year between 2005 and 2009. Black-dashed lines indicate day of the year that the 2006 and 2007 eruptions occurred.</figcaption></figure>

be occurring at depth. In particular, we focus on DRZ and FWVZ station, where the velocity is observed to decrease earlier between 1 and 1.5 Hz around the time of eruptive activity.

In reality, the coda may contain both surface and body waves (Obermann et al., 2013, 2016; Yuan et al., 2021), potentially complicating the location of changes with depth. Specifically, early parts of the coda are thought to be dominated by surface waves while body waves begin to dominate at later lag times (Obermann et al., 2013, 2016). This transition has been observed in wavefield simulations to occur at approximately 6–8 times the scattering mean free time (Obermann et al., 2016). While the mean free path length is unknown for Ruapehu, estimations of this parameter at 18 active volcanoes in Japan suggest a value of 1–3 km (at 1–4 Hz) could be expected (Hirose et al., 2022). Assuming an approximate Rayleigh wave velocity of 1 km/s at the surface—similar to that measured by Godfrey et al. (2017)—this transition would then be expected between 10 and 30 s lag time. Thus it is possible body waves sampling deeper portions of the volcano contribute to our measurements, where the maximum lag here is 25 s at 1 Hz and 15 s at 2 Hz. In any case, the assumption of Rayleigh waves can still provide a first-order estimate of the location of velocity changes. In the presence of body waves, this then represents a minimum estimate of the depth.

For the sensitivity kernel, the codes of Haney and Tsai (2017) are modified to output the kernel for a given velocity model. In this case, we use a four-layered model of shear-wave velocities proposed by Godfrey et al. (2017) based on ambient noise tomography (Table S4 in Supporting Information S1), with 20 m spacing. A constant $V_p/V_s$ ratio of 1.7 and density of 2000 kg/m³ are used, satisfying the conditions of Equation 2. The kernels corresponding to the Rayleigh wave fundamental mode of select frequencies are shown in Figure S13 in Supporting Information S1.
From Equation 3, we can resolve the apparent shear-wave velocity changes with depth by performing a simple weighted-damped least squares inversion of the form:

$$ \mathbf{m}_{est} = \mathbf{C}_m \mathbf{G}^T (\mathbf{G}\mathbf{C}_m \mathbf{G}^T + \mathbf{C}_d)^{-1} \mathbf{d}, \quad (4) $$

where $C_d$ is the data covariance matrix and $C_m$ is the model covariance matrix. In this study, we define the model covariance matrix following Haney and Tsai (2017) as:

$$ C_m(i,j) = \sigma_m^2 \exp\left(-\frac{|z_i - z_j|}{\lambda}\right), \quad (5) $$

where $\sigma_m = \gamma\bar{\sigma}_d$ is the standard deviation of the model parameters, with $\sigma_d$ the uncertainty in phase velocity measurement from the weighted linear regression (and $\bar{\sigma}_d$ the average) and $\gamma$ a user-defined tuning factor. $z_i$ and $z_j$ are the depth values at the i-th and j-th model points, and $\lambda$ is the correlation length. The data covariance matrix is defined as a diagonal matrix with measurement uncertainties on the diagonal ($C_d = \sigma_d I$).

Appropriate choices for $\lambda$ and $\gamma$ are identified using the L-curve criterion (Hansen & O'Leary, 1993). For simplicity we choose the same values for DRZ and FWVZ station, setting $\lambda = 20$ m and $\gamma = 60$. This choice lies in the region of the maximum curvature of the L-curve for both stations (Figure S14 in Supporting Information S1). Examples of the fit between observed and predicted apparent velocity changes are shown in Figure S15 in Supporting Information S1 for DRZ station and Figure S16 in Supporting Information S1 for FWVZ station.

Following the inversion, apparent velocity changes with depth are plotted for stations DRZ and FWVZ in Figure 6. Alongside these measurements we compute the restitution index (Vergely et al., 2010), calculated as the sum of the elements in each row of the resolution matrix (Backus & Gilbert, 1967), given as:

$$ R = C_m G^T (G C_m G^T + C_d)^{-1} G $$

Values closer to, or above, one indicate that the model estimates (velocity changes) are well resolved, while values less than 0.8 indicate that the model is not well constrained (Lesage et al., 2014; Vergely et al., 2010). At both stations, we see that the restitution index decays rapidly from approximately 1,500 m depth. This can be related to the minimum frequency used in the inversion (0.25 Hz), which reduces in sensitivity from approximately the same depth (Figure S13 in Supporting Information S1). Of interest are the high values between approximately 750–1,300 m. This likely reflects reduced measurement error from the weighted linear regression at the lowest used frequencies, which are more sensitive to this depth range. In any case, all values of the restitution index are greater than 0.8 above ~1,500 m, thus we consider the inversion capable of resolving model estimates above this depth.

The results of the inversion suggest that annual velocity changes measured above 1 Hz at summit stations are located within the upper few hundred meters (Figure 7). At DRZ station, the high velocity measured from July (>1 Hz) is resolved down to approximately 300 m depth (Figure 7c). Closer to the surface (<50–100 m), velocities are approximately zero or slightly negative. At FWVZ station, the velocity high measured during the same time period is resolved to 200–250 m depth (Figure 7d). Unlike at DRZ station, however, velocity changes are mostly positive closer to the surface (<50–100 m) during winter/spring months, albeit reduced in amplitude. Note, an earlier increase in velocity (from March) is observed at DRZ station (dashed-region in Figure 7a) that, after inverting for depth, would appear connected to the velocity increase from June (at 300–400 m depth). We consider this reflects a separate processes however, with their separation more apparent for FWVZ station (Figure 7d). This is discussed further in Section 4.2.3. Finally, at both stations, we observe an earlier decrease of velocities at 200–300 m depth—relative to slightly shallower depths of ~50–200 m—from approximately early September (~day 245). This reflects the earlier decreases in apparent velocities observed between 1.0 and 1.5 Hz (Figures 7a and 7b), coinciding with the relative timing of eruptive activity in 2006 and 2007. This feature is also clearly resolved in the predicted velocities based on the model solution (Figures 7e and 7f).
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Resolved velocity changes with depth for stations DRZ (a) and FWVZ (c) and the respective values of restitution index (b) and (d). Restitution index values less than 0.8 indicate that the model is not well constrained (Lesage et al., 2014; Vergely et al., 2010).</figcaption></figure>

# 4. Discussion

The discussion focuses on three key aspects. First, the validity of the results is considered in relation to the use of tremor as the dominant source above 1.0 Hz. We then consider the possible mechanisms that could drive the velocity changes observed. Finally, we discuss the possibility of seasonality in eruptive activity at Ruapehu in relation to observations made in this study.

## 4.1. Reliability of Velocity Changes Computed From Tremor

Many of the key observations presented come from the analysis of cross-correlation functions within the frequency band of volcanic tremor at Ruapehu. The application of cross-wavelet transform is intended to minimize the influence of varying seismic sources on our results. Specifically, measurements can be retained at frequencies where there is a strong enough similarity between individual cross-correlation functions and the reference function while unreliable measurements can be excluded. This exclusion is based on quality control criteria, where we define a minimum wavelet coherence required for individual delay time measurements. Similarly, we require that at least 25% of possible measurements within the coda are of non-zero weight (i.e., are above this threshold).

One approach to test the effectiveness of the wavelet coherence during quality control is to examine the influence of increasing the threshold. Thus, the minimum coherence threshold is increased to 0.8, and 0.9 (from 0.7) and the average yearly velocity changes re-computed for stations DRZ, FWVZ, and MTVZ (Figure S17 in Supporting Information S1). As the threshold is increased, we observe that progressively more data is excluded. This mostly occurs at the stations closer to the summit, where almost all data above 2.3 Hz is excluded at DRZ station at a
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Summary of inversion solution for upper 500 m at DRZ and FWVZ stations. Observed apparent velocity changes at DRZ (a) and FWVZ (b). Modeled velocity changes with depth at DRZ (c) and FWVZ (d). Predicted velocity changes, based on modeled solutions, at DRZ (e) and FWVZ (f). Dashed-region in (a) and (c) (DRZ) correspond to velocity changes considered distinct from velocity increase starting in June. The same feature at FWVZ is outlined in (b) and (d). Orange arrows in c and d highlight the decreasing velocity trend at depth (relative to surface) during time of year associated with eruptive activity (vertical black dashed-lines indicating day of year for 2006 and 2007 eruptions).</figcaption></figure>

threshold of 0.8 (Figure S17g in Supporting Information S1). A large amount of data above 1.5 Hz is also excluded at this station between July and November. These observations may relate to a change in scattering properties resulting from the addition of a compacting snow layer (Guillemot et al., 2021; Obermann et al., 2013). Of significance though is that increasing the threshold does not appear to influence the primary observations. For example, the trend of velocity changes decreasing first at lower frequencies—following the annual increase in winter—is still clearly present at both DRZ (Figure S17g in Supporting Information S1) and FWVZ station (Figure S17h in Supporting Information S1) when the threshold is increased to 0.8. Similarly, velocity changes at
MTVZ remain highly consistent (Figure S17i in Supporting Information S1). At the maximum tested threshold (0.9), there is significant data loss at all three stations about 1 Hz (Figures S17g–S17i in Supporting Information S1). However, the earlier decrease in velocity from September/October between 1.0 and 1.5 Hz—relative to higher frequencies (e.g., 1.5–2.5 Hz)—is still partially resolved. Thus, we do not believe our observations are sensitive to low quality measurements as a result of variable tremor characteristics.

Another factor to be considered is the influence that velocity changes computed during a single year could have on the yearly averages. We therefore examine the stability of the results when individual years are removed for the stations located on the volcano (Figure S18 in Supporting Information S1), thus including TRVZ station. Once again the main features of the velocity changes are well-preserved for stations DRZ and FWVZ no matter which year is excluded from the data. At TRVZ station, there is slightly more suggestion of a similar earlier velocity decrease between 1.0 and 1.5 Hz (following the winter maximum) when excluding either the year 2005 or 2006. However, the opposite could be argued if excluding either 2007 or 2008 from the average; thus, the average trend at this station is unclear with this data set. The influence of the eruptions themselves is also checked by removing the years 2006 and 2007 from the average (Figure S19 in Supporting Information S1). In this case the earlier decrease of velocities at lower frequencies is slightly less obvious at DRZ station, but still visible (Figure S19d in Supporting Information S1). At FWVZ station the decrease is still clearly visible when removing eruption years (Figure S19e in Supporting Information S1). Thus there is no evidence that anomalous velocity changes recorded during a single year are responsible for the main features we observe.

We also verify the consistency of the results when examining the yearly trends recovered using only a single component-pair (EN, EZ, or NZ), rather than using the averaged result (Figure S20 in Supporting Information S1). This choice also shows minimal influence on the results, with the same key features observed at both DRZ (Figures S20a–S20c in Supporting Information S1) and FWVZ (Figures S20d–S20f in Supporting Information S1) stations regardless of component-pair used. Interestingly, the greatest variability appears for different component-pairs at MTVZ station (Figures S20g–S20i in Supporting Information S1). While all component-pairs generally show higher velocities from January to July (and lower from July to December), variation between the different frequencies is observed. Velocity changes using the East-Vertical (EZ) component, for example, show an elevated velocity between approximately 0.5–1.0 Hz throughout the calendar year. Changes measured using the North-Vertical (NZ) components, however, show reduced velocities below 1 Hz throughout the year. In this case, the transition between the velocity high and low occurs at approximately 1 Hz; that is, the transition between tremor and non-tremor dominated cross-correlation functions. This could relate to differences in the paths sampled by the scattered waves originating from different sources, with tremor arriving at MTVZ station from the North-East. In contrast, the highly attenuative Taupo Volcanic Zone—similarly located to the North-East—has been observed to suppress Rayleigh waves originating from oceanic noise sources arriving from this direction (Behr et al., 2013). Varied directionality of different seismic sources, combined with non-uniformly distributed velocity changes, could thus explain this result.

Evidence of differences in recovered velocity changes depending on the dominant seismic source raises questions regarding the lowest frequency that the velocity high during winter is observed. At all three stations located on the volcano, the velocity increase during winter is observed down to 1 Hz (Figures 5a–5c). This similarly coincides with the lower limit of tremor (Figure 3a). Since it has been observed that tremor dominates above 1 Hz across the network (e.g., Figure 3g), and stations off the volcano show relatively consistent velocity changes across different frequencies, we do not believe that the change in the dominant seismic source has artificially defined this lower limit. However, given the consistency between the lower limit of tremor and increased velocity changes during winter, an influence cannot be conclusively ruled out. That said, minimal evidence is observed to suggest that the decreasing velocity at approximately 1.0–1.5 Hz from September/October is related to a change in tremor characteristics.

Combining these observations, we do not believe that the use of tremor as the dominant seismic source for cross-correlation functions above ~1 Hz has adversely impacted our results. Note, however, that this interpretation is only valid for the time-scales of features discussed (months rather than days). Shorter-term changes (days) could still be influenced by subtle variation in the character of tremor and would need to be scrutinized more closely before interpreting as real changes within the medium.
## 4.2. Mechanism of Velocity Changes

Differences between the trends of annual velocity changes at stations on and off the volcano suggest differences in the processes driving velocity changes. In this section we consider three processes that could influence seasonal seismic velocity changes. These are (a) Surface-loading due to snow mass, (b) Seasonal freeze-thaw cycles, and (c) Changes in groundwater level.

### 4.2.1. Surface-Loading Due To Snow Mass

The presence of a snowpack acts to apply an additional surface-load to the sub-surface (otherwise known as snow-loading). The increased pressure of such a load results in pore space reduction and crack closing in the medium (Silver et al., 2007). The effect of this is to increase seismic velocities (Nur, 1971). Thus, multiple studies have attributed winter increases in velocity to the presence of snow and subsequent snow-loading (Donaldson et al., 2019; Hotovec-Ellis et al., 2014; Makus et al., 2023; Wang et al., 2017).

We examine this further by applying a simple model of an axisymmetric surface load acting on a semi-infinite elastic medium. For a surface load $P(r) = \rho_c g h(r)$, the pressure perturbation with depth can be resolved following:

$$P(r, z) = -\frac{\sigma_{rr} + \sigma_{\theta\theta} + \sigma_{zz}}{3}$$

where $\sigma_{rr}$, $\sigma_{\theta\theta}$, and $\sigma_{zz}$ the diagonal components of the stress tensor (positive in tension) with cylindrical co-ordinates $r$, $\theta$, and $z$. These are computed using a MATLAB package (Le Contellec, 2024) based on the solution provided by Sneddon (1951), with $\sigma_{rr}$, $\sigma_{\theta\theta}$, and $\sigma_{zz}$ defined following:

$$\sigma_{rr} = \int_{0}^{\infty} (1 - z\xi) Z(\xi) e^{-\xi z} J_{0}(\xi r) d\xi + \frac{1}{r} \int_{0}^{\infty} Z(\xi) \left(z - \frac{1 - 2\nu}{\xi}\right) e^{-\xi z} J_{1}(\xi r) d\xi \quad (6)$$

$$\sigma_{zz} = z \int_{0}^{\infty} Z(\xi) e^{-\xi z} J_{0}(\xi r) d\xi + \int_{0}^{\infty} Z(\xi) e^{-\xi z} J_{0}(\xi r) d\xi \quad (7)$$

$$\sigma_{\theta\theta} = 2\nu \int_{0}^{\infty} Z(\xi) e^{-\xi z} J_{0}(\xi r) d\xi - \frac{1}{r} \int_{0}^{\infty} Z(\xi) \left(z - \frac{1 - 2\nu}{\xi}\right) e^{-\xi z} J_{1}(\xi r) d\xi \quad (8)$$

where

$$Z(\xi) = -\xi \int_{0}^{\infty} rp(r)J_{0}(\xi r)dr \quad (9)$$

and $J_0$ and $J_1$ are Bessel functions of order zero and one respectively, $\mu$ the shear modulus and $\lambda$ the Lamé coefficient (Le Contellec et al., 2024). The stress components are considered positive in tension, hence the negative sign in $P(r, z)$, and $z$ is orientated downwards.

Our preferred model for the snow-load is that of a sigmoid-shape, based on the logistic function, where the thickness of snow ($h$) is defined as a function of radius ($r$) following:

$$h(r) = h_{\max} \left( 1 - \frac{1}{1 + e^{-k(r/R - r_0)}} \right) \quad (10)$$

where $h_{max}$ is the maximum thickness of the snow, $R$ the radius of the load, $k$ the logistic growth rate (set to 10), or steepness of the curve, and $r_0$ the function midpoint. This is set to $r_0 = 0.5$, meaning that the midpoint of the function occurs at $r = R/2$. With this, the load is represented by an initial—approximately exponential— increase of $h$ beginning from $r = R$, slowing to a more linear increase around $r = R/2$ and finally increasing more slowly toward $h_{max}$ at $r = 0$ (Figure 8a). This choice reflects the observation that mountain regions often show an initial strongly positive elevation gradient of snow thickness which then levels out at a certain altitude,
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Model inputs for axisymmetric snow load (with sigmoid shape) applied at the surface and corresponding dimensionless pressures calculated for semi-infinite elastic medium. (a) Dimensionless load with sigmoid shape. (b) Dimensionless 2D pressure calculations plotted against radius and depth (both also dimensionless). (c) Surface load at different times of year using *R* inferred from MODIS satellite data and *h*<sub>max</sub> inferred from ski-station snow thickness measurements.</figcaption></figure>

even decreasing at the highest elevations (Grünewald & Lehning, 2011; Grünewald et al., 2014; Kirchner et al., 2014).

Modeling the temporal influence of the snow-load is achieved through varying the values of *R* and *h*<sub>max</sub> (Figure 8c). Values of *R* are derived using MODIS satellite data (Figure 1b). Values of *h*<sub>max</sub> are taken as the average snow thickness recorded by ski stations at Whakapapa and Turoa. No adjustment is made to these values due to their reduced elevation relative to the summit (~700 m difference) given uncertainties about the distribution of snow. Finally, the density of snow ρ<sub>c</sub> is set to 0.4 kg/m<sup>3</sup>—based on ski station measurements between 2005 and 2009 (average 0.42 kg/m<sup>3</sup>)—and Poissons ratio of the underlying rock ν set to 0.2. This is similar to values recorded by Mordensky et al. (2018) and consistent with andesitic rocks more generally (Heap et al., 2020). The final solution (Figure 8b) is 2D pressures (dimensionalized by ρ<sub>c</sub>gh<sub>max</sub>) in *z* — *r* space (dimensionalized by *R*).

Models suggest a pressure increase on the order of ~3–10 kPa from the snow-load (Figures 9a–9c). These are largest at axis of the load (r = 0), with the months of August and September associated with larger changes and greater depth penetration (Figure 9a). At all distances shown, however, the temporal pressure evolution follows a similar pattern. The month of July is associated with a more rapid velocity increase as the snow thickness increases, peaking mid-August. September then marks a transition to decreasing pressures and reduced depth penetration as the load reduces (Figure 9). This transition occurs slowly closer to the axis of the load (r < 2000 m,
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Pressure change calculations at different radius *r* between June 20 and October 21 for sigmoid-shaped snow-load and the estimated velocity changes from Equation 11. Note velocity changes shown down to 1,000 m only. (a) Pressure change at *r* = 0 m. (b) Pressure change at *r* = 2,000 m. (c) Pressure change at *r* = 4,000 m. (d) Velocity changes at *r* = 0 m. (e) Velocity changes at *r* = 2,000 m. (f) Velocity changes at *r* = 4,000 m. (g) Rate of pressure change at 1,000 and 2,000 m depth for *r* = 0 m. (h) Rate of pressure change at 1,000 and 2,000 m depth for *r* = 2,000 m. (i) Rate of pressure change at 1,000 and 2,000 m depth for *r* = 4,000 m.</figcaption></figure>

Figures 9a and 9b) relative to greater distances where the pressure decrease is more rapid (*r* > = 4000 m, Figure 9c).

The depth of pressure changes modeled as a result of the snow load (a few kilometers) are much deeper than velocity changes suggested by our inversion (upper few hundred meters). This can be explained by considering how the amplitude of shear-wave velocity changes as a function of depth. Based on relationships between seismic velocity and induced stresses derived by Tromp and Trampert (2018) and Fokker et al. (2021) show that seismic velocity changes (dβ/β) due to a pressure change (ignoring deviatoric stresses) can be approximated as:

$$ \frac{d\beta}{\beta} = \frac{\mu'}{2\mu} P \qquad (11) $$

where μ is the shear-wave modulus and μ' is pressure derivative (∂μ/∂P). Since μ is expected to increase with depth (and μ' expected to decrease) the magnitude of velocity change at the surface is expected to be larger than at depth for the same magnitude of pressure change. This idea is shown in Figures 9d–9f, with seismic velocity changes mostly localized within the upper few hundred meters following a transformation from pore pressure changes, consistent with inversion results (Figure 7). This was achieved using P- (Vp) and S- (Vs) wave velocities based on those from the original depth inversion, with full details found in Text S3 in Supporting Information S1.
The predicted velocity changes show some consistency with our observed results. For example, the changes are predominantly localized within the upper 400 m, similar to the results of the inversion (Figures 7c and 7d). Furthermore, the reduced penetration depth—and therefore earlier velocity decrease at approximately 200–300 m—is also recovered. However, there are also some important differences. For example, the observed velocities increase rapidly from June/July (Figures 7a and 7b), whereas the snow-loading model predicts only a modest increase during month of July (e.g., Figure 9d). There is also a strong mismatch is the magnitude of changes predicted using Equation 11 (0.1%). These are far less than those modeled from the inversion of measured velocity changes (>1.0%).

Several factors may be responsible for the mismatch between estimated and observed seismic velocity changes from a snow-load; for example, uncertainties in the near-surface values of the shear modulus and pressure derivative. Similarly, the distribution of the snow-load over the volcano is poorly constrained, where no modification is made to account for snow thicknesses measured at lower depths than the max elevation. Additionally, the model does not include topography, which could also be expected to influence the distribution of pressure. However, these uncertainties alone are not likely to explain an order of magnitude difference in modeled and observed velocity changes. One further consideration is the simplicity of modeling the medium as homogeneous. For example, any reduction in porosity due to the snow load would ordinarily act to increase the Elastic modulus and, therefore, seismic velocities (Heap et al., 2020; Yu et al., 2016). This process could expected to be most efficient at the near-surface, especially in the presence of porous volcanic rock. Thus, while the order of magnitude difference suggests that snow-loading alone may not be enough to explain the seasonal velocity changes observed on Ruapehu, we consider that the modeled changes represent a minimum estimate and, therefore, underestimate the true impact of the snow load on near-surface velocities.

### 4.2.2. Seasonal Freeze-Thaw Cycles

Frozen ground can have a significant effect on seismic velocity changes through strong increases of the shear modulus relative to unfrozen ground (Miao et al., 2019; Steinmann et al., 2021). While these changes are typically shallow (less than a few meters), large changes in the shear velocity can influence measurements of velocity at frequencies usually expected to be sensitive to deeper parts of the crust (Gassenmeier et al., 2014; Steinmann et al., 2021).

The similarity between seismic velocity changes expected due to freezing of the ground and that of snow-loading, or ground-water changes, raises questions on how to differentiate between them. For example, conflicting arguments have been provided for the process driving winter increases in velocity in eastern Hokkaido, Japan (Miao et al., 2019; Wang et al., 2017). Wang et al. (2017) originally interpret these increases as due to the combined effect of snow-loading and reduced infiltration of groundwater due to the snow cover. However, Miao et al. (2019) interpret these change as related to seasonal frozen soil, noting that eastern Hokkaido gets much less snow than western Hokkaido, which showed a weaker and opposite trend (Wang et al., 2017). They suggest that the absence of seasonal increases in western Hokkaido could be due to the thermal insulation effects of snow cover.

Snow cover has been widely found to have a strong influence on the thermal regime of the ground (e.g., Oliva et al., 2014; T. Zhang, 2005; W. Zhang et al., 2021). Very thin snow cover (<20cm) reflects incoming solar radiation due to its high albedo, and thus can act to cool the ground. However, as snow thickness increases it provides an insulating effect, limiting the transfer of heat between the soil and the atmosphere (T. Zhang, 2005). This decoupling occurs for thicknesses greater than approximately 40 cm (T. Zhang, 2005; W. Zhang et al., 2021). Similarly, it is often found that soil temperatures increase with increasing snow thickness (e.g., Goncharova et al., 2019; T. Zhang, 2005; W. Zhang et al., 2021). Thus, thick snow cover—as observed on Ruapehu—may prevent the development of frozen ground in winter months.

Considering the local geology can also help to distinguish between the different processes. Two of the three stations showing a winter velocity increase (FWVZ and TRVZ) are located on the flanks of the volcano, composed of lava flows and autobreccias mixed with glacial and pyroclastic deposits (Conway et al., 2016). Both stations are located on outcropping andesitic lava flows associated with the Whakapapa formation (post-glacial lava flows erupted since 15 ka) (Conway et al., 2016). It can be expected that a shear-wave velocity increase in andesitic rock following freezing would be, proportionally, far less than measured in soil as in Steinmann
et al. (2021). Thus, combined with the insulating effect of a snow-pack, we consider it unlikely that frozen ground would have as strong an influence on velocity changes in this study.

### 4.2.3. Groundwater Changes

Seasonal variations in precipitation, and subsequently the groundwater level, have regularly been observed to correlate with seismic velocity changes (e.g., Clements & Denolle, 2018; Mao et al., 2022; Obermann et al., 2014; Rivet et al., 2015; Sens-Schönfelder & Wegler, 2006; Wang et al., 2017). This is also observed off the volcano (at MTVZ station, for example), where seismic velocity changes are inversely correlated with pore pressure changes estimated from precipitation records (Figure 4b). This inverse relationship is thought to reflect changes in the forces between grain contacts as rock becomes saturated. Thus, an increase in pore pressure due to an increase in the groundwater table reduces the contact area between grains, thus promoting a seismic velocity decrease (Fokker et al., 2021; Grêt et al., 2006). These changes can be on the order of 0.1%–1.0%, similar to those observed at Ruapehu.

The presence of snow can have a significant influence on groundwater recharge (Earman et al., 2006; Nygren et al., 2020; Somers & McKenzie, 2020). At lower elevations, where minimal snow is present, rain may play a dominant role in recharging the groundwater system (Somers & McKenzie, 2020). At higher elevations, with lower temperatures, water infiltration is expected to decrease as precipitation falls instead as snow. This would act to lower groundwater levels, promoting a velocity increase in colder months through reduced pore pressures. Greater infiltration is then expected through snow-melt once temperatures increase above zero, resulting in a velocity decrease through increasing groundwater levels and pore pressures (Nygren et al., 2020). This phenomena may therefore also explain the opposite trends observed at stations on the volcano and those off the volcano (Figure 3).

The magnitude of changes expected from groundwater changes, and therefore pore pressure changes, can be described similarly to that of pressure changes in Equation 11. From Fokker et al. (2021), a change in pore pressure u via effective stress induces a shear-wave velocity change following:

$$ \frac{d\beta}{\beta} = -\frac{\mu'}{2\mu}u. \qquad (12) $$

Note that with the negative sign, a positive change in pore pressures results in a decrease in velocities. For the shear modulus μ and its derivative μ' used earlier in modeling the snow-loading, velocity increases of ~0.05% were predicted for a pressure change of 5–6 kPa within the upper few hundred meters. Thus, for velocity changes on the order of ±0.5% (i.e., a range of 1%) a pore pressure change of approximately 100 kPa would be required. Similarly, for ±1.0% change, pore pressure changes of 200 kPa would be required. This corresponds to groundwater changes equivalent to a level change of 10–20 m. This is not outside the realm of possibility for groundwater changes in mountainous regions (e.g., Ofterdinger et al., 2014); however, it is uncertain if changes of this magnitude are present on Ruapehu at the highest elevations. Similarly, there is a high uncertainty with regards to the location of the groundwater. Simulations suggest that the groundwater table may be located as much as a few hundred meters below the surface in mountainous regions (Forster & Smith, 1988; Somers & McKenzie, 2020; Somers et al., 2019). This does not take into account, however, that volcanoes can host perched aquifers within zones of high permeability (e.g., D'Ozouville et al., 2008; Nakano et al., 2022), thus could have groundwater located closer to the surface. The location or presence of such aquifers at Ruapehu is similarly poorly known.

One approach to distinguish between ground-water and snow-loading related changes might have been to study the shorter-term correlation between snow thicknesses and seismic velocity changes. If groundwater level changes due to reduced infiltration dominate, the rate of velocity increase should be less dependent on the amount of snow falling at a given moment. However, given the possible influence of seismic source changes, studying the short-term relationship between snow thickness and velocity changes is challenging. That said, perhaps the most convincing argument toward groundwater induced changes relates to the initial rapid increase in velocities in June/July (Figures 7a and 7b). These changes occur during a time period when snow accumulation is relatively modest relative to later (August/September) months. Thus, modeled velocities from snow-loading during this time period are small relative to later months (Figures 9d–9f). Thus, lower groundwater levels due to reduced
water infiltration represents a more likely explanation for the initial rapid increase in velocities. It is interesting to note, however, that the maximum velocity recorded in the winter of 2005—during which less snow accumulated —is lower than those of later winters (Figure 4a). At DRZ station, for example, maximum velocities recorded in winter 2005 are 0.2%–0.3% lower than in later winters. Similarly, they are 0.1%–0.2% lower at FWVZ station. Thus, while groundwater-induced changes still offer a more plausible explanation for a rapid velocity increase in the early winter months, differences in the maximum velocity recorded for each year may support an influence from snow-loading, as suggested in modeling of the snow-load (Section 4.2.1).

At stations at lower elevations, the seasonal trend of velocity changes appear to be consistent across a broad frequency range (Figure 5). At MTVZ, for example, the higher velocities beginning in January—and lasting until July—are observed at all frequencies between 0.5 and 4.0 Hz (Figures 3d and 5c). This is similar for all stations aside from those at higher elevations (DRZ, FWVZ, and TRVZ). We do, however, observe slightly higher velocities between approximately 0.5–1.5 Hz at these stations between January and July (Figures 5a, 5b, and 5f) prior to larger increases of velocities in winter. These higher velocities earlier in the year are clearly seen in the depth inversions for DRZ and FWVZ, at depths between approximately 200–600 m (Figures 6a, 6c, 7c, and 7d). This could be indicative of the same seasonal trend as observed at lower elevations. However, it is unclear why different behavior would be observed at different frequencies/depths on the volcano if ground-water changes were to dominate.

One consideration relates to the use of a dynamic lag time window (based on 20 cycles of a given period to define the end of the window). This could influence the sampling radius for different frequencies, with lower frequencies (longer periods) sensitive to a wider area than higher frequencies (shorter periods) due to longer travel-paths of scattered waves. Lower frequencies may then sample a wider portion of the medium down-slope of the station, also noting that the propagation direction of the initial tremor source is away from the crater lake/summit area. Thus scattered waves may preferentially sample down-slope regions of the volcano. To test the influence of the dynamic lag time window, the results of DRZ station are compared with those obtained using a static coda window (Figure S22 in Supporting Information S1). From this, elevated velocities earlier in the year—below approximately 1.25 Hz—are still recovered using static lag times of 5–15 s and 5–25 s (Figures S22b and S22c in Supporting Information S1). The earlier decrease between approximately 1.0–1.5 Hz from late-September/early October is also still present. This suggests that the choice of a dynamic coda window is not a significant factor in resolving these features. However, this does not take into account that seismic velocities of surface waves tend to increase with decreasing frequency. Thus, lower frequencies may still have wider areas of sensitivity even with a static coda window for all frequencies. Furthermore, laterally varying scattering properties can have a significant impact on the spatial sensitivity of coda waves (van Dinther et al., 2021). Thus, a strongly directional seismic source, combined with a non-uniform scattering medium, could also contribute to complex observations.

Without further understanding of the groundwater distribution on Ruapehu, it is difficult to determine whether groundwater changes can explain all the various observations. Similarly, limited understanding on the true distribution of snow on the volcano (in both space and time) makes it difficult to fully model this process. Both groundwater-level and snow-loading appear able to explain some of the observations, while also exhibiting limitations in their ability to account for certain phenomena. We consider the most likely scenario is a combination of both, with ground-water level changes more dominant (noting that velocity changes modeled from snow-loading are an order of magnitude smaller than observed). However, it is also likely that the influence of the snow-load is under-estimated when considering a homogeneous medium. Determining the true relative contribution of each process would require a more comprehensive understanding of spatial velocity changes, snow distribution, and sub-surface hydrology, on and around the volcano. Future application of geophysical surveys toward this purpose (e.g., Kereszturi et al., 2021; Miller et al., 2020; Peterson et al., 2021) would be beneficial toward this goal. Similarly, modeling a poroelastic medium would be beneficial to better constrain the true impact of a snow-load on near-surface seismic velocities.

### 4.3. Possible Links to Eruptive Activity

The study of seasonal changes at Ruapehu is motivated by evidence of seasonality in the eruptive record, with over 50% of eruptions—within the de-clustered catalog (Figure 2b)—occurring within the months of Spring (see Section 2.1). The most recent eruptions in 2006 and 2007 occurred in early October and late-September
respectively, just after seasonal velocities begin to decrease at depth (between 200 and 300 m). It is possible that, if a seasonal trigger exists, the same process driving this velocity decrease may be related to increased eruptive activity (directly or indirectly). However, as previously discussed (Section 4.2), it is unclear whether this specific feature is related to snow-loading/unloading or groundwater (melt-water) recharge. Thus we evaluate both mechanisms as potential triggers of increased spring-time activity.

The impact of melt-water on volcanic systems can be considered similar to the impact of rainfall. Multiple studies have linked increased volcanic activity to precipitation events (e.g., Elsworth et al., 2004; Gaete et al., 2020; Hort et al., 2003; Matthews et al., 2002; Uchida, 2024; Violette et al., 2001). Commonly invoked mechanisms are (a): percolation of meteoric water into the hydrothermal/magmatic system, where it interacts with hot rock and fluids (e.g., Gaete et al., 2020; Hort et al., 2003), and (b) pressure build-up due to infiltrating water blocking the escape of magmatic gases (e.g., Elsworth et al., 2004; Matthews et al., 2009). At Ruapehu, the main vents responsible for gas emissions are within the crater lake (Christenson, 2000; Werner et al., 2008). Thus additional groundwater is not expected to block pre-existing gas pathways. Furthermore, observations of snow coverage (and thickness) would suggest more significant snow-melt on the summit plateau occurs in late-spring/early Summer (December–February). This is consistent with observations at the neighboring volcano Ngauruhoe (500 m lower in elevation), with repeating seismicity—suggested as due to snow-melt infiltration (Jolly et al., 2012; Park et al., 2019)—detected between November and February in 2005–2009. At Ruapehu, however, these months show the lowest activity in the eruptive record (Figure 2b). Thus, a melt-water related trigger at the start of Spring would appear unlikely.

Seasonal snow-loading (or unloading) has been proposed as a mechanism driving eruptive activity at both Katla volcano in Iceland (Albino et al., 2010) and Villarrica volcano in Chile (De Groote, 2022). In the latter case, such a mechanism was similarly used to explain a high proportion of spring-time eruptions relative to other seasons (De Groote, 2022). In these cases it is considered that unloading can promote failure within magma chamber wall, promoting dyke initiation (Albino et al., 2010; De Groote, 2022). At Ruapehu, however, it seems unlikely that significant movements of magma toward shallow depths are involved since most of the spring-time activity consists of smaller, non-magmatic, events (Figure S3 in Supporting Information S1). This is reflected in the latest eruptions of 2006 and 2007, both considered to have resulted from rapid failure of a hydrothermal seal rather than any magmatic intrusion (Ardid et al., 2023; Christenson et al., 2010; Girona et al., 2018).

Alternatively, it is interesting to consider the relationship between externally induced stresses and magmatic degassing. For example, increases in gas emissions have been identified at a number of volcanoes associated with tidal patterns (e.g., Bredemeyer & Hansteen, 2014; Cigolini et al., 2009; Dinger et al., 2019). This is interpreted as due to magma displacement within the conduit induced by cycles of compression and extension (Bredemeyer & Hansteen, 2014; Girona et al., 2018; Michael & Christoffel, 1975). The resultant differential flow of magma can increase the bubble coalescence rate and thus the separation of gas from magma, enhancing volatile degassing (Bredemeyer & Hansteen, 2014; Dinger et al., 2019). The effectiveness of this process is suggested to be more significant at depths of several kilometers, with tide-induced bubble coalescence possibly negligible at shallow depths relative to classical transport mechanisms (Dinger et al., 2019). Furthermore, such a mechanism has been proposed previously toward explaining a correlation between fortnightly earth tides and eruptions at the neighboring volcano Ngauruhoe (Michael & Christoffel, 1975).

While the pressure changes due to a snow load (1–8 kPa, Figures 9a–9c) are generally larger than those of tidal stresses (0.1–10 kPa, (Dinger et al., 2019)), their timescales are vastly different. However, based on modeled pressures, the monthly rate of pressure change due to a snow load (at $z = 1$ and 2 km, Figures 9g–9i) is $\pm 1$ kPa or approximately 0.5 kPa per fortnight. This is within the range of fortnightly tidal stresses, suggesting that snow-loading/unloading may similarly be able to induce greater degassing through encouraging the differential flow of magma in response to an applied pressure gradient. Furthermore, the deeper propagation of these pressures (> 2 km, as suggested by our models), would imply some interaction with magma storage zones beneath the volcano (estimated between ~2 and 9 km depth (G. Kilgour et al., 2013)).

We therefore pose the question of whether snow-loading/unloading could enhance the separation of gas from deeper magma at Ruapehu (as illustrated in Figure 10). At the start of winter, the accumulation of snow leads to an increase in the surface load, applying a positive pressure change on the system and a corresponding velocity increase at the near-surface (Figure 10b). At the same time, increased pressure on the deeper magmatic system may promote increased degassing. As snow continues to accumulate through winter the pressure increase
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Conceptual model of possible unloading-induced degassing process related to loading/unloading of a snow-load. (a) Pre-winter (e.g., January–May): minimal snow-load present on the volcano. (b) Early winter (e.g., June–July): snow-load begins to accumulate, inducing a positive pressure change. (c) Late-winter (e.g., July–September): maximum snow-load present on the volcano. (d) Spring (e.g., September–November): snow-load begins to melt, with pressures decreasing in the subsurface.</figcaption></figure>

propagates to a greater depth, while velocities reach a maximum at the near surface (Figure 10c). In the presence of a hydrothermal seal, as inferred for previous eruptions (Ardid et al., 2023; Christenson et al., 2010), gas accumulates within the shallow hydrothermal system (approximately uppermost 400 m beneath summit plateau; Jones et al., 2008). Finally, as snow begins to melt in spring, the depth of penetration of the snow-load decreases and, in-turn, velocities increase at approximately 200–300 m depth (Figure 10d). Similarly, the reduction in pressure at depth may induce further degassing. Combined with decompression of the surrounding rock, and bubble expansion in gases accumulated within the shallow hydrothermal system, this subtle pressure modulation could act as the final trigger if the volcanic system is in a highly critical state.

Such a model also provides a potential explanation for the relative absence of eruptive activity in December through February (Figure 2b). If degassing processes are most significant during the loading/unloading process (e.g., from July to October), it follows that the likelihood of reaching a critical state within the hydrothermal system is elevated during this time period. However, if no eruption is triggered during removal of the load, it would imply that the system was not in (or did not reach) a critical state. Then, with no further mechanism for loading-induced degassing, the system remains in a state of low activity in the following months (e.g., Summer), with reduced likelihood of gas-driven eruptions.

It is clear that further investigation would be needed to confirm whether snow-loading/unloading plays a definitive role in triggering eruptions at Ruapehu. One approach could be to examine the yearly trends of measured gas fluxes during periods where the volcano is in a state of low activity (thus, measurements are not biased by stronger volcanic activity). However, measurements of gas flux at the volcano (Figure S24 in Supporting Information S1) at the volcano are relatively sparse. During the 13-year period between 2009 and 2022 (thus excluding 2006–2007 volcanic activity) the average sampling interval is approximately 60 days. Thus, while there is some evidence of increased degassing (SO₂) between July–November (relative to December–March), sample issues complicate the ability to interpret this further. Similarly, May and June are also associated with higher SO₂ fluxes between 2009 and 2022. Finally, any relationship may also be influenced by other factors for
example, the presence of a seal or whether new gas-rich magma is present. It is therefore difficult to assess the reliability of these trends without greater temporal resolution of gas-flux measurements. Toward this, the recent installation of four scanning DOAS (Differential Optical Absorption Spectrometry) sensors in 2023 should provide valuable measurements of SO₂ present in the volcanic plume for future comparison.

In the meantime it is useful to consider why eruptions do not occur more regularly if snow-loading/unloading is a potential trigger. No eruption has occurred within the last 15 years, with the last recorded activity in the catalog being volcanic earthquake and lake level increase in July 2009 due to rapid fluid release into the lake (Ardid et al., 2023; Scott, 2013). Examining the de-clustered catalog shows that many of the eruptions that occurred in spring months are clustered in time (Figure S23 in Supporting Information S1). For example, the time period from mid-1960s to early 1980s contains the majority of events. Similarly, we observe multiple instances where, just like in 2006 and 2007, consecutive spring eruptions occur with minimal activity in between; for example, in 1972–1973, 1980–1981, and 1997–1998. Thus the effectiveness of any external triggering mechanism likely requires that the volcanic system is already in an unstable state. With regards to the most recent eruptive activity, both the 2006 and 2007 phreatic eruptions were preceded by deep seismic swarms (5–15 km) in the 5–7 months before the eruption (G. N. Kilgour et al., 2014). These swarms are suggested to reflect the injection of new, gas-rich, magma into the magma storage region, possibly then making the system more prone to increased degassing and, thus, external triggering. Additionally, the presence of a seal—as previously discussed—may also be a necessary pre-requisite to allow pressure to buildup within the shallow hydrothermal system.

In summary, our analysis indicates a potential link between seasonal snow cycles and the likelihood of eruptive activity at Ruapehu. Similarly, the timing of the last two eruptions, which occurred during a period associated with a decrease in velocity at 200–300 m depth, could offer important insights into potential triggers. Investigating these changes further—both at Ruapehu and other volcanoes that similarly show increased activity following the removal of a seasonal snow-load (e.g., Villarica, De Groote, 2022)—would help clarify whether they are directly associated with the unloading mechanism or an alternative process. In any case, the observation that seasonality may exist in the eruptive record has important implications in terms of risk/hazard assessment. If this is the case, the time of year could be a valuable inclusion as a parameter in automated probabilistic forecasting model already in use at Ruapehu (Christophersen et al., 2022). Similarly, seismic velocity changes themselves could serve as a valuable monitoring parameter, provided a better understanding of the relative contribution of different external processes can be established.

**5. Conclusions**

Evidence of preferential timing in the eruption chronology at Ruapehu suggests the possibility of seasonal triggering of small-scale eruptions at the volcano. Spring months show the greatest activity, with over 50% of eruptions recorded in a declustered eruption catalog having occurred in the months of September, October, and November. The timing of these eruptions (in spring) suggests there may be some connection to the seasonal snow-cycles, with snow unloading previously suggested as a possible trigger for eruptive activity.

We study single-station seismic velocity changes over a 4 year period (2005–2009) that contains two spring-time, gas-driven, eruptions occurring approximately 1 year apart. To better handle the influence of variable tremor at higher frequencies (>1 Hz), the cross-wavelet transform approach is used to measure delay times. In addition to providing greater frequency resolution, the method allows data to be excluded in the time-frequency space when incoherent, while keeping data during time periods where tremor shows relative stability.

Substantial seasonal changes are observed that differ on and off the volcano. At stations on the volcano, velocities vary by ±0.5%, reaching a maximum during greater snow cover in winter months. At stations off the volcano, velocities closely follow modeled pore pressure changes associated with precipitation events. Of interest, we observe a pattern of an earlier velocity decrease at approximately 200–300 m depth (relative to the surface), which coincides with the timing of the last two eruptions. Ground-water changes are considered a likely mechanism responsible for seasonal velocity changes on the volcano, owing to reduced ground-water infiltration during winter. A contribution from snow-loading/unloading is also predicted from modeling. Further understanding of both the distribution of snow and groundwater on the volcano, and improved modeling of the loading process (e.g., poro-elastic modeling), would be beneficial toward determining the relative contribution of each process.
While ground-water related changes are suggested to have a stronger contribution to seasonal velocity changes, we consider snow-loading/unloading a more likely mechanism for greater spring-time activity. Modeling suggests snow-loading/unloading can induce stress changes at several kilometers depth at rates similar to those of tidal stresses (previously suggested to induce increased degassing). Thus we hypothesize that increased degassing associated with loading and unloading of the volcanic system may contribute to increased likelihood of gas-driven eruptions. However, long gaps without eruptions (many years), combined with evidence that spring-time eruptions cluster in time (consecutive years), would indicate that specific conditions are required for such a trigger to be effective. These interpretations clearly require further exploration, both at Ruapehu and other volcanoes that show similar preferential spring-time activity. Doing so could have important implications for risk/hazard assessment at volcanoes experiencing seasonal snow cover.

## Data Availability Statement

Seismic data are acquired through the data center GeoNet (GNS Science, 2021). The python package MSNoise (Lecocq et al., 2014) was used to process cross-correlation functions, and the python package Covseisnet (Covseisnet, 2020) used to compute the covariance matrix spectral width. The wavelet transform was performed using the python package PyCWT (Krieger & Freij, 2023) and velocity changes measured using python routines following methodology developed by Mao et al. (2020). Precipitation data were acquired via New Zealand's National Climate Database (NIWA, 2020). Snow data was provided by Ruapehu Alpine Lifts Ltd (RAL) on request, and MODIS satellite data (Hall et al., 2006) used to map snow coverage. Sensitivity kernels for depth inversion were prepared using modified codes of Haney and Tsai (2017). Modeling an axisymmetric surface load was achieved using the python package of Le Contellec (2024). Gas flux data are acquired through Geonet Field Time Series (FITS) database (GNS Science, 1954). Figures prepared using python package Matplotlib (Hunter, 2007).

## Acknowledgments

A. Yates was supported by funding from the Ministère de l'Enseignement supérieur, de la Recherche et de l'Innovation (MESRI). We also further acknowledge the New Zealand GeoNet programme and its sponsors EQC, GNS Science, LINZ, NEMA and MBIE for making the seismic data used in this study available. Finally, we thank the editor Dr. Satoshi Ide and reviewers (Manuela Köpfli and one anonymous reviewer) for their comments towards improving this manuscript.

## References

Albino, F., Pinel, V., & Sigmundsson, F. (2010). Influence of surface load variations on eruption likelihood: Application to two Icelandic subglacial volcanoes, Grímsvötn and Katla. *Geophysical Journal International, 181*(3), 1510–1524. https://doi.org/10.1111/j.1365-246X.2010.04603.x

Ardid, A., Dempsey, D., Caudron, C., Cronin, S. J., Miller, C. A., Melchor, I., et al. (2023). Using template matching to detect hidden fluid release episodes beneath crater lakes in Ruapehu, Copahue, and Kawah Ijen Volcanoes. *Journal of Geophysical Research: Solid Earth*, 128(10). https://doi.org/10.1029/2023JB026729

Backus, G. E., & Gilbert, J. F. (1967). Numerical applications of a formalism for geophysical inverse problems. *Geophysical Journal of the Royal Astronomical Society*, 13(1–3), 247–276. https://doi.org/10.1111/j.1365-246X.1967.tb02159.x

Ballmer, S., Wolfe, C. J., Okubo, P. G., Haney, M. M., & Thurber, C. H. (2013). Ambient seismic noise interferometry in Hawaii reveals long-range observability of volcanic tremor. *Geophysical Journal International*, 194(1), 512–523. https://doi.org/10.1093/gji/ggt112

Battaglia, J., Métaxian, J. P., & Garaebiti, E. (2012). Earthquake-volcano interaction imaged by coda wave interferometry. *Geophysical Research Letters*, 39(11), 4–7. https://doi.org/10.1029/2012GL052003

Behr, Y., Townend, J., Bowen, M., Carter, L., Gorman, R., Brooks, L., & Bannister, S. (2013). Source directionality of ambient seismic noise inferred from three-component beamforming. *Journal of Geophysical Research: Solid Earth*, 118(1), 240–248. https://doi.org/10.1029/2012JB009382

Bredemeyer, S., & Hansteen, T. H. (2014). Synchronous degassing patterns of the neighbouring volcanoes Llaima and Villarrica in south-central Chile: The influence of tidal forces. *International Journal of Earth Sciences*, 103(7), 1999–2012. https://doi.org/10.1007/s00531-014-1029-2

Brenguier, F., Shapiro, N. M., Campillo, M., Ferrazzini, V., Duputel, Z., Coutant, O., & Nercessian, A. (2008). Towards forecasting volcanic eruptions using seismic noise. *Nature Geoscience*, 1(2), 126–130. https://doi.org/10.1038/ngeo1014

Budi-Santoso, A., & Lesage, P. (2016). Velocity variations associated with the large 2010 eruption of Merapi volcano, Java, retrieved from seismic multiplets and ambient noise cross-correlation. *Geophysical Journal International*, 206(1), 221–240. https://doi.org/10.1093/gji/ggw145

Cannata, A., Cannavò, F., Montalto, P., Ercoli, M., Mancinelli, P., Pauselli, C., & Leto, G. (2017). Monitoring crustal changes at volcanoes by seismic noise interferometry: Mt. Etna case of study. *Journal of Volcanology and Geothermal Research*, 337, 165–174. https://doi.org/10.1016/j.jvolgeores.2017.03.023

Caricchi, L., Townsend, M., Rivalta, E., & Namiki, A. (2021). The build-up and triggers of volcanic eruptions. *Nature Reviews Earth & Environment*, 2(7), 458–476. https://doi.org/10.1038/s43017-021-00174-8

Caudron, C., Aoki, Y., Lecocq, T., De Plaen, R., Soubestre, J., Mordret, A., et al. (2022). Hidden pressurized fluids prior to the 2014 phreatic eruption at Mt Ontake. *Nature Communications*, 13(1), 1–9. https://doi.org/10.1038/s41467-022-32252-w

Caudron, C., Girona, T., Jolly, A., Christenson, B., Savage, M. K., Carniel, R., et al. (2021). A quest for unrest in multiparameter observations at Whakaari/White Island volcano, New Zealand 2007–2018. *Earth Planets and Space*, 73(1), 195. https://doi.org/10.1186/s40623-021-01506-0

Caudron, C., Lecocq, T., Syahbana, D. K., Mccausland, W., Watlet, A., Camelbeeck, T., et al. (2015). Stress and mass changes at a “wet” volcano: Example during the 2011–2012 volcanic unrest at Kawah Ijen volcano (Indonesia). *Journal of Geophysical Research: Solid Earth*, 120(7), 5117–5134. https://doi.org/10.1002/2014JB011590

Chouet, B. A., & Matoza, R. S. (2013). A multi-decadal view of seismic methods for detecting precursors of magma movement and eruption. *Journal of Volcanology and Geothermal Research*, 252, 108–175. https://doi.org/10.1016/j.jvolgeores.2012.11.013
Christenson, B. W. (2000). Geochemistry of fluids associated with the 1995–1996 eruption of Mt. Ruapehu, New Zealand: Signatures and processes in the magmatic-hydrothermal system. *Journal of Volcanology and Geothermal Research*, 97(1–4), 1–30. https://doi.org/10.1016/S0377-0273(99)00167-5

Christenson, B. W., Reyes, A. G., Young, R., Moebis, A., Sherburn, S., Cole-Baker, J., & Britten, K. (2010). Cyclic processes and factors leading to phreatic eruption events: Insights from the 25 September 2007 eruption through Ruapehu Crater Lake, New Zealand. *Journal of Volcanology and Geothermal Research*, 191(1–2), 15–32. https://doi.org/10.1016/j.jvolgeores.2010.01.008

Christophersen, A., Behr, Y., & Miller, C. (2022). Automated eruption forecasting at frequently active volcanoes using Bayesian networks learned from monitoring data and expert elicitation: Application to Mt Ruapehu, Aotearoa, New Zealand. *Frontiers in Earth Science*, 10, 1–21. https://doi.org/10.3389/feart.2022.905965

Cigolini, C., Poggi, P., Ripepe, M., Laiolo, M., Ciamberlini, C., Delle Donne, D., et al. (2009). Radon surveys and real-time monitoring at Stromboli volcano: Influence of soil temperature, atmospheric pressure and tidal forces on $^{222}$Rn degassing. *Journal of Volcanology and Geothermal Research*, 184(3–4), 381–388. https://doi.org/10.1016/j.jvolgeores.2009.04.019

Clarke, D., Zaccarelli, L., Shapiro, N. M., & Brenguier, F. (2011). Assessment of resolution and accuracy of the moving window cross spectral technique for monitoring crustal temporal variations using ambient seismic noise. *Geophysical Journal International*, 186(2), 867–882. https://doi.org/10.1111/j.1365-246X.2011.05074.x

Clements, T., & Denolle, M. A. (2018). Tracking groundwater levels using the ambient seismic field. *Geophysical Research Letters*, 45(13), 6459–6465. https://doi.org/10.1029/2018GL077706

Conway, C. E., Leonard, G. S., Townsend, D. B., Calvert, A. T., Wilson, C. J., Gamble, J. A., & Eaves, S. R. (2016). A high-resolution $^{40}$Ar/$^{39}$Ar lava chronology and edifice construction history for Ruapehu volcano, New Zealand. *Journal of Volcanology and Geothermal Research*, 327, 152–179. https://doi.org/10.1016/j.jvolgeores.2016.07.006

Conway, C. E., Pure, L. R., & Ishizuka, O. (2023). An assessment of potential causal links between deglaciation and eruption rates at arc volcanoes. *Frontiers in Earth Science*, 11(April). https://doi.org/10.3389/feart.2023.1082342

Covseisnet. (2020). [Software]. Retrieved from https://covseisnet.gricad-pages.univ-grenoble-alpes.fr/covseisnet/

De Groote, J. (2022). Intra-annual snow cover variations as a potential explanation for eruption seasonality at Villarrica volcano Chile Masters Thesis.

De Plaen, R. S. M., Lecocq, T., Caudron, C., Ferrazzini, V., & Francis, O. (2016). Single station monitoring of volcanoes using seismic ambient noise. *Geophysical Research Letters*, 43(16), 8511–8518. https://doi.org/10.1002/2016GL070078

Di Leo, G., & Sardanelli, F. (2020). Statistical significance: P value, 0.05 threshold, and applications to radiomics—Reasons for a conservative approach. *European Radiology Experimental*, 4(1), 18. https://doi.org/10.1186/s41747-020-0145-y

Dinger, F., Bredemeyer, S., Arellano, S., Bobrowski, N., Platt, U., & Wagner, T. (2019). On the link between Earth tides and volcanic degassing. *Solid Earth*, 10(3), 725–740. https://doi.org/10.5194/se-10-725-2019

Donaldson, C., Caudron, C., Green, R. G., Thelen, W. A., & White, R. S. (2017). Relative seismic velocity variations correlate with deformation at Kīlauea volcano. *Science Advances*, 3(6), e1700219. https://doi.org/10.1126/sciadv.1700219

Donaldson, C., Winder, T., Caudron, C., & White, R. S. (2019). Crustal seismic velocity responds to a magmatic intrusion and seasonal loading in Iceland's Northern Volcanic Zone. *Science Advances*, 5(11). https://doi.org/10.1126/sciadv.aax6642

D'Ozouville, N., Auken, E., Sorensen, K., Violette, S., de Marsily, G., Deffontaines, B., & Merlen, G. (2008). Extensive perched aquifer and structural implications revealed by 3D resistivity mapping in a Galapagos volcano. *Earth and Planetary Science Letters*, 269(3–4), 518–522. https://doi.org/10.1016/j.epsl.2008.03.011

Earman, S., Campbell, A. R., Phillips, F. M., & Newman, B. D. (2006). Isotopic exchange between snow and atmospheric water vapor: Estimation of the snowmelt component of groundwater recharge in the southwestern United States. *Journal of Geophysical Research: Atmospheres*, 111(9), 1–18. https://doi.org/10.1029/2005JD006470

Elsworth, D., Voight, B., Thompson, G., & Young, S. R. (2004). Thermal-hydrologic mechanism for rainfall-triggered collapse of lava domes. *Geology*, 32(11), 969–972. https://doi.org/10.1130/G20730.1

Farge, M. (1992). Wavelet transforms and their applications to turbulence. *Annual Review of Fluid Mechanics*, 24(1), 395–458. https://doi.org/10.1007/s10556-016-0140-6

Farquharson, J. I., & Amelung, F. (2020). Extreme rainfall triggered the 2018 rift eruption at Kīlauea Volcano. *Nature*, 580(7804), 491–495. https://doi.org/10.1038/s41586-020-2172-5

Fisher, R. A. (1935). The logic of inductive inference (with discussion). *Journal of the Royal Statistical Society*, 98(1), 39–82. https://doi.org/10.2307/2342435

Fokker, E., Ruigrok, E., Hawkins, R., & Trampert, J. (2021). Physics-based relationship for pore pressure and vertical stress monitoring using seismic velocity variations. *Remote Sensing*, 13(14), 1–20. https://doi.org/10.3390/rs13142684

Forster, C., & Smith, L. (1988). Groundwater flow systems in mountainous terrain: 2. Controlling factors. *Water Resources Research*, 24(7), 1011–1023. https://doi.org/10.1029/WR024i007p01011

Gaete, A., Walter, T. R., Bredemeyer, S., Zimmer, M., Kujawa, C., Franco Marin, L., et al. (2020). Processes culminating in the 2015 phreatic explosion at Lascar volcano, Chile, evidenced by multiparametric data. *Natural Hazards and Earth System Sciences*, 20(2), 377–397. https://doi.org/10.5194/nhess-20-377-2020

Gassenmeier, M., Sens-Schönfelder, C., Delatre, M., & Korn, M. (2014). Monitoring of environmental influences on seismic velocity at the geological storage site for CO$_2$ in Ketzin (Germany) with ambient seismic noise. *Geophysical Journal International*, 200(1), 524–533. https://doi.org/10.1093/gji/ggu413

Girona, T., Huber, C., & Caudron, C. (2018). Sensitivity to lunar cycles prior to the 2007 eruption of Ruapehu volcano. *Scientific Reports*, 8(1), 1476. https://doi.org/10.1038/s41598-018-19307-z

GNS Science. (1954). GeoNet Aotearoa New Zealand Manually Collected Volcano Data [Dataset]. https://doi.org/10.21420/PSP7-KW60?x=y

GNS Science. (2021). GeoNet Aotearoa New Zealand Seismic Digital Waveform Dataset [Dataset]. https://doi.org/10.21420/G19Y-9D40

Godfrey, H. J., Fry, B., & Savage, M. K. (2017). Shear-wave velocity structure of the Tongariro Volcanic Centre, New Zealand: Fast Rayleigh and slow Love waves indicate strong shallow anisotropy. *Journal of Volcanology and Geothermal Research*, 336, 33–50. https://doi.org/10.1016/j.jvolgeores.2017.01.019

Gómez-García, C., Brenguier, F., Boué, P., Shapiro, N. M., Droznin, D. V., Droznina, S. Y., et al. (2018). Retrieving robust noise-based seismic velocity changes from sparse data sets: Synthetic tests and application to Klyuchevskoy volcanic group (Kamchatka). *Geophysical Journal International*, 1218–1236. https://doi.org/10.1093/gji/ggy190

Goncharova, O. Y., Matyshak, G. V., Epstein, H. E., Sefilian, A. R., & Bobrik, A. A. (2019). Influence of snow cover on soil temperatures: Meso- and micro-scale topographic effects (a case study from the northern West Siberia discontinuous permafrost zone). *Catena*, 183, 104224. https://doi.org/10.1016/j.catena.2019.104224
Grêt, A., Snieder, R., & Scales, J. (2006). Time-lapse monitoring of rock properties with coda wave interferometry. *Journal of Geophysical Research, 111*(3), 1–11. https://doi.org/10.1029/2004JB003354

Grünewald, T., Bühler, Y., & Lehning, M. (2014). Elevation dependency of mountain snow depth. *The Cryosphere*, 8(6), 2381–2394. https://doi.org/10.5194/tc-8-2381-2014

Grünewald, T., & Lehning, M. (2011). Altitudinal dependency of snow amounts in two small alpine catchments: Can catchment-wide snow amounts be estimated via single snow or precipitation stations? *Annals of Glaciology*, 52(58), 153–158. https://doi.org/10.3189/172756411797252248

Guillemot, A., Van Herwijnen, A., Larose, E., Mayer, S., & Baillet, L. (2021). Effect of snowfall on changes in relative seismic velocity measured by ambient noise correlation. *The Cryosphere*, 15(12), 5805–5817. https://doi.org/10.5194/tc-15-5805-2021

Hall, D. K., Riggs, G. A., & Salomonson, V. V. (2006). MODIS/Terra Snow Cover 5-Min L2 Swath 500m. Version 5 [Dataset]. *Boulder, Colorado USA: NASA National Snow and Ice Data Center Distributed Active Archive Center*. https://doi.org/10.5067/ACYTYZB9BEOS

Haney, M. M., & Tsai, V. C. (2017). Perturbational and nonperturbational inversion of Rayleigh-wave velocities [Software]. *Geophysics*, 82(3), F15–F28. https://doi.org/10.1190/GEO2016-0397.1

Hansen, P. C., & O'Leary, D. P. (1993). The Use of the L-Curve in the Regularization of Discrete Ill-Posed problems. *SIAM Journal on Scientific Computing*, 14(6), 1487–1503. https://doi.org/10.1137/0914086

Heap, M. J., Villeneuve, M., Albino, F., Farquharson, J. I., Brothelande, E., Amelung, F., et al. (2020). Towards more realistic values of elastic moduli for volcano modelling. *Journal of Volcanology and Geothermal Research*, 390, 106684. https://doi.org/10.1016/j.jvolgeores.2019.106684

Hirose, T., Ueda, H., & Fujita, E. (2022). Scattering and intrinsic absorption parameters of Rayleigh waves at 18 active volcanoes in Japan inferred using seismic interferometry. *Bulletin of Volcanology*, 84(3), 34. https://doi.org/10.1007/s00445-022-01536-w

Hobiger, M., Wegler, U., Shiomi, K., & Nakahara, H. (2014). Single-station cross-correlation analysis of ambient seismic noise: Application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake. *Geophysical Journal International*, 198(1), 90–109. https://doi.org/10.1093/gji/ggu115

Hort, M., Seyfried, R., & Vöge, M. (2003). Radar Doppler velocimetry of volcanic eruptions: Theoretical considerations and quantitative documentation of changes in eruptive behaviour at Stromboli volcano, Italy. *Geophysical Journal International*, 154(2), 515–532. https://doi.org/10.1046/j.1365-246X.2003.01982.x

Hotovec-Ellis, A. J., Gomberg, J., Vidale, J. E., & Creager, K. (2014). A continuous record of intereruption velocity change at Mount St. Helens from coda wave interferometry. *Journal of Geophysical Research: Solid Earth*, 119(3), 2199–2214. https://doi.org/10.1002/2014JB011307

Hunter, J. D. (2007). Matplotlib: A 2D graphics environment [Software]. *Computing in Science and Engineering*, 9(3), 90–95. https://doi.org/10.1109/MCSE.2007.55

Hurst, A. W. (1998). Shallow seismicity beneath Ruapehu Crater Lake: Results of a 1994 seismometer deployment. *Bulletin of Volcanology*, 60(1), 1–9. https://doi.org/10.1007/s004450050212

Hurst, A. W., & Steven, S. (1993). Volcanic tremor at ruapehu: Characteristics and implications for the resonant source. *New Zealand Journal of Geology and Geophysics*, 36(4), 475–485. https://doi.org/10.1080/00288306.1993.9514593

Jaggar, T., Finch, R., & Emerson, O. (1924). The lava tide, seasonal tilt, and the volcanic cycle. *Monthly Weather Review*, 52(3), 142–145. https://doi.org/10.1175/1520-0493(1924)52<142:tltsta>2.0.co;2

Jellinek, A. M., Manga, M., & Saar, M. O. (2004). Did melting glaciers cause volcanic eruptions in eastern California? Probing the mechanics of dike formation. *Journal of Geophysical Research*, 109(9), 1–10. https://doi.org/10.1029/2004JB002978

Jolly, A. D., Neuberg, J., Jousset, P., & Sherburn, S. (2012). A new source process for evolving repetitious earthquakes at Ngauruhoe volcano, New Zealand. *Journal of Volcanology and Geothermal Research*, 215–216, 26–39. https://doi.org/10.1016/j.jvolgeores.2011.11.010

Jones, K. A., Ingham, M. R., & Bibby, H. M. (2008). The hydrothermal vent system of Mount Ruapehu, New Zealand—A high frequency MT survey of the summit plateau. *Journal of Volcanology and Geothermal Research*, 176(4), 591–600. https://doi.org/10.1016/j.jvolgeores.2008.05.006

Journeau, C., Shapiro, N. M., Seydoux, L., Soubestre, J., Ferrazzini, V., & Peltier, A. (2020). Detection, classification, and location of seismomovilcanic tremors with multi-component seismic data, example from the Piton de la Fournaise volcano (La Réunion, France). *Journal of Geophysical Research: Solid Earth*, 125, (8), 1–38. https://doi.org/10.1029/2019JB019333

Kereszturi, G., Schaefer, L., Mead, S., Miller, C., Procter, J., & Kennedy, B. (2021). Synthesis of hydrothermal alteration, rock mechanics and geophysical mapping to constrain failure and debris avalanche hazards at Mt. Ruapehu (New Zealand). *New Zealand Journal of Geology and Geophysics*, 64(2–3), 421–442. https://doi.org/10.1080/00288306.2021.1885048

Kilgour, G., Blundy, J., Cashman, K., & Mader, H. M. (2013). Small volume andesite magmas and melt-mush interactions at Ruapehu, New Zealand: Evidence from melt inclusions. *Contributions to Mineralogy and Petrology*, 166(2), 371–392. https://doi.org/10.1007/s00410-013-0880-7

Kilgour, G. N., Saunders, K. E., Blundy, J. D., Cashman, K. V., Scott, B. J., & Miller, C. A. (2014). Timescales of magmatic processes at Ruapehu volcano from diffusion chronometry and their comparison to monitoring data. *Journal of Volcanology and Geothermal Research*, 288, 62–75. https://doi.org/10.1016/j.jvolgeores.2014.09.010

Kirchner, P. B., Bales, R. C., Molotch, N. P., Flanagan, J., & Guo, Q. (2014). LiDAR measurement of seasonal snow accumulation along an elevation gradient in the southern Sierra Nevada, California. *Hydrology and Earth System Sciences*, 18(10), 4261–4275. https://doi.org/10.5194/hess-18-4261-2014

Krieger, S., & Freij, N. (2023). PyCWT: Wavelet spectral analysis in Python. V. 0.4.0-beta [Software]. https://pycwt.readthedocs.io/en/latest/Latter, J. H., & Balsillie, F. H. (1984). Volcano-seismic activity at ruapehu. In *New Zealand volcanological record* (pp. 31–37).

Lecocq, T., Caudron, C., & Brenguier, F. (2014). MSNoise, a Python package for monitoring seismic velocity changes using ambient seismic noise [Software]. *Seismological Research Letters*, 85(3), 715–726. https://doi.org/10.1785/0220130073

Le Contellec, A. (2024). Computing the stress field resulting from a surface topographic load on top of a semi-infinite medium in axisymmetric geometry [Software]. *Zenodo*. https://doi.org/10.5281/zenodo.8348775

Le Contellec, A., Michaut, C., Maccaferri, F., Pinel, V., Chambat, F., & Smrekar, S. (2024). Insights into Venus' crustal plateaus from dyke trajectories below craters. *Journal of Geophysical Research: Planets*, 129(5). https://doi.org/10.1029/2023JE008189

Leonard, G. S., Cole, R. P., Christenson, B. W., Conway, C. E., Cronin, S. J., Gamble, J. A., et al. (2021). Ruapehu and Tongariro stratovolcanoes: A review of current understanding. *New Zealand Journal of Geology and Geophysics*, 64(2–3), 389–420. https://doi.org/10.1080/00288306.2021.1909080

Lesage, P., Reyes-Dávila, G., & Arámbula-Mendoza, R. (2014). Large tectonic earthquakes induce sharp temporary decreases in seismic velocity in Volcán de Colima, Mexico. *Journal of Geophysical Research: Solid Earth*, 3782–3803. https://doi.org/10.1002/2015JB012608
Makus, P., Sens-sch, C., Illien, L., Walter, T., Yates, A., & Tilmann, F. (2023). Deciphering the whisper of volcanoes by combining conventional and Al-based approaches: Monitoring Kamchatka's Klyuchevskoy group with fluctuating noise fields. *Journal of Geophysical Research: Solid Earth*, 128(4). https://doi.org/10.1029/2022JB025738

Manville, V., Hodgson, K. A., Houghton, B. F., Keys, J. R., & White, J. D. (2000). Tephra, snow and water: Complex sedimentary responses at an active snow-capped stratovolcano, Ruapehu, New Zealand. *Bulletin of Volcanology*, 62(4–5), 278–293. https://doi.org/10.1007/s004450000096

Mao, S., Lecointre, A., van der Hilst, R. D., & Campillo, M. (2022). Space-time monitoring of groundwater fluctuations with passive seismic interferometry. *Nature Communications*, 13(1), 1–9. https://doi.org/10.1038/s41467-022-32194-3

Mao, S., Mordret, A., Campillo, M., Fang, H., & van der Hilst, R. D. (2020). On the measurement of seismic traveltime changes in the time-frequency domain with wavelet cross-spectrum analysis. *Geophysical Journal International*, 221(1), 550–568. https://doi.org/10.1093/gji/ggz495

Mardia, K. V. (1975). Directional statistics. *Journal of the Royal Statistical Society - Series B: Statistical Methodology*, 37(3), 349–371. https://doi.org/10.1111/j.2517-6161.1975.tb01550.x

Mason, B. G., Pyle, D. M., Dade, W. B., & Jupp, T. (2004). Seasonality of volcanic eruptions. *Journal of Geophysical Research*, 109(4), 1–12. https://doi.org/10.1029/2002JB002293

Matthews, A. J., Barclay, J., Carn, S., Thompson, G., Alexander, J., Herd, R., & Williams, C. (2002). Rainfall-induced volcanic activity on Montserrat. *Geophysical Research Letters*, 29(13), 22-1-22-4. https://doi.org/10.1029/2002GL014863

Matthews, A. J., Barclay, J., & Johnstone, J. E. (2009). The fast response of volcano-seismic activity to intense precipitation: Triggering of primary volcanic activity by rainfall at Soufrière Hills Volcano, Montserrat. *Journal of Volcanology and Geothermal Research*, 184(3-4), 405–415. https://doi.org/10.1016/j.jvolgeores.2009.05.010

McNutt, S. R., & Beavan, R. J. (1987). Eruptions of Pavlof Volcano and their possible modulation by ocean load and tectonic stresses. *Journal of Geophysical Research*, 92(B11), 11509–11523. https://doi.org/10.1029/jb092ib11p11509

Miao, Y., Shi, Y., Zhuang, H. Y., Wang, S. Y., Liu, H. B., & Yu, X. B. (2019). Influence of seasonal frozen soil on near-surface shear wave velocity in eastern Hokkaido, Japan. *Geophysical Research Letters*, 46(16), 9497–9508. https://doi.org/10.1029/2019GL082282

Michael, M. O., & Christoffel, D. A. (1975). Triggering of eruptions of Mt Ngauruhoe by fortnightly earth tide maxima, January 1972–June 1974. *New Zealand Journal of Geology and Geophysics*, 18(2), 273–277. https://doi.org/10.1080/00288306.1975.10418198

Miller, C. A., Schaefer, L. N., Kereszturi, G., & Fournier, D. (2020). Three-dimensional mapping of Mt. Ruapehu volcano, New Zealand, from aeromagnetic data inversion and hyperspectral imaging. *Journal of Geophysical Research: Solid Earth*, 125(2), 1–24. https://doi.org/10.1029/2019JB018247

Mordensky, S. P., Villeneuve, M. C., Kennedy, B. M., Heap, M. J., Gravley, D. M., Farquharson, J. I., & Reuschlé, T. (2018). Physical and mechanical property relationships of a shallow intrusion and volcanic host rock, Pinnacle Ridge, Mt. Ruapehu, New Zealand. *Journal of Volcanology and Geothermal Research*, 359, 1–20. https://doi.org/10.1016/j.jvolgeores.2018.05.020

Mordret, A., Courbis, R., Brenguier, F., Chmiel, M., Garambois, S., Mao, S., et al. (2020). Noise-based ballistic wave passive seismic monitoring—Part 2: Surface waves. *Geophysical Journal International*, 221(1), 692–705. https://doi.org/10.1093/gji/gga016

Mordret, A., Jolly, A. D., Duputel, Z., & Fournier, N. (2010). Monitoring of phreatic eruptions using interferometry on retrieved cross-correlation function from ambient seismic noise: Results from Mt. Ruapehu, New Zealand. *Journal of Volcanology and Geothermal Research*, 191(1–2), 46–59. https://doi.org/10.1016/j.jvolgeores.2010.01.010

Morlet, J., Arens, G., Fourgeau, E., & Giard, D. (1982a). Wave propagation and sampling theory—Part I. Complex signal and scattering in multilayered media. *Geophysics*, 47(2), 203–221. https://doi.org/10.1190/1.1441328

Morlet, J., Arens, G., Fourgeau, E., & Giard, D. (1982b). Wave propagation and sampling theory—Part II. Sampling theory and complex waves. *Geophysics*, 47(2), 222–236. https://doi.org/10.1190/1.1441329

Nakano, T., Yamada, Y., & Shin, K. C. (2022). Effects of snow and land modification on an andesite lava aquifer in Chokai volcano, northwestern Japan. *Journal of Hydrology*, 612(PB), 128191. https://doi.org/10.1016/j.jhydrol.2022.128191

Neuberg, J. (2000). External modulation of volcanic activity. *Geophysical Journal International*, 142(1), 232–240. https://doi.org/10.1046/j.1365-246X.2000.00161.x

Nishimura, T. (2017). Triggering of volcanic eruptions by large earthquakes. *Geophysical Research Letters*, 44(15), 7750–7756. https://doi.org/10.1002/2017GL074579

NIWA. (2020). CliFlo: NIWA's National Climate Database on the web [Dataset]. Retrieved from http://cliflo.niwa.co.nz/

Nur, A. (1971). Effects of stress on velocity anisotropy in rocks with cracks. *Journal of Geophysical Research*, 76(8), 2022–2034. https://doi.org/10.1029/JB076i008p02022

Nygren, M., Giese, M., Kløve, B., Haaf, E., Rossi, P. M., & Barthel, R. (2020). Changes in seasonality of groundwater level fluctuations in a temperate-cold climate transition zone. *Journal of Hydrology*, 8, 100062. https://doi.org/10.1016/j.hydroa.2020.100062

Obermann, A., Froment, B., Campillo, M., Larose, E., Planès, T., Valette, B., et al. (2014). Seismic noise correlations to image structural and mechanical changes associated with the Mw 7.9 2008 Wenchuan earthquake. *Journal of Geophysical Research: Solid Earth*, 119(4), 3155–3168. https://doi.org/10.1002/2013JB010932

Obermann, A., Planès, T., Hadzioannou, C., & Campillo, M. (2016). Lapse-time-dependent coda-wave depth sensitivity to local velocity perturbations in 3-D heterogeneous elastic media. *Geophysical Journal International*, 207(1), 59–66. https://doi.org/10.1093/gji/ggw264

Obermann, A., Planès, T., Larose, E., & Campillo, M. (2013). Imaging preeruptive and coeruptive structural and mechanical changes of a volcano with ambient seismic noise. *Journal of Geophysical Research: Solid Earth*, 118(12), 6285–6294. https://doi.org/10.1002/2013JB010399

Oferdinger, U. S., Renard, P., & Loew, S. (2014). Hydraulic subsurface measurements and hydrodynamic modelling as indicators for groundwater flow systems in the Rotondo granite, Central Alps (Switzerland). *Hydrological Processes*, 28(2), 255–278. https://doi.org/10.1002/hyp.9568

Oliva, M., Gómez Ortiz, A., Salvador, F., Salvà, M., Pereira, P., & Geraldes, M. (2014). Long-term soil temperature dynamics in the Sierra Nevada, Spain. *Geoderma*, 235–236, 170–181. https://doi.org/10.1016/j.geoderma.2014.07.012

Park, I., Jolly, A., Kim, K. Y., & Kennedy, B. (2019). Temporal variations of repeating low frequency volcanic earthquakes at Ngauruhoe Volcano, New Zealand. *Journal of Volcanology and Geothermal Research*, 373, 108–119. https://doi.org/10.1016/j.jvolgeores.2019.01.024

Peterson, D. E., Finn, C. A., & Bedrosian, P. A. (2021). Airborne geophysical imaging of weak zones on Iliamna volcano, Alaska: Implications for slope stability. *Journal of Geophysical Research: Solid Earth*, 126(3), 1–21. https://doi.org/10.1029/2020JB020807

Poupinet, G., Ellsworth, W. L., & Frechet, J. (1984). Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California. *Journal of Geophysical Research*, 89(B7), 5719–5731. https://doi.org/10.1029/JB089iB07p05719

Rivet, D., Brenguier, F., & Cappa, F. (2015). Improved detection of preeruptive seismic velocity drops at the Piton de La Fournaise volcano. *Geophysical Research Letters*, 42(15), 6332–6339. https://doi.org/10.1002/2015GL064835.Received
Scott, B. J. (2013). *A revised catalogue of Ruapehu volcano eruptive activity: 1830–2012*. GNS Science.

Sens-Schönfelder, C., & Wegler, U. (2006). Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia. *Geophysical Research Letters*, *33*(21), 1–5. https://doi.org/10.1029/2006GL027797

Seydoux, L., Shapiro, N. M., De Rosny, J., Brenguier, F., & Landès, M. (2016). Detecting seismic activity with a covariance matrix analysis of data recorded on seismic arrays. *Geophysical Journal International*, *204*(3), 1430–1442. https://doi.org/10.1093/gji/ggv531

Seydoux, L., Shapiro, N. M., de Rosny, J., & Landès, M. (2016). Spatial coherence of the seismic wavefield continuously recorded by the USArray. *Geophysical Research Letters*, *43*(18), 9644–9652. https://doi.org/10.1002/2016GL070320

Sherburn, S., Bryan, C. J., Hurst, A. W., Latter, J. H., & Scott, B. J. (1999). Seismicity of Ruapehu volcano, New Zealand, 1971–1996: A review. *Journal of Volcanology and Geothermal Research*, *88*(4), 255–278. https://doi.org/10.1016/S0377-0273(99)00014-1

Sigmundsson, F., Pinel, V., Lund, B., Albino, F., Pagli, C., Geirsson, H., & Sturkell, E. (2010). Climate effects on volcanism: Influence on magmatic systems of loading and unloading from ice mass variations, with examples from Iceland. *Philosophical Transactions of the Royal Society A: Mathematical, Physical & Engineering Sciences*, *368*(1919), 2519–2534. https://doi.org/10.1098/rsta.2010.0042

Sigvaldason, G. E., Annertz, K., & Nilsson, M. (1992). Effect of glacier loading/deloading on volcanism: Postglacial volcanic production rate of the Dyngjufjoll area, central Iceland. *Bulletin of Volcanology*, *54*(5), 393–404.

Silver, P. G., Daley, T. M., Niu, F., & Majer, E. L. (2007). Active source monitoring of cross-well seismic travel time for stress-induced changes. *Bulletin of the Seismological Society of America*, *97*(1 B), 281–293. https://doi.org/10.1785/0120060120

Sneddon, I. N. (1951). *Fourier transforms*. Mc. Grow-Hill Book Company, Inc.

Somers, L. D., & McKenzie, J. M. (2020). A review of groundwater in high mountain environments. *Wiley Interdisciplinary Reviews: Water*, *7*(6), 1–27. https://doi.org/10.1002/wat2.1475

Somers, L. D., McKenzie, J. M., Mark, B. G., Lagos, P., Ng, G. H. C., Wickert, A. D., et al. (2019). Groundwater Buffers decreasing glacier melt in an Andean Watershed—But not forever. *Geophysical Research Letters*, *46*(22), 13016–13026. https://doi.org/10.1029/2019GL084730

Soubestre, J., Shapiro, N. M., Seydoux, L., de Rosny, J., Droznin, D. V., Droznina, S. Y., et al. (2018). Network-based detection and classification of seismovolcanic tremors: Example from the Klyuchevskoy volcanic group in Kamchatka. *Journal of Geophysical Research: Solid Earth*, *123*(1), 564–582. https://doi.org/10.1002/2017JB014726

Steinmann, R., Hadzioioannou, C., & Larose, E. (2021). Effect of centimetric freezing of the near subsurface on Rayleigh and Love wave velocity in ambient seismic noise correlations. *Geophysical Journal International*, *224*(1), 626–636. https://doi.org/10.1093/gji/ggaa406

Talwani, P., Chen, L., & Gahalaut, K. (2007). Seismogenic permeability, ks. *Journal of Geophysical Research*, *112*(7), 1–18. https://doi.org/10.1029/2006JB004665

Torrence, C., & Compo, G. P. (1998). A practical guide to wavelet analysis. *Bulletin of the American Meteorological Society*, *79*(1), 61–78. https://doi.org/10.1175/1520-0477(1998)079(0061:APGTWA)2.0:CO:2

Tromp, J., & Trampert, J. (2018). Effects of induced stress on seismic forward modelling and inversion. *Geophysical Journal International*, *213*(2), 851–867. https://doi.org/10.1093/gji/ggy020

Uchida, N. (2024). Effects of extreme rainfall on phreatic eruptions: A case study of Mt. Ontake in Japan. *AIMS Geosciences*, *10*(2), 208–227. https://doi.org/10.3934/geosci.2024012

van Dinther, C., Margerin, L., & Campillo, M. (2021). Implications of laterally varying scattering properties for subsurface monitoring with coda wave sensitivity kernels: Application to volcanic and fault zone setting. *Journal of Geophysical Research: Solid Earth*, *126*(12), 1–21. https://doi.org/10.1029/2021JB022554

Vergely, J.-L., Valette, B., Lallement, R., & Raimond, S. (2010). Spatial distribution of interstellar dust in the Sun's vicinity. *Astronomy and Astrophysics*, *518*, A31. https://doi.org/10.1051/0004-6361/200913962

Violette, S., De Marsily, G., Carbonnel, J. P., Goblet, P., Ledoux, E., Tijani, S. M., & Vouille, G. (2001). Can rainfall trigger volcanic eruptions? A mechanical stress model of an active volcano: “Piton de la Fournaise”, Reunion island. *Terra Nova*, *13*(1), 18–24. https://doi.org/10.1046/j.1365-3121.2001.00297.x

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., et al. (2020). SciPy 1.0: Fundamental algorithms for scientific computing in Python. *Nature Methods*, *17*(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2

Wang, Q. Y., Brenguier, F., Campillo, M., Lecointre, A., Takeda, T., & Aoki, Y. (2017). Seasonal crustal seismic velocity changes throughout Japan. *Journal of Geophysical Research: Solid Earth*, *122*(10), 7987–8002. https://doi.org/10.1002/2017JB014307

Wasserstein, R. L., & Lazar, N. A. (2016). The ASA's Statement on p-values: Context, process, and purpose. *The American Statistician*, *70*(2), 129–133. https://doi.org/10.1080/00031305.2016.1154108

Werner, C., Hurst, T., Scott, B., Sherburn, S., Christenson, B. W., Britten, K., et al. (2008). Variability of passive gas emissions, seismicity, and deformation during crater lake growth at White Island Volcano, New Zealand, 2002–2006. *Journal of Geophysical Research*, *113*(B1). https://doi.org/10.1029/2007JB005094

Yates, A., Caudron, C., Lesage, P., Lecoocq, T., & Soubestre, J. (2023). Assessing similarity in continuous seismic cross-correlation functions using hierarchical clustering: Application to Ruapehu and Piton de la Fournaise volcanoes. *Geophysical Journal International*, *233*(1), 472–489. https://doi.org/10.1093/gji/ggac469

Yates, A. S., Savage, M. K., Jolly, A. D., Caudron, C., & Hamling, I. J. (2019). Volcanic, coseismic, and seasonal changes detected at White Island (Whakaari) Volcano, New Zealand, using seismic ambient noise. *Geophysical Research Letters*, *46*(1), 99–108. https://doi.org/10.1029/2018GL080580

Yokoyama, I. (1971). Volcanic eruptions triggered by tectonic earthquakes. *Geophysical bulletin of the Hokkaido University*, *25*, 129–139.

Yu, C., Ji, S., & Li, Q. (2016). Effects of porosity on seismic velocities, elastic moduli and Poisson's ratios of solid materials and rocks. *Journal of Rock Mechanics and Geotechnical Engineering*, *8*(1), 35–49. https://doi.org/10.1016/j.jrmge.2015.07.004

Yuan, C., Bryan, J., & Denolle, M. A. (2021). Numerical comparison of time-frequency-and wavelet-domain methods for coda wave interferometry. *Geophysical Journal International*, *226*(2), 828–846. https://doi.org/10.1093/gji/ggab140

Zhang, T. (2005). Influence of the seasonal snow cover on the ground thermal regime: An overview. *Reviews of Geophysics*, *43*(4). https://doi.org/10.1029/2004RG000157

Zhang, W., Shen, Y., Wang, X., Kang, S., Chen, A., Mao, W., & Zhong, X. (2021). Snow cover controls seasonally frozen ground regime on the southern edge of Altai Mountains. *Agricultural and Forest Meteorology*, *297*, 108271. https://doi.org/10.1016/j.agrformet.2020.108271
Supporting Information for "Seasonal snow cycles and their possible influence on seismic velocity changes and eruptive activity at Ruapehu volcano, New Zealand"

Alexander S., Yates<fnref n="1" />,<fnref n="2" />,<fnref n="3" />. Corentin, Caudron<fnref n="1" />,<fnref n="2" />,<fnref n="3" />. Aurélien, Mordret <fnref n="1" />,<fnref n="4" />.

Philippe, Lesage<fnref n="1" />. Virginie, Pinel<fnref n="1" />. Thomas, Lecocq<fnref n="5" />. Craig A., Miller<fnref n="6" />.

Oliver D., Lamb.<fnref n="6" />. Nicolas, Fournier.<fnref n="6" />.

Grenoble, France

<footnote n="1">¹Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTerre, 38000</footnote>

<footnote n="2">²Laboratoire G-Time, Department of Geosciences, Environment and Society, Université libre de Bruxelles, Belgium</footnote>

<footnote n="3">³WEL Research Institute, Avenue Pasteur, 6, 1300 Wavre, Belgium</footnote>

<footnote n="4">⁴Department of Geophysics and Sedimentary Basins, GEUS, Geological Survey of Denmark and Greenland, Copenhagen, Denmark</footnote>

<footnote n="5">⁵Seismology and Gravimetry Department, Royal Observatory of Belgium, Uccle, Belgium</footnote>

<footnote n="6">⁶GNS Science, Wairakei Research Centre, 114 Karetoto Road, RD4, Taupo 3384, New Zealand</footnote>
16

S1 Declustering of eruption catalog

17

18

19

20

Eruptive events at Mount Ruapehu are extracted from the revised catalogue of Scott (2013), which consists of 586 number of events dating back to the 1830s. Events are assigned a size between 0 and 5, following Table S1. We plot these by month in Figure 2a, excluding events that do not have a listed start month (all of which are pre-1900).
<table><caption>**Table S1.** Semi-quantitative estimate of Ruapehu eruptions based on observed effects, sourced from the revised catalogue of Ruapehu eruptive activity (Scott, 2013) (previously modified from Houghton et al. (1987), Sherburn et al. (1999), Wilson (2009))</caption><thead><tr><td>Eruption Scale</td><td>Observed Effects</td></tr></thead><tbody><tr><td>0</td><td>Crater Lake steaming, hotter than normal (i.e. above 30-35 °C) creating additional interest, but no observations (or confirmation) of activity in lake.</td></tr><tr><td>1</td><td>Small phreatic eruptions confined to Crater Lake.</td></tr><tr><td>2</td><td>Phreatic or phreatomagmatic eruption accompanied by surges; material deposited outside Crater Lake, but still confined to the crater basin. May produce larger flows/floods in Whangaehu Valley.</td></tr><tr><td>3</td><td>Deposition of material outside the crater basin; possible remobilisation/lahars in upper catchments and Whangaehu valley; OR small scale explosive eruptions/intermittent ash emission when no lake is present.</td></tr><tr><td>4</td><td>Material deposited well outside the crater basin onto the summit plateau and outer flanks. Lahars possible in several catchments; OR explosive ash eruptions when no lake is present producing columns up to 10000 ft.</td></tr><tr><td>5</td><td>Large scale explosive eruption displacing moderate volumes of the lake, lahars in all/most major valleys. The summit and slopes covered, with ashfall off the cone; OR explosive eruptions when no lake is present producing tall (> 10000 ft) eruption columns and ashfall off the cone.</td></tr></tbody></table>

It is important to ensure that the dataset has a good level of completeness before exploring patterns in the eruption history. Scott (2013) consider only data from 1940 onwards towards determining eruption recurrence rates. They also do not consider scale 0 and 1 events, as these events are severely under-represented. Further they observe that
25

26

27

28

29

30

31

the eruptive episodes of 1945 and 1995/96 have a significant impact on the eruptive frequency accounting for 65% of observations. Thus, they also remove all eruptions from 1945 and the 1995/1996 eruption sequence in their analysis. This was similarly done by Christophersen et al. (2022) when using a Bayesian network model to forecast eruptions at the volcano, where they exclude eruptions during the October 1945 to August 1946 and July 1995 to December 1997 magmatic episodes. Unlike Scott (2013), however, they do not exclude data pre-1940.

32

33

34

35

36

37

38

39

40

We similarly clean the catalog, including only size 2+ and excluding eruptive activity associated with the 1945 and 1995/96 eruptive episodes. Following Scott (2013) all eruptions during 1945 are removed, consisting of 133 events starting March 3 and ending December 31. Unlike Christophersen et al. (2022), we do not exclude any further data between January and August 1946, as this period contains only four events (with the first beginning April 17). Similarly, for the 1995/1996 eruptive sequence, we exclude activity using the period defined by Scott (2013) from January 1995 to September 1996, where the pause between the end of this period and the start of eruptive activity in 1997 is over one year. Thus, we prefer to include the 1997 activity in the dataset.

41

42

43

44

45

46

47

48

49

50

51

Size 2 events are also said to be under-represented in the catalogue (Scott, 2013). Given the small size of stresses associated with external triggering mechanisms (1-10 kpa) (Caricchi et al., 2021), it is preferable to include these events in the catalog however. Examining the cumulative event count for all sizes suggests that, while the 1940 cut-off considered by Scott (2013) appears to be a suitable towards ensuring good completeness for size 3 events, a later cut-off may be more suitable for size 2 events. For example, 80% of all size 2 events in the catalog are recorded after 1965. In contrast, the number of size 3 events increases steadily until 1945 (with 80% of events recorded after this year). Thus we prefer instead a cut-off year of 1960 in order to target a higher completeness of size 2 events. The updated catalog of size 2+ events, excluding pre-1960 events and the 1995/96 eruptive sequences is shown in Figure 2b.
<figure><img src="image_14.png" /></figure>

**Figure S1.** Reverse cumulative count of events of different size, after removing 1945 and 1995/96 eruptive activity following Scott (2013). Black-dashed line represented cut-off year considered by Scott (2013) towards analysing eruption recurrence rates. Red-dashed line represents cut-off year preferred here.

Finally, since we are interested in whether eruptions may be triggered by external processes at the volcano, it is necessary to de-cluster the catalog. Specifically, events that are causatively linked to a previous event should be excluded. Christophersen et al. (2022) use a de-clustering interval of 14 days when examining eruptions of size 3 or larger post-1830. This was chosen based on examining the change of slope when comparing the number of eruptions at different de-clustering intervals. They then remove all eruptions occurring within 14-days of a prior eruption.

We perform the same analysis with the remaining size 2+ events, after removing both events pre-1960 and the 1995/96 eruption sequence (Figure S2). From this, we find that any interval between 10- and 20- days appears a suitable choice, thus opt for the same 14-day interval as used by Christophersen et al. (2022) for size 3+ events. The final catalog after de-clustering, consisting of 49 eruptive events, is shown in Figure 2b.
<figure><img src="image_15.png" /></figure>

**Figure S2.** Number of eruptions plotted against de-clustering interval for size 2+ events in Scott (2013) eruption catalogue, after removing both events pre-1960 and the 1995/96 eruption sequence. The black-dashed line represents the 14-day interval, and intervals between 10 and 20 shaded grey.

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

**Figure S3.** Classification of de-clustered eruptions using descriptions from catalog of Scott (2013).

## S2 Statistical Analysis against non-uniformity

If there was no external triggering of eruptive activity at the volcano, it could be expected that the chance of an eruption on any given day is equal. Thus, we evaluate the statistical significance of the eruptive activity at Mt Ruapehu against a null hypothesis of uniformly distributed eruptions throughout the year. In Figure S4, the observed event counts are compared with those expected from a uniform distribution. We show this for events with sizes 2+ (Fig. S4a), events with sizes 2-3 (Fig. S4b), and also just size 2 events (Fig. S4c). This follows that we might expect that larger events are less
likely influenced by external triggering processes. Note, however, that there is less certainty regarding the completeness of the catalog of size 2 events.

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

**Figure S4.** De-clustered catalog of events at Mt Ruapehu from catalog of Scott (2013) (see Section S1), with only events post-1960, compared with expected events assuming a uniform distribution. (a) Size 2-5 events. (b) Size 2-3 events. (c) Size 2 events only.

We apply two statistical tests to assess the significance of variations in event count by month. The first of these is the Fisher's exact test, which calculates the exact probability of obtaining the observed distribution of frequencies under the null hypothesis of independence. This was preferred to the chi-squared test of independence, which is known to be unreliable for small sample sizes (Tate and Hyer, 1973). The second test we apply is the Rayleigh test, which is more suited to circular data. This test assesses whether the data exhibit significant clustering around a certain direction by computing the mean length of the data represented as vectors on a unit circle. Note therefore that the test assumes that data are unimodal.
The Fisher's exact test requires that we split the data into two independent categories. The probability of obtaining the observed distribution of these two variables against a null hypothesis of a uniform distribution is calculated. We separate the data in two ways; (1) comparing the frequency of eruptions in October to those outside of October, and (2) comparing the frequency of eruptions in September-November with those outside of this period. In the case of the former, we focus on October as this is the month with the largest number of eruptions. In the latter case, we group the time period from September to November, acknowledging that the assumption of independent variables may not be satisfied if a seasonal process is responsible for the higher than expected eruption counts during these three months (i.e. Fig. S4). An example of a contingency table used to calculate probabilities for this scenario is shown in Table S2.

<table><caption>**Table S2.** Contingency table used for Fisher's exact test, sizes 2+, grouping events between September-November. The expected count assumes that eruptions are equally likely on any day throughout the year</caption><thead><tr><th></th><th>Observed</th><th>Expected</th></tr></thead><tbody><tr><td>Sept-Nov Eruption</td><td>25.0</td><td>12.4</td></tr><tr><td>Dec-Aug Eruption</td><td>25.0</td><td>36.6</td></tr></tbody></table>

We compute p values following the Fisher's exact test using all events above size 2 and also just smaller events between size 2 and 3, with the final p values shown in Table S3. When using all events, we observed that the distribution is only significant (at 5% level) towards rejecting the null hypothesis when we group the eruptions between September and November (p=0.01). Using only smaller events, however, we achieve statistical significance regardless of whether the data are split into October versus non-October eruptions or Sept-Nov versus Dec-Aug eruptions (p=0.04 and p<0.01 respectively). Thus, the Fisher's exact test suggests there may be a significant difference in eruption rates during certain months. However, while the Fisher's exact test is well-suited to smaller sample sizes, the assumption of independent variables is not easily satisfied with temporal datasets. While the time period between September-November is associated with melting of snow, there is also melting of snow continuing into the new year. Thus, we
cannot rule out the possibility of a correlation between these months, which is why a circular representation of the dataset may be more appropriate (Figure 2c–d).

<table><caption>**Table S3.** Summary of p-value results for different tests and subsets of the data.</caption><thead><tr><th>Test</th><th>Size 2-5</th><th>Size 2-3</th><th>Size 2</th></tr></thead><tbody><tr><td>Fisher's exact test (October)</td><td>0.091</td><td>0.039</td><td>0.023</td></tr><tr><td>Fisher's exact test (September-November)</td><td>0.012</td><td>0.008</td><td>0.006</td></tr><tr><td>Rayleigh test</td><td>0.027</td><td>0.013</td><td>0.002</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Figure S5.** Normalized spectral width measurement between 2005–2015. Lower values indicate a more coherent seismic wavefield dominated by fewer seismic sources. Computed using stations DRZ, TRVZ, TUVZ, WNVZ. Eruption timing indicated by red-dashed lines.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

**Figure S6.** Example of delay-time measurements using wavelet transform for a cross-correlation function during coherent tremor (2007-03-01) at DRZ station (EN components). Dashed-lines correspond to minimum and maximum lag time for the coda window used in apparent velocity measurements. (a) Time-delay measurements in time-frequency space. (b) Coherence between the current and reference cross-correlation function. (c) Amplitude of wavelet cross-spectrum used in determining the weighting function for linear-regression. (d) Weighting function. (e) Velocity changes computed for DRZ (EN components), with the day corresponding to measurements in a–d indicated by the black arrow. Gray areas correspond to measurements that did not meet quality control criteria.
<figure><img src="image_20.png" /></figure>

**Figure S7.** Example of delay-time measurements using wavelet transform for a cross-correlation function during period of reduced tremor (2006-07-09) at DRZ station (EN components). Dashed-lines correspond to minimum and maximum lag time for the coda window used in apparent velocity measurements. (a) Time-delay measurements in time-frequency space. (b) Coherence between the current and reference cross-correlation function. (c) Amplitude of wavelet cross-spectrum used in determining the weighting function for linear-regression. (d) Weighting function. (e) Velocity changes computed for DRZ (EN components), with the day corresponding to measurements in a–d indicated by the black arrow. Gray areas correspond to measurements that did not meet quality control criteria.
<figure><img src="image_21.png" /></figure>

**Figure S8.** Example of delay-time measurements using wavelet transform for a cross-correlation function during period where tremor characteristics differ from longer-term average (2008-01-01) at DRZ station (EN components). Dashed-lines correspond to minimum and maximum lag time for the coda window used in apparent velocity measurements. (a) Time-delay measurements in time-frequency space. (b) Coherence between the current and reference cross-correlation function. (c) Amplitude of wavelet cross-spectrum used in determining the weighting function for linear-regression. (d) Weighting function. (e) Velocity changes computed for DRZ (EN components), with the day corresponding to measurements in a–d indicated by the black arrow. Gray areas correspond to measurements that did not meet quality control criteria.
<figure><img src="image_22.png" /><figcaption>**Figure S9.** Snow coverage and snow thickness datasets plotted against months of the year. (a) Snow coverage for individual years between 2005–2009. Vertical grey dashed-lines mark months displayed in (b) Average snow thicknesses recorded by Whakapapa and Turoa ski stations between 2005–2009. Black line corresponds to snow radius in (a).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure S10.** Rainfall measurements recorded at the Chateau station and corresponding mean pore pressure changes (black line) computed following Equation 1</figcaption></figure>
<figure><img src="image_24.png" /></figure>

**Figure S11.** Wavelet coherence between 10-day and reference cross-correlation functions for all nine stations plotted against calendar day, stacking (averaging) the coherence of each year between 2005–2009). Black-dashed lines indicate day of the year that the 2006 and 2007 eruptions occurred.
<figure><img src="image_25.png" /></figure>

**Figure S12.** Comparison of the timing of seasonal velocity high (observed at 2 Hz, black lines) against station elevation. Measurements correspond to times when apparent velocities are above 0%. Uncertainty in timing is indicated by red lines. Labels within each subplot correspond to those in Figure 5.

<table><thead><tr><td>Depth (m)</td><td>Vs (m/s)</td></tr></thead><tbody><tr><td>200</td><td>900</td></tr><tr><td>1500</td><td>1200</td></tr><tr><td>3000</td><td>1400</td></tr><tr><td>5000</td><td>1700</td></tr></tbody></table>

**Table S4.** Simplified velocity model of Godfrey et al. (2017) at Ruapehu.
<figure><img src="image_26.png" /></figure>

**Figure S13.** Rayleigh wave depth sensitivity kernels (fundamental mode) for selected frequencies, computed using codes of Haney and Tsai (2017) using velocity model of Godfrey et al. (2017).
<figure><img src="image_27.png" /></figure>

**Figure S14.** L-curve plots for various choices of regularization parameters $\lambda$ and $\gamma$ for DRZ and FWVZ stations. Different points are labeled with corresponding regularization parameters $(\lambda, \gamma)$. Red point corresponds to point of maximum positive curvature. (a) DRZ station, with points within dashed rectangle shown in (b). (c) FWVZ station, with points within dashed rectangle shown in (d).
<figure><img src="image_28.png" /><figcaption>**Figure S15.** Comparison between observed (black points) and predicted (red curve) velocity changes (based on resolved model parameters) for select days at DRZ station.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure S16.** Comparison between observed (black points) and predicted (red curve) velocity changes (based on resolved model parameters) for select days at FWVZ station.</figcaption></figure>
<figure><img src="image_30.png" /></figure>

**Figure S17.** Comparison of yearly velocity changes at DRZ, FWVZ, and MTVZ for different minimum values of yearly coherence. Gray areas correspond to measurements that did not meet quality control criteria. (a-c) minimum coherence = 0.7 for DRZ, FWVZ and MTVZ respectively. (d-f) minimum coherence = 0.8 for DRZ, FWVZ and MTVZ respectively. (g-i) minimum coherence = 0.9 for DRZ, FWVZ and MTVZ respectively
<figure><img src="image_31.png" /><figcaption>**Figure S18.** Comparison of yearly velocity changes at DRZ, FWVZ and TRVZ after removing individual years from average (a-c) year 2005 removed for DRZ, FWVZ, and TRVZ respectively. (d-f) year 2006 removed for DRZ, FWVZ, and TRVZ respectively. (g-i) year 2007 removed for DRZ, FWVZ, and TRVZ respectively. (j-l) year 2008 removed for DRZ, FWVZ, and TRVZ respectively.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>**Figure S19.** Comparison of year velocity changes when removing years of volcanic activity. (a-c) Average yearly velocity changes at DRZ, FWVZ, and TRVZ stations respectively during years 2005-2009. (d-f) Average yearly velocity changes at DRZ, FWVZ, TRVZ respectively after removing 2006 and 2007.</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>**Figure S20.** Yearly velocity changes using only a single component-pair at three stations. (a) DRZ, EN component-pair. (b) DRZ, EZ. (c) DRZ, NZ. (d) FWVZ, EN. (e) FWVZ, EZ. (f) FWVZ, NZ. (g) MTVZ, EN. (h) MTVZ, EZ. (i) MTVZ, NZ.</figcaption></figure>

## S3 Converting modeled pore pressures to velocity changes

The transformation from pore pressure changes to seismic velocity changes can be performed using P- ($V_p$) and S- ($V_s$) wave velocities based on the original depth inversion (Table S4). The density $\rho$ is then computed as $\rho = (V_p + 2.37)/2.8$ (Fig. S21b) (following the velocity-density relationship of Gebrande (1982) for volcanic rocks), the shear modulus $\mu$ as $\mu = V_s^2\rho$ (Fig. S21c), and finally the confining pressure $P(z)$ as $P(z) = \int_0^z \rho(z)gdz$ (Fig. S21d). From this the derivative of the shear modulus with respect to pressure ($\mu'$) can be resolved (Fig. S21e). To ensure smooth variations of these parameters, the velocity model of Godfrey et al. (2017) was fit with the analytical function of Lesage et al. (2018) with varied initial parameters (Fig. S21a). The obtained values of shear modulus vary between approximately 1 GPa at the surface and 3.5 GPa at 2000 m (corresponding to elastic modulus between approximately 2.5 GPa to 8.5 GPa, for isotropic and homogeneous materials). These are on the smaller side relative to laboratory mea-
surements of andesitic rock at Ruapehu (elastic modulus between approximately 1–40 GPa, (Mordensky et al., 2018; Heap et al., 2020)). Obtained values of the pressure derivative $\mu'$ (unitless) vary between approximately 250 at the surface and <100 beyond a few hundred meters depth. Few experiments have studied the influence of confining pressure on the elastic or shear modulus (Wu et al., 2019), thus this value is relatively unconstrained. However, approximating this value using data from multiple studies that show changes in the elastic, or shear, modulus against confining pressure suggests a range of values ranging between 1–1000 for both soil and rock samples (Sahaphol and Miura, 2005; Sawang-suriya et al., 2006; Arslan et al., 2008; Wang et al., 2009; Wu et al., 2019; Fokker et al., 2021, 2023). Thus, the values computed here are within the range of values that might be expected.

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

**Figure S21.** Model parameters used in converting pressure changes with depth to shear velocity changes. (a) P- and S- wave velocities ($V_p$ and $V_s$ respectively). These are fit to model of (Godfrey et al., 2017) from Table S4 using $V(z) = V_0[(z + a)^\alpha - a^\alpha + 1]$. [$V_{p0} = 1530$ m/s, $\alpha_p = 0.1$, $a_p = 110$] and [$V_{s0} = 900$ m/s, $\alpha_s = 0.1$, $a_s = 100$] are used respectively for smoothed P- and S- wave velocities. (b) Density profile computed from $V_p$ profile. (c) Shear modulus with depth, computed using $V_s$ and density profile. (d) confining pressure computed using density profile. (e) shear modulus derivative with respect to pressure.
## S4 Dynamic vs static window

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

**Figure S22.** Comparison between measurements made using (a) dynamic coda window and static lag times windows of (b) 5–15 seconds and (c) 5–25 seconds.

## S5 Eruptive activity discussion extra plots

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

**Figure S23.** Eruptions from de-clustered catalog through time. Orange bars correspond to spring eruptions, other eruptions shown as black bars. Red shaded region in 1995/96 corresponding to removed activity during de-clustering process (Section S1).
<figure><img src="image_37.png" /></figure>

**Figure S24.** (a) Violin plot of SO2 flux by month recorded at Ruapehu between 2009–2022. This shows vertically oriented kernel density plots representing the full distribution of the data, combined with box and whisker plots. Note, negative values in the kernel density plots are artifacts arising from flux values being close to zero. For days where multiple measurements were made, an average was taken. (b) Number of gas flux measurements used in a for each month.

References

Arslan, A. T., Koca, M. Y., Aydogmus, T., Klapperich, H., and Yılmaz, H. R. (2008). Correlation of unconfined compressive strength with young's modulus and poisson's ratio in gypsum from Sivas (Turkey). *Rock Mechanics and Rock Engineering*, 41(6):941–950.

Caricchi, L., Townsend, M., Rivalta, E., and Namiki, A. (2021). The build-up and triggers of volcanic eruptions. *Nature Reviews Earth and Environment*, 2(7):458–476.

Christophersen, A., Behr, Y., and Miller, C. (2022). Automated Eruption Forecasting at Frequently Active Volcanoes Using Bayesian Networks Learned From Monitoring Data and Expert Elicitation: Application to Mt Ruapehu, Aotearoa, New Zealand. *Frontiers in Earth Science*, 10(July):1–21.

Fokker, E., Ruigrok, E., Hawkins, R., and Trampert, J. (2021). Physics-based relationship for pore pressure and vertical stress monitoring using seismic velocity variations. *Remote Sensing*, 13(14):1–20.

Fokker, E., Ruigrok, E., Hawkins, R., and Trampert, J. (2023). 4D Physics-Based Pore Pressure Monitoring Using Passive Image Interferometry Geophysical Re-
search Letters. *Geophysical Research Letters*, 50.

Gebrande, H. (1982). 3.1. 2.3 Velocity-density relations: 3.1. 2 Elastic wave velocities and constants of elasticity at normal conditions. *Subvolume B*, pages 17–20.

Godfrey, H. J., Fry, B., and Savage, M. K. (2017). Shear-wave velocity structure of the Tongariro Volcanic Centre, New Zealand: Fast Rayleigh and slow Love waves indicate strong shallow anisotropy. *Journal of Volcanology and Geothermal Research*, 336:33–50.

Haney, M. M. and Tsai, V. C. (2017). Perturbational and nonperturbational inversion of Rayleigh-wave velocities. *Geophysics*, 82(3):F15–F28. [Software].

Heap, M. J., Villeneuve, M., Albino, F., Farquharson, J. I., Brothelande, E., Amelung, F., Got, J. L., and Baud, P. (2020). Towards more realistic values of elastic moduli for volcano modelling. *Journal of Volcanology and Geothermal Research*, 390:106684.

Houghton, B. F., Latter, J. H., and Hackett, W. R. (1987). Volcanic hazard assessment for Ruapehu composite volcano, taupo volcanic zone, New Zealand. *Bulletin of Volcanology*, 49(6):737–751.

Lesage, P., Heap, M. J., and Kushnir, A. (2018). A generic model for the shallow velocity structure of volcanoes. *Journal of Volcanology and Geothermal Research*, 356:114–126.

Mordensky, S. P., Villeneuve, M. C., Kennedy, B. M., Heap, M. J., Gravley, D. M., Farquharson, J. I., and Reuschlé, T. (2018). Physical and mechanical property relationships of a shallow intrusion and volcanic host rock, Pinnacle Ridge, Mt. Ruapehu, New Zealand. *Journal of Volcanology and Geothermal Research*, 359:1–20.

Sahaphol, T. and Miura, S. (2005). Shear moduli of volcanic soils. *Soil Dynamics and Earthquake Engineering*, 25(2):157–165.

Sawangsuriya, A., Bosscher, P. J., and Edil, T. B. (2006). Application of Soil Stiffness Gauge in Assessing Small-Strain Stiffness of Sand with Different Fabrics and Densities. *Geotechnical Testing Journal*, 29(3):207–216.

Scott, B. J. (2013). A revised catalogue of Ruapehu volcano eruptive activity: 1830-2012. *GNS Science*.

Sherburn, S., Bryan, C. J., Hurst, A. W., Latter, J. H., and Scott, B. J. (1999). Seismicity of Ruapehu volcano, New Zealand, 1971-1996: A review. *Journal of Vol-
canology and Geothermal Research, 88(4):255-278.

Tate, M. W. and Hyer, L. A. (1973). Inaccuracy of the X2 test of goodness of fit when expected frequencies are small. *Journal of the American Statistical Association*, 68(344):836-841.

Wang, X. S., Jiang, X. W., Wan, L., Song, G., and Xia, Q. (2009). Evaluation of depth-dependent porosity and bulk modulus of a shear using permeability-depth trends. *International Journal of Rock Mechanics and Mining Sciences*, 46(7):1175-1181.

Wilson, G. M. (2009). Monitoring at Ruapehu volcano ; can eruptions be predicted ?
*BSc Hon thesis, University of Canterbury*, page 79.

Wu, X., Jiang, Y., Guan, Z., and Gong, B. (2019). Influence of confining pressure-dependent Young's modulus on the convergence of underground excavation. *Tunnelling and Underground Space Technology*, 83(October 2018):135-144.