Multiphysics Modeling of Volcanic Unrest at Mt. Ruapehu (New Zealand)

Pre‐eruptive signals at the crater lake‐bearing Mt. Ruapehu (New Zealand) are either absent or hard to identify. Here, we report on geophysical anomalies arising from hydrothermal unrest (HTU) and magmatic unrest (MU) using multiphysics numerical modeling. Distinct spatio‐temporal anomalies are revealed when jointly solving for ground displacements and changes in gravitational and electrical potential fields for a set of subsurface disturbances including magma recharge and anomalous hydrothermal flow. Protracted hydrothermal injections induce measurable surface displacements (>0.5 cm) at Ruapehu's summit plateau, while magmatic pressurization (5–20 MPa) results in ground displacements below detection limits. Source density changes of 10 kg/m3 (MU simulations) and CO2 fluxes between 2,150 and 3,600 t/d (HTU simulations) induce resolvable residual gravity changes between +8 and −8 μGal at the plateau. Absolute self‐potential (SP) anomalies are predicted to vary between 0.3 and 2.5 mV for all unrest simulations and exceed the detection limit of conventional electric surveying. Parameter space exploration indicates that variations of up to 400% in the Biot‐Willis coefficient produce negligible differences in surface displacement in MU simulations, but strongly impact surface displacement in HTU simulations. Our interpretation of the findings is that monitoring of changes in SP and gravity should permit insights into MU at Ruapehu, while HTU is best characterized using ground displacements, residual gravity changes and SP anomalies. Our findings are useful to inform multiparameter monitoring strategies at Ruapehu and other volcanoes hosting crater lakes.


Multiphysics Modeling of Volcanic Unrest at Mt. Ruapehu (New Zealand)
2 of 23 Todesco & Berrino, 2005). Multi-parameter geophysical studies help to identify driving mechanisms and source properties behind volcanic unrest, especially when interpretations of field observations are combined with data modeling (e.g., Gottsmann, Flynn, & Hickey, 2020;Gottsmann et al., 2008;Hickey et al., 2016;Rinaldi et al., 2011;Wauthier et al., 2016;Zhan et al., 2021). Joint ground displacement and gravity change time series have, for example, been used at several volcanoes to interrogate enigmatic unrest processes Gottsmann, Biggs, et al., 2020;Zhan et al., 2019). While ground deformation monitoring (Sparks et al., 2012) is common at many restless volcanoes, monitoring of gravimetric and electrical potential field changes is scarce, despite joint inversion of multiphysics data sets such as from seismic and SP investigations providing useful information on the timing and evolution of different source mechanisms (Mahardika et al., 2012;Zlotnicki, 2015 and references therein). Here we present a suite of multiphysics models which jointly and simultaneously solve for ground displacements and gravitational and electrical potential field changes arising from MU and HTU processes. We test for the detectability of unrest signals and focus our study on Mt. Ruapehu in New Zealand, a volcano with a recent history of enigmatic unrest episodes which might herald renewal of eruptive activity.

Geological Background and Motivation
Mt. Ruapehu is a large stratovolcano of dominantly andesite composition and one of New Zealand's most active volcanoes. This volcano is North Island's highest peak at 2,797 m a.m.s.l. and it hosts three ski fields, which during winter months, hosts thousands of recreational users. Ruapehu is located in the Tongariro National Park (TNP), alongside two other active andesitic volcanoes (Ngauruhoe and Tongariro), the Tongariro Volcanic Centre forms the southwestern edge of the Taupo Volcanic Zone (TVZ; G. N. Kilgour et al., 2013;C. A. Miller et al., 2020;Rowlands et al., 2005). The TVZ is a NNE-trending rifted arc basin resulting from oblique, westward subduction of the Pacific Plate beneath the Australian Plate (Cole, 1990).
Volcanism at Ruapehu has been active for the past ∼250 ka (Gamble et al., 2003) with eruptive activity resulting from hydrothermal or magmatic perturbations, or a combination of both. HTU is thought to be provoked by the pulsating ascent of heat and magmatic fluids through the active HTS which feeds Ruapehu's acid crater lake (Te Wai -Moe in M ori; Hurst et al. (1991), Christenson and Wood (1993), Jones et al. (2008), Leonard et al. (2021)). Beneath the lake, geophysical (Ingham et al., 2009;Jones et al., 2008;Rowlands et al., 2005) and petrological (G. N. Kilgour et al., 2013) studies highlight a transcrustal mush zone within which distinct compositional magma batches are believed to reside (G. N. Kilgour et al., 2013;Nakagawa et al., 1999Nakagawa et al., , 2002. It has been proposed that MU might be triggered by the interaction of recharge from deeper reservoirs with remnant magmas stored in the crustal mush zone (Conway et al., 2020;Gamble et al., 1999;G. N. Kilgour et al., 2014;Nakagawa et al., 1999) with the potential to culminate in an eruption.
3 of 23 changes trigger ground displacements and gravity changes. HTU is simulated by injecting hot multi-phase and multi-component fluids (CO 2 , H 2 O) at the bottom of Ruapehu's HTS. Pore pressure and temperature changes trigger thermo-poroelastic responses in the HTS and edifice observable by ground displacements and changes in the gravitational and electrical potential fields. In this study, processes triggering MU and HTU are modeled in isolation to study geophysical fingerprints resulting from either unrest with the aim to identify key geophysical observables. The simulations solve different equations described next.

Strain-Induced Fluid Flow
For a single-phase fluid, strain-induced flow through a water-saturated porous rock can be described by Darcy's law: with v being the Darcy velocity, κ being the permeability of the porous rock, η f being the fluid's viscosity, p f being the pore pressure and p 0 the initial pore pressure distribution. Driving forces for fluid flow are temporal strain changes ( ) in the subsurface (Biot, 1962;Wang, 2000): where q is the mass source/sink, ρ f is the fluid density, and ϵ vol is the volumetric strain. S PE is the poroelastic with Φ the porosity of the medium, χ f the fluid's compressibility and K the bulk modulus. Equations 1 and 2 denote the solid-to-fluid coupling and are solved for MU simulations solely.

Hydrothermal Model
The simulation of HTU is based on the flow of fluid and heat in a porous medium by solving a set of mass and energy balance equations which are described as follows (see Pruess et al. (1999); Xu et al. (2004) for further reading): with F being the flux, q being the source/sink term and Q being the accumulation term for m mass components (H 2 O and CO 2 , α = 1, 2 hence m = 2) and the energy equation (α = N). Accumulation term (Q α ) and fluid fluxes for the mass balance equation are defined as follows: with the subscript β = l or g characterizing the liquid and gas phase, respectively, with the permeability κ β , the density ρ β , the saturation S β and the mass fraction of component m in phase β. Mass fluxes F α can be calculated by using the extended Darcy's law for multi-component and multi-phase fluid flow: with the Darcy's velocity v β in phase β, the relative permeability k rβ and the gravitational acceleration vector g. All other parameters are equivalent to the single-phase fluid, but accounting for different phases (β = l or g).
For the energy equation, the accumulation term (Q N ) and heat flux (F N ) are defined as: where e β is the specific internal energy and h β the specific enthalpy in phase β, T is the temperature and ρ, cp, and λ are the density, heat capacity and thermal conductivity of the porous medium, respectively. 4 of 23

Observables
Our study focuses on ground displacement, gravity changes and SP signals as described below. We test for the detectability of modeled signals, with common GNSS surveys at Ruapehu resolving ground displacements of 0.5 cm horizontally to 1 cm vertically (Mordret et al., 2010). Gravity changes at the 5 μGal level are resolvable with carefully executed standard survey protocols (Battaglia et al., 2008) and the detectability of SP observations ranges between a few and 100 microvolts (μV), with most field equipment resolving 0.1 mV (Crespy et al., 2008;Revil & Jardani, 2013;Zlotnicki, 2015).

Ground Displacement
In areas where rocks are in quasi-static equilibrium (deformation processes occur slowly), displacement resulting from thermo-poroelastic responses can be derived from Hooke's law coupled with pressure and temperature effects (Fung, 1965;McTigue, 1986;Rice & Cleary, 1976): with the stress σ b and strain ε tensor, the displacement vector u and the identity matrix I. The elasticity matrix C = C(E,ν) and bulk modulus ( = 3(1−2 ) ) are represented by the Young's modulus (E) and the Poisson's ratio (ν). The first term on the right hand side of Equation 7 represents Hooke's law of linear elasticity, while the second and third term account for stress and strain variations resulting from temperature (ΔT) and pore pressure (Δp f ) changes, respectively. Key parameters for thermo-poroelastic response are the volumetric thermal expansion coefficient α T and Biot-Willis coefficient α BW . For MU simulations, we fully couple poroelastic responses with stress and strain changes affecting fluid flow and vice versa, while HTU simulations represent a one-way-coupling, where temperature and pressure changes control deformation process but not vice versa.
In volcanic areas where magmatic reservoirs heat surrounding rocks, viscoelastic behavior most appropriately characterizes time-dependent deformation processes (Del Negro et al., 2009). Therefore, we invoke a temperature-dependent viscoelastic rheology (see Text S1 and Figure S1 in Supporting Information S1) of the crust in the MU model by solving stress-strain relations using a Standard Linear Solid (SLS) parameterization (Del Negro et al., 2009;Hickey & Gottsmann, 2014;Hickey et al., 2016), which is most representative for crustal material (Head et al., 2019(Head et al., , 2021. The SLS model consists of an elastic branch controlled by the shear modulus G and a viscoelastic branch characterized by the relaxation time τ 0 : with the shear viscosity η r . Both branches are split equally using the fractional components (μ 1 = μ 0 = 0.5) of G. The shear viscosity is derived using the Arrhenius approximation: where A is the Dorn parameter (A = 10 9 Pa s), H is the activation energy (H = 120,000 J mol −1 ), R is the gas constant (R = 8.314 J mol −1 K −1 ) and T is the temperature. In our parameterization, near-elastic behavior of rocks (over timescales relevant for the study) occur in volumes of low temperature such as the edifice.

Gravity Changes
Gravity changes at volcanoes are attributed to subsurface density changes Δρ(x, y, z) resulting from the redistribution of hydrothermal fluids (e.g., Rinaldi et al., 2011;Todesco et al., 2010) or magma, and shifting of density boundaries by concurrent ground deformation. (Bonafede & Mazzanti, 1998; 5 of 23 Currenti, 2014;Gottsmann, Biggs, et al., 2020). Gravity changes Δg are calculated by solving the Poisson's equation for the gravitational potential ϕ g (Cai & Wang, 2005): where G is the gravitational constant. By imposing Dirichlet boundary conditions of zero at infinity the mathematical problem is closed.
Subsurface density changes for MU simulations consist of three source terms and can be calculated as follows (Bonafede & Mazzanti, 1998;Currenti, 2014;Zhang et al., 2004): with the density of the medium ρ r and the source density change Δρ m . The first term on the right-hand side results from the displacements of subsurface density boundaries. The second term quantifies density variations in the transcrustal magma reservoir due to influx of new mass, controlled by the contraction of resident magma and the reservoir expansion. The third term accounts for the compressibility of the surrounding rock (Bonafede & Mazzanti, 1998).
Density variations from fluid redistribution in HTU simulations are calculated with respect to the initial fluid density distribution (ρ 0 ; ; Currenti and Napoli (2017)): for each time step (k), using the fluid density and saturation in phase β.
We derive residual gravity changes from ( ) = − − , where − is the free-air effect, with the vertical displacement w and the theoretical Free-Air gradient (−308.6 μGal/m).

Self-Potential
Self-potential (SP) anomalies in porous media arise from the drag of excess charge with the fluid flow (electrokinetic processes, e.g., Revil and Florsch, 2010;Revil et al., 2012). Here, we couple SP signals to pore pressure changes in response to strain-induced fluid flow (MU) or the injection of a hot multi-phase and multi-component fluid (HTU). The total current density (j) resulting from electrokinetic processes is calculated as follows (see Bolève et al., 2011;Revil, Pezard, and Glover, 1999;Revil, Schwaeger, et al., 1999;Sill, 1983; for details): where σ is the electrical conductivity of the medium, φ the electrical potential, p f the pore pressure and L SP the streaming current coupling coefficient. The latter is related to the streaming-potential coupling coefficient (C SP ) via L SP = −C SP σ, whereby C SP is a key parameter to quantify hydro-electric mechanisms. Applying the continuity equation for electrical charge (∇ ⋅ j = 0) to Equation 13 yields the Poisson's equation: where ℑ is the volumetric current source density defined as ℑ = −∇ ⋅ ( ∇ ) . An electrical reference potential is set to zero at an arbitrary point, as the electrical potential is a relative measure.

Model Implementation
We develop a suite of 2D axisymmetric forward models to simultaneously solve for ground displacement, SP and gravity changes at Mt. Ruapehu by MU and HTU. All numerical models incorporate topography as well as subsurface mechanical and hydro-electric heterogeneity (see Table 2). 6 of 23

Magmatic Unrest Model
We simulate MU using the commercial Finite-Element Analysis package COMSOL Multiphysics (Version 5.3). Figure 1 shows the model geometry and domain size with a radial (r) and vertical extent (z) of 50 and 75 km, respectively. The crust (z < 0 km) and the edifice (z ≥ 0 km) make up the main domains. The edifice above 1.5 km a.m.s.l. (Figure 1b) is divided into three sub-domains representing the HTS (r = 0-50 m), the transition zone (TZ, r = 50-150 m) and the edifice (r > 150 m), respectively.
A vertically extended mush zone is proposed to be located between 2 and 9 km depth below mean sea level (Ingham et al., 2009;Jolly et al., 2010;G. N. Kilgour et al., 2013;Rowlands et al., 2005). We represent this mush zone as a prolate ellipsoidal domain embedded in a viscoelastic crust. The semi-axes of the ellipsoid are derived to match reservoir volume estimates from previous eruptive volumes following Browning et al. (2015) (see Text S2 in Supporting Information S1). See Table S2 in Supporting Information S1 for the reservoir geometry of a maximum dense rock equivalent eruptive volume of 3 ⋅ 10 7 m 3 (G. Kilgour et al., 2010) with a reference tensile strength of 10 MPa. The injection of new magma into a transcrustal reservoir can trigger pressurization Figure 1. Illustration of the 2D asymmetric model setup for magmatic unrest (MU; upper panels) and hydrothermal unrest (HTU; lower panels) simulations. (a) The mush zone (red ellipsoid) is embedded in a viscoelastic crust located at a center depth of z = −5.5 km on the symmetry axis. Boundary conditions for solid mechanics are also shown. The edifice above z ≥ 1.5 km (b) is divided into the hydrothermal system, the transition zone and the edifice, in which poroelastic and electric processes are simulated. A no-flow and electric insulation boundary surrounds these domains, while internal boundaries are treated as continuous. An electrical reference potential (V Ref = 0 (V)) is applied at r = 5.5 km (yellow circle). Ruapehu's crater lake is shown for representation but is excluded in the numerical models. The lower panels show the meshes for HTU simulations, with the TOUGHREACT model being confined to the edifice (c). The electro-mechanical model dimension (d) is extended radially and vertically, with the red line representing Ruapehu's topography. Model domains for HTU simulations are equivalent to the MU model setup. 7 of 23 and density changes (e.g., Browning et al., 2015;Gottsmann et al., 2003;Gudmundsson, 2006). In the absence of precise data, we allocate a source pressure change (ΔP) of 10 MPa to the boundaries of the mush zone with ΔP matching the tensile strength of the crust (T 0 = 10 MPa) as proposed by Gudmundsson (2012). Furthermore, we assign a density change of 10 kg/m 3 to the mush zone resulting from the intrusion of relatively high-density magma into Ruapehu's mush zone as proposed by G. Kilgour et al. (2010), Nakagawa et al. (1999). To account for instantaneous source pressurization and density changes, we stepped ΔP and Δρ m at t = 10 −6 d, while ΔP and Δρ m are kept constant thereafter. Solid mechanics and gravity changes are modeled in the crust and edifice, while an additional domain above the free surface is required to simulate gravity changes ( Figure 1a). Boundary conditions for the solid mechanics solver are a free surface along the edifices topography, roller conditions (free of vertical displacement) at the right boundary, a fixed bottom boundary and a symmetry axis on the left side. Dirichlet boundary conditions for gravity changes are set to zero at the outer boundaries.
We solve poroelastic responses and strain-induced SP anomalies in the sub-domains (HTS, TZ, edifice) above 1.5 km a.m.s.l. using the approaches proposed in (Arens et al., 2020;Strehlow et al., 2015). The initial pore pressure distributions for strain-induced fluid flow at t = 0 is taken from the background HTU simulation (see Section 3.3.2). Boundary conditions are no-flow and an electric insulation with an electrical reference potential of φ = 0 V applied at r = 5.5 km and z = 1.52 km as the electrical potential is relative to a reference point. All internal boundaries for solid mechanics, gravity changes, poroelasticity and electrokinetics are continuous.

Hydrothermal Unrest Model
We use the TOUGHREACT code (Pruess et al., 1999;Xu et al., 2004, EOS2 module) to simulate HTU at Ruapehu by solving for heat and mass transport in porous media, but neglect reactive transport. The HTU model geometry ( Figure 1c) is confined to the edifice (z ≥ 1.5 km a.m.s.l and r < 5.5 km) with model domains being equivalent to the MU model setup. The HTS dimension has been envisaged by seismicity (Hurst, 1998), geochemistry (Christenson et al., 2010) and hyperspectral imaging (C. A. Miller et al., 2020).
We first perform a background simulation to establish the baseline condition prior to unrest at Ruapehu by injecting hydrothermal fluids (a mix of H 2 O and CO 2 ) at the bottom of the HTS (0 < r < 50 m) over a time span of 3,000 years. We then use the resulting distribution of pressure, temperature and gas saturation as the initial condition at t = 0 for the subsequent simulations of five unrest scenarios from anomalous injections each lasting for a period of 1 year (see Table 1 for fluid fluxes). We prescribe atmospheric boundary conditions along the ground surface (P = 0.101325 MPa, T = 10°C, p CO2 = 39 Pa), impermeable and adiabatic boundary conditions at the sides, and a basal heat flux of 0.086 W/m 2 (Stern et al., 1987) and impermeability at the base of the model, except for the points of fluid injection.
During multi-phase flow, gas pressures might differ from liquid pressures due to interfacial curvature and capillary forces, where the pressure difference between gas and liquid phases is equal to the capillary pressure . For HTU simulations, we define relative permeability and capillary pressure as function of liquid saturation (S l ) following Todesco et al. (2010). We calculate relative permeability κ rβ for phase β = l,g using the Corey function (Brooks & Corey, 1964;Pruess et al., 1999 and references within): with the effective saturation S e = (S l −S lr )/(1 − S lr −S gr ), the residual liquid saturation S lr = 0.33 and the residual gas saturation S gr of 0.05 Todesco et al., 2010). Capillary pressure is calculated as a linear function of S l for S > S lr , while capillary pressure is set to 0.01 MPa for S l < S lr (Todesco et al., 2010).
For HTU simulations, we test three different and distinct unrest scenarios (unrest I-III) using a set of injection rates given the absence of accurate observations. Scenario I represents the lowest injection rate while scenario III representing the highest (Table 1). Total injection rates are calculated from heat output for background and unrest activity given by Giggenbach and Glover (1975) using a fluid enthalpy of 3,000 kJ/kg (Hurst et al., 1991 Note. Total fluxes are calculated from heat outputs of Ruapehu's crater lake using a fluid enthalpy of 3 ⋅ 10 6 J/kg (Hurst et al., 1991). Heat outputs range between 200 MW for quiescence and 1000 MW during unrest (Giggenbach & Glover, 1975). Injection rates for H 2 O and CO 2 are derived from total fluxes using the fluids' molar ratios.   background and unrest simulations, respectively, while higher molar ratios are common during periods of unrest (e.g., Rinaldi et al., 2011;Todesco et al., 2010). Molar ratios for the hydrothermal fluids are chosen in accordance with CO 2 field observations at Ruapehu (see Text S5 in Supporting Information S1). While the CO 2 flux from unrest I is below the lower bound of recent emission records for unrest II it represents the upper bound. Note that the CO 2 flux for unrest III exceeds recent records (see Figure S3 in Supporting Information S1). However, by including such a flux allows us to study the detectability of geophysical anomalies from slightly higher CO 2 injection rates than recently observed, and identify whether certain anomalies become exclusively detectable at highest injection rates (unrest III).
Resultant deformation, gravity changes and SP anomalies are solved using the finite-difference method presented by Coco and Russo (2013) using the coordinate transformation method (Coco et al., 2014). The hydrothermal model dimension is extended radially and vertically for the electro-mechanical HTU simulations as shown in Figure 1c. Boundary conditions for displacement, gravity changes and SP simulations are set to zero at infinity, with an additional free-stress boundary conditions along the ground surface for deformation processes.

Parameterization
We parameterize our models with best-estimate or known values of subsurface conditions at Ruapehu as reported in Table 2. For model domains z > 0 km (edifice, HTS and TZ), we choose rock properties according to Ruapehu's andesitic deposits (Heap, Kushnir, et al., 2020;Mordensky et al., 2018). The HTS is represented by an altered, porous, permeable and water-saturated andesite, whereas the edifice is a stiff, dense, less permeable and less porous andesite. Hydraulic and electric rock properties for the TZ fall between values of the HTS and the edifice. Mechanical and thermal properties (e.g., E, λ) for the HTS, TZ and edifice are assigned in accordance with their porosity and water-saturation.
Thermal properties for the crust are chosen to match a greywacke composition (Mielke et al., 2016), while mechanical parameters such as the rock density (ρ c ) and the dynamic Young's modulus (E d ) are derived from 2D seismic P-wave velocities (Rowlands et al., 2005) using the Brocher (2005) relationships (Equation S3-S4 in Supporting Information S1) and a Poisson's ratio of 0.25 ( Figure S2 in Supporting Information S1). We convert E d to static modulus (E s ) using a conversion of E s = 0.5*E d (e.g., Cheng & Johnston, 1981;Gudmundsson, 1983). We fit crustal density and static Young's modulus (E c ) by a third-order polynomial to obtain a continuous function of depth (z).
We set the Biot-Willis coefficient equivalent to the domains rock porosity. In the absence of precise data, we vary α BW in the parameter study (see below) between the rock porosity and 1 for soft materials according to  Calculated (see Text S2 in Supporting Information S1) Note. HTS-hydrothermal system, TZ-transition zone, e-edifice, c-crust, f-fluid, MZ-mush zone, m-magma, SP-Self-potential.
10 of 23 Wang (2000). The electric conductivities of the sub-domains HTS, TZ and edifice are taken from magnetotelluric studies by Ingham et al. (2009) and Jones et al. (2008). In the absence of direct measurements of C SP for Ruapehu, we derive C SP after Revil and Pezard (1998) (Equation S4 in Supporting Information S1) with the fluid conductivity being calculated according to Byrdina et al. (2018) using Ruapehu's crater lake chemistry (see Text S4 in Supporting Information S1).
We test the influence of selected parameters on modeled unrest anomalies by exploring plausible value ranges of parameters for which either only sparse or no data exist for Ruapehu. For all unrest simulations we investigate the effect of α BW (α BW × 2, ×4) and C SP (±10 −8 −±10 −10 V/Pa) on geophysical anomalies individually by varying these parameters in all sub-domains above z = 1.5 km. Additionally, we test different reservoir volumes (8.9-35.7 km 3 , Table S1 in Supporting Information S1) with the reservoir strengths (VxΔP) being equivalent across all volumes tested and source density changes (Δρ m = 10-300 kg/m 3 ) for MU simulations.

Results
Here we present the results of the unrest simulations. We report the solutions for the temporal evolution of unrest observables on Ruapehu's summit plateau with coordinates r = 500 m and z = 2,640 m. We choose the summit plateau due to its flat topography and the opportunity to capture near-field effects from unrest whilst also accounting for operational safety (V. Miller et al., 2003). Figure 2 shows geophysical anomalies for the reference parameterization of simulated MU. Along the ground surface, we find peak vertical displacements at a distance of 3.5 km from the symmetry axis, with a maximum uplift of 0.87 cm attained 100 days after source pressurization. A localized maximum occurs above the HTS (at r < 50 m) at the bottom of the crater lake with values 0.9 times the maximum amplitude (Figure 2a). Horizontal displacements peak at a distance of ∼6.75 km from the deformation center with amplitudes of ∼0.70 cm (see Figure S9a in Supporting Information S1). Residual gravity changes are at their maximum above the HTS (8.5 μGal at t = 1 d, Figure 2e) with a linear decrease in signal magnitude with distance from the HTS. Concurrent SP anomalies peak above the HTS (SP max ∼ 1.3 mV at t = 10 days) with anomalies rapidly falling off to negative values at distances r < 400m from the HTS. The temporal evolution of the unrest observables at the plateau (r = 500 m and z = 2,640 m) are illustrated in Figure 2 (lower panels). We find a maximum amplitude change of 0.10 cm for vertical displacement (w) in the edifice within the first 30 days after source pressurization, with peak magnitudes of 0.76 cm at t = 30 days. Horizontal displacement (u) shows an initial minimum of 0.04 cm followed by a continuous increase in magnitude with time. After 350 days (referred to as 1 year hereafter) maximal u at the plateau is 0.07 times smaller than the peak magnitude (u ∼ 0.7 cm) at r = 6.75 km (Figure 2d and Figure S9b in Supporting Information S1). Residual gravity changes (δg r ) decrease rapidly within the first 30 days with a maximum change of 0.5 μGal. SP anomalies decrease linearly with time showing an absolute amplitude change of 0.7 mV.

Parameter Exploration
For large reservoir volumes (Figure 3 left panels) >7.2 km 3 /MPa, we observe magnitudes of 1.3, 1.4, and 2.1 times the initial reference values for vertical displacement, horizontal displacement and residual gravity changes, respectively. SP anomalies remain broadly unchanged for the explored reservoir strength variations.
An increase in Biot-Willis coefficient (α BW , Figure 3 right panels) reduces vertical displacements relative to reference values throughout time, while horizontal displacements are reduced for t < 75 days but increased thereafter. Residual gravity changes and SP anomalies are amplified with respect to the initial reference values for the highest α BW tested (α BW × 4). Figure 4a shows the impact of varying the source density change (Δρ m ) on residual gravity changes (δg r ). We find a correlation between peak δg r and Δρ m , with greatest δg r of 35 times reference amplitudes for Δρ m = 300 kg/m 3 . Polarity of δg r corresponds to Δρ m polarities.
The influence of the streaming-potential coupling coefficient (C SP ) on SP anomalies is shown in Figure 4b.
For the highest C SP tested (10 −8 V/Pa), SP magnitudes are amplified up to 10 times the initial reference value, while SP amplitudes for ±C SP = 10 −10 V/Pa show negligible changes. SP time series for negative and positive C SP polarities are inverted. 12 of 23

Hydrothermal Injection
Panels a, e, and i in Figure 5 depict the initial pore pressure, temperature and gas saturation distributions after background injection over 3,000 years. We find peak pore pressures around the injection area (r < 50 m, z = 1.5 km) whereas the far-field pore pressure distribution mirrors the topography. Temperature and gas saturation are elevated around the HTS, with highest temperatures of ∼300°C at the bottom of the HTS. Maximal gas saturation are simulated below the crater lake (z ∼ 2,440 m, r < 200 m).
Relative to background distributions, variations in pore pressure (Δp f ), temperature (ΔT) and gas saturation (ΔS) correlate with anomalous fluid fluxes ( Figure 5). We observe maximum (Δp f ) of 9 MPa, ΔT of 60°C and ΔS of 0.6 for unrest III after 1 year of anomalous injection ( Figure 5 panels d, h and l). For unrest I, concurrent amplitudes changes of (Δp f ), ΔT, and ΔS are between 0.1 and 0.18 times maximum amplitudes of unrest III. Variations in pore pressure, temperature and gas saturation are confined to the HTS and its proximity, whereas the far-field is broadly undisturbed. For increasing injection rates, the gas plume propagates further toward the surface and additionally dispersing into the TZ at z = 1.5 km.

Hydrothermal Unrest Anomalies
Figure 6 (panels a-h) shows the results of simulated HTU anomalies along the ground surface at different times since anomalous injections. We find peak (positive/negative) anomalies except for horizontal displacements (maximum u at r > 50 m) directly above the HTS with magnitudes falling off rapidly with distance to the HTS.
After 1 year of anomalous injection, maximum vertical uplift of ∼3 cm and horizontal displacement (u) of ∼1.3 cm is observed for unrest III, while unrest I induces magnitudes 14% and 9% of maximum w and u, respectively. Residual gravity changes are of negative polarity with minimal values above the HTS ranging between −83 μGal (unrest III) and −18 μGal (unrest I) at t = 350 days and t = 100 days, respectively. SP anomalies ( Figure 6, panels d and h) peak above the HTS with the maximum of 6.5 mV (1 year) corresponding to the highest injection rate, while SP anomalies resulting from unrest I are only of 0.63 mV.
Time series of simulated HTU anomalies at the plateau are illustrated in Figure 6 (panels i-l). Vertical and horizontal displacements exhibit similar temporal evaluations with increasing magnitude with time and a maximum ARENS ET AL.

10.1029/2022GC010572
13 of 23 of ∼0.86 cm for w and u for unrest III at t = 350 days. Residual gravity changes decrease monotonically with time for unrest I, whereas δg r for unrest ≥II show time-delayed minima at 100 (δg r = −7.8 μGal greatest minimum) and 200 days for unrest III and II, respectively. SP time series fluctuates throughout time with an overall increase of SP amplitudes for unrest ≥II with time, while SP anomalies for unrest I remain broadly unchanged. Maximum SP amplitudes of 2.6 mV are observed for unrest III, which is 8.5 times the SP magnitude of unrest I.

Parameter Exploration
The upper panels in Figure 7 show the influence of Biot-Willis coefficient (α BW ) on temporal ground displacements at the plateau. Displacements correlate positively with α BW . For the largest values tested (α BW × 4), we obtain vertical displacements up to 3.5 times (unrest III) higher than the maximum reference amplitude, with similar changes in magnitude for horizontal displacements. Residual gravity changes decrease with increasing α BW with minimum values of −13 μGal for unrest III, which is 1.6 times the reference minimum. Figure 7 (lower panels) shows the effect of the streaming-potential coupling coefficient (C SP ) on temporal SP anomalies. While C SP values strongly control the SP amplitudes and polarities, the temporal evolution mirrors the reference time series (Figure 6l) for positive C SP , but is inverted for negative C SP polarities. We find that SP anomalies vary up to a factor of 10 (unrest III) smaller or greater than that of the reference SP amplitude for negative or positive C SP , respectively. , temperature (f-h) and gas saturation (j-l) with respect to the background simulation are illustrated for 1 year of anomalous injection. Unrest III represents the highest injection rate (use Table 1

Discussion
Our multiphysical modeling approach is the first study investigating multi-parametric anomalies from MU and HTU processes at Ruapehu. We have shown that magmatic and hydrothermal perturbations induce markedly different spatio-temporal observables. Simulation results depend strongly on underpinning model assumptions and parameterization which in our study are constrained by geophysical, geological and petrological data.

Magmatic Unrest Anomalies
While spatial ground displacement patterns from MU simulations are broadly similar to predictions from time-independent elastic half-space solutions for a prolate magma reservoir at Ruapehu (V. Miller et al., 2003), we note several key differences: (a) a non-linear evolution of ground displacements due to poroelasticity in the edifice and crustal viscoelasticity (see Figure 2 and Figure S4 in Supporting Information S1), (b) poroelasticity in the HTS and TZ reduces the magnitude of vertical displacement w for r < 200 m and (c) reduced (by up to 50%) magnitudes of ground displacements in our study. The latter corroborates results reported in Males and Gottsmann (2021) where subsurface heterogeneity and volcano prominence control the stress and strain partitioning and hence the displacement magnitudes. Additionally, the displacement magnitude is controlled, as expected, by elastic parameters, source pressure (ΔP) and the location and dimension of the magmatic reservoir (e.g., Hickey et al., 2013).
Subsurface displacements influence residual gravity changes through the gravity contributions from host rock compression and shifting density boundaries (Equation 11). However, these contributions are of minor importance in our study as δg r is predominantly governed by source density changes (see Figure S5 in Supporting Information S1), corroborating findings reported in Gottsmann, Biggs, et al. (2020). As Δρ m remain constant  Figure S10 in Supporting Information S1). throughout time, the temporal evolution of δg r is opposite to the temporal evolution of w due to the free-air effect.
In terms of spatial patterns, we find an agreement of δg r along the ground surface between our study and findings in Currenti (2014), with peak amplitudes directly above the HTS. The temporal evolution of δg r is similar to that of the vertical displacement governed by poroelastic responses of the edifice and viscoelastic processes in the crust. Visco-poroelastic processes appear to dominate ground deformation at the beginning of the perturbation (see Figure 2) with ground subsidence following initial uplift. This compares to subsidence only in simulations accounting for poroelastic effects (see Figure S4 in Supporting Information S1). However, given the resolution limit of geodetic observations neither process is detectable.
Self-potential anomalies from strain-induced fluid flow peak above the HTS, where pore pressure variations are at their largest (Equation 13). We find an absolute SP change of 0.7 mV after 1 year of magmatic perturbation. The continuous decrease in SP magnitude with time indicates a drop in pore pressure as C SP and pore pressure are positively correlated in this study, as opposed to the inverse relationship for non-acidic waters described elsewhere (Revil, Saracco, & Labazuy, 2003;Rizzo et al., 2004).

Parameter Exploration
The parametric study revealed minor variations in ground displacement magnitudes with changing reservoir volumes and Biot-Willis coefficients (α BW ). Since reservoir strength is kept constant in all simulations, resultant ground deformations are controlled by visco-poroelastic responses of the surrounding media to induced pressure perturbations compared to reference solutions. In a one-way coupling approach ground displacement correlates with pore pressure changes and α BW (Currenti & Williams, 2014;Raziperchikolaee et al., 2020). However, in our two-way coupling approach where Δp f affects stresses and strains and vice versa, the effect of α BW on radial displacements in particular is more complex. Ground displacements are generally controlled by stress and pore pressure changes in response to subsurface heterogeneities (see Hickey and Gottsmann, 2014;Strehlow et al., 2015). Residual gravity changes are strongly influenced by changes in source density and pressure (Δρ m and V, respectively). We show that Δρ m and δg r correlate in terms of magnitude. In our MU models, δg r are primarily governed by the increase in source density as a result of the injection of new magma, a common assumption behind episodes of unrest at Ruapehu (G. N. Kilgour et al., 2013;Nakagawa et al., 1999). Source density changes of 10 kg/m 3 in greater reservoir volumes results in larger δg r values compared to the reference simulation. The dependency of δg r on α BW is negligible across the range of tested values, as the effect of α BW on displacements of density boundaries is much smaller than the effect of source density changes (see Section 5.1.1).
Similar to findings in Arens et al. (2020), our study finds that SP anomalies are governed by electrokinetic processes arising from poroelastic responses to subsurface perturbations and are hence primarily controlled by α BW . Furthermore, we show that SP magnitudes are markedly controlled by C SP matching results reported by Arens et al. (2020), where C SP is categorized as an influential parameter.

Hydrothermal Injection
The injection of hydrothermal fluids disturbs the physicochemical conditions in the subsurface and manifests as variations in pore pressure, temperature and gas saturation in the subsurface ( Figure 5). Comparing our results with findings reported in Christenson et al. (2010) we note differences in model parameterization (e.g., injection rates, parameters), model setup (e.g., flat surface, initial conditions) and HTS volume compared to our study. Although simulated background gas saturation and temperature distributions in our study broadly resemble the background conditions of Ruapehu considered in Christenson et al. (2010), temporal changes in the gas and temperature distribution in our study are predicted over a much larger space. This might result from a wider injection area and a longer-lasting injection period compared to the study of Christenson et al. (2010). Unlike the linear pore pressure evolution in other studies (e.g., Christenson et al., 2010;Stissi et al., 2021), we simulate elevated initial pore pressures around the HTS and its proximity after protracted background injection. As a result pore pressures reach ∼10 MPa around the injection area and are similar to the pore pressure parameterization in . The overall pattern of pore pressure distribution mirrors topography, indicating that topographic effects must be taken into account when investigating fluid flow in a volcanic edifice (Figure 5a).
Transient variations of pore pressure, temperature and gas saturation caused by anomalous injection are confined to the injection area (HTS) and its proximity as observed by . Similar to findings reported in Christenson et al. (2010), pressure and temperature pulses (relative to the background) propagate toward the crater lake bottom over time ( Figure 5, and Figures S7 and S8 in Supporting Information S1). Note though, that in contrast to the simulated intrusion of gases into the crater lake in Christenson et al. (2010), in our study a deep-seated single-phase gas plume develops in the HTS (see Figure 5 and Figure S6 in Supporting Information S1). Furthermore, the drop in ΔS below the crater ( Figure 5 panels j-l, r < 200 m) indicates that liquid H 2 O enters previously gas-enriched areas as it migrates quickly through permeable domains (e.g., Todesco et al., 2010), such as the HTS and TZ. The positive correlations between injection rates and magnitudes of Δp f , ΔT, and ΔS match results reported in .

Hydrothermal Unrest Anomalies
We find similarities in the spatial displacement patterns from HTU simulations to findings reported in Stissi et al. (2021) and . Here, highest vertical displacements correlate with largest pore pressure and temperature variations ( Figure 5). As peak (Δp f ) and ΔT values are encountered around the HTS, vertical displacements fall off rapidly with distance from the HTS. We find that ground displacements increase with both injection rate and time due to the thermo-poroelastic response caused by protracted pore pressure and temperature variations. That is to say that ground displacements evolve in unison with the severity of HTU matching findings in .

of 23
In contrast to ground displacements, magnitudes of residual gravity changes correlate negatively with fluid fluxes. The spatio-temporal behavior of δg r is controlled by fluid density variations (e.g., Todesco, 2009;Todesco & Berrino, 2005). That is to say, in areas where H 2 O replaces gases (e.g., in the TZ; increase in Δρ f ) positive δg r are expected, while negative δg r arise where gas-rich fluids ascend (e.g., in the HTS; drop in Δρ f ). Subsurface heterogeneities strongly govern the distribution of H 2 O and CO 2 (Todesco et al., 2010). For instance, the permeable HTS favors the upwards migration of H 2 O and CO 2 due to influx of new fluids at its base, which might prompt the discharge of H 2 O at the surface causing an overall decrease in Δρ f and δg r . This behavior could explain the negative δg r values directly above the HTS.
Spatio-temporal H 2 O and CO 2 fluctuations govern (Δp f ) and hence electrokinetic processes. SP magnitudes correlate with (Δp f ) for the strongest HTU, with an overall increase in SP amplitude over time for protracted HTU (>200 days). We find that the spatial SP pattern matches observations at other volcanoes (volcano-electric effect after Revil, Saracco, and Labazuy (2003)) with peak SP anomalies directly above zones of hydrothermal upflow (Zlotnicki & Nishida, 2003).

Parameter Exploration
We show that the Biot-Willis coefficient influences geodetic anomalies from hydrothermal perturbations. Ground displacements are governed by the poroelastic response of the one-way coupling approach (Equation 7) and are hence controlled by Δp f and α BW . That is to say, that uplift correlates with Δp f and α BW (see also Raziperchikolaee et al. (2020)). The choice of poroelastic coupling (one-way vs. two-way) could explain the different effect of α BW on displacements between HTU and MU simulations. The influence of α BW on δg r is predominantly caused by the gravity contributions from the free-air effect and hence w; that is, δg r magnitudes decrease for increasing α BW (and w). SP anomalies are not governed by α BW in HTU simulations due to the one-way coupling approach. Like in our MU parametric simulations, SP magnitudes from hydrothermal perturbation correlate with the key parameter C SP .

Implications for Geophysical Unrest Monitoring at Ruapehu
While changes in ground elevation are routinely monitored at Ruapehu, monitoring of SP and gravity changes is absent. It is interesting to note that prior to the most recent magmatic eruption at Ruapehu in 2007, no ground displacements were observed (Mordret et al., 2010). To explain this and to identify geophysical anomalies indicative of MU or HTU, we compare the simulated magnitudes of ground displacements as well as SP and gravity changes with detection levels of conventional surveying techniques. Our analysis is focused on the near-field of the crater lake at the plateau (r = 500 m, z = 2,640 m) where instrumentation could be deployed and maintained.
Ground displacements from MU simulations remain below the detectability limits of 1 cm vertically and 0.5 cm horizontally by GNSS surveys (Mordret et al., 2010) on the plateau. However, horizontal displacements become detectable at a distance of 6.75 km from the HTS (see Figure S9 in Supporting Information S1). Our parametric investigations from MU simulations show that ground displacements at the plateau remain below detection levels even at the largest magmatic perturbation explored in this study. The detectability of ground displacements from HTU simulations is complex and differs for horizontal and vertical displacement. Although w displacements in reference simulations are below conventional detection limits, u displacements are detectable for unrest ≥II. For unrest III conditions, u exceed detection limits after a much shorter period of time (at t ∼ 100 days) compared to unrest II (at t = ∼300 days). As such the geodetic detectability of unrest depends on the magnitude of subsurface perturbations. For the largest Biot-Willis coefficient explored in this study, the peak u displacement becomes detectable after a shorter time compared to the reference simulations (e.g., unrest III at t ∼ 50 days vs. unrest II at t = ∼90 days). Additionally, w during unrest ≥II exceeds detection limits for all α BW values tested. The absence of pre-eruptive displacement anomalies in the most recent phreatic eruption might be explained by anomalous hydrothermal injection at rates similar to conditions simulated in unrest I.
Residual gravity changes from MU reference and most parametric simulations are above detection levels of ±5 μGal (Battaglia et al., 2008), but their temporal variations remain undetectable. For HTU reference simulations, injection rates ≥unrest II induce measurable δg r after >40 days of anomalous injection, while δg r from fluid fluxes I remain undetectable throughout. Higher α BW values result in higher δg r values and favor their detectability.

of 23
The SP anomaly from protracted unrest reaches an absolute change of 0.7 mV (MU simulations), while SP anomalies range between <0.5 to a maximum of ∼2.5 mV for unrest I and III in the HTU simulations, respectively. SP magnitudes from subsurface perturbations fall within the detectability levels of standard field observations (0.1 mV; Grobbe and Barde-Cabusson, 2019;Revil and Jardani, 2013). Parametric studies for both unrest scenarios have shown that SP magnitudes increase significantly with increasing the streaming-potential coupling coefficient.
Although some simulations predict ground displacements, gravity changes and perturbations in SP above detectability limits, the temporal evolutions of the signals are predicted to be difficult to resolve. However, some combinations of observables are indicative of source processes. For example, a temporal decrease in w displacements and simultaneous increase in u displacements might indicate magma pressurization and the time-dependent visco-poroelastic response of the surrounding media. The temporal evolution of δg r is similar in both HTU and MU simulations whereby the signal amplitude decreases initially followed by an increase. However, δg r values are positive in MU simulations and negative in HTU simulations. Fluid density distribution from HTU simulations depends on the spatio-temporal distribution of gas and liquid in the subsurface and fluctuates as a result of fluid injections and redistribution. Therefore the change in magnitude of δg r with time is more pronounced in HTU simulations compared to MU simulations (see Figures 2 and 6). Self-potential anomalies decrease in MU simulations with time but increase in HTU simulations (see Figures 2 and 6).
We identify distinct sets of detectable geophysical anomalies at Ruapehu's plateau which could be used to interrogate the nature of volcanic unrest. We find that density changes in the crustal mush zone and electrokinetic processes from strain-induced fluid flow in the volcanic edifice (z > 1.5 km) induce measurable gravitational and electrical potential field anomalies at the plateau and hence are indicative of MU. Horizontal displacements in the far-field might act as additional indicators of source pressurization (MU simulations), but ground displacements in the proximity of the HTS are not detectable. Protracted HTU is identifiable by SP anomalies for all HTU simulations and ground displacements for unrest ≥II. In addition, residual gravity changes become a distinctive fingerprint of HTU for CO 2 fluxes matching those during the 2007 unrest (i.e., unrest II). This implies that protracted HTU at the higher end of CO 2 fluxes explored in our models (unrest III) yields detectable residual gravity changes. We therefore recommend the implementation of continuous gravity and SP monitoring at Ruapehu, which in combination with existing monitoring techniques (e.g., seismicity, fluid chemistry) at Ruapehu could significantly improve interpretations of source processes during unrest periods. The summit plateau would be suitable to safely locate monitoring instrumentation (V. Miller et al., 2003); based on our findings a combination of the three geophysical signals from either magmatic or hydrothermal perturbation is detectable. As continuous GNSS sites at Ruapehu are located >500 m from the HTS (http://www.geonet.org. nz), we suggest the implementation of GNSS sites at the summit plateau to allow for signal detectability (e.g., HTU). Most signals fall off rapidly with distance from the HTS; locating monitoring sites at r > 500 m drastically reduces signal detectability. At the same time, installing and maintaining monitoring stations closer to the HTS could be challenging due to the steep topography and potential impact of ballistics during eruptions (G. Kilgour et al., 2010;Strehlow et al., 2017).

Model Limitations
We use a simplified model geometry (2D axisymmetrical) to keep simulations computationally cost-efficient, but sufficiently complex to gain first-order insights into geophysical anomalies caused by MU and HTU at Mt. Ruapehu. Both unrest processes are studied in isolation, while in reality magmatic and hydrothermal perturbations might superimpose. Furthermore, we do not account for the interaction of magma with the HTS. All models presented in this study incorporate subsurface mechanical, electrical and hydraulic heterogeneity and account for a topography representative of the volcano. All of the multi-parametric data sets that helped constrain our models are either 1D or 2D. Should 3D variations of these parameters become available, the models can be adapted to provide 3D solutions.
Inherent model limitations for MU simulations have been described in detail in Arens et al. (2020). Our HTU simulations do not account for super-critical conditions, although there is evidence for pressures and temperatures in HTSs at active volcanoes exceeding the critical point of water (Reinsch et al., 2017). The thermo-poroelastic coupling approach used in this study is most representative of short-term hydrothermal perturbations  with applications to many volcanoes (e.g., Currenti & Napoli, 2017; Fournier &

of 23
Chardot, 2012; Todesco & Berrino, 2005). However, it has been shown that a two-way coupling approach is more applicable for temporally protracted perturbations (Neuzil, 2003;Rutqvist et al., 2002), where subsurface strain affects hydraulic rock properties (e.g., κ, ϕ) which in turn govern the flow behavior and in turn stresses and strains. Similar to studies by for example, Hutnak et al. (2009);; Fournier and Chardot (2012); Rinaldi et al. (2011), we neglect the effect of (a) shifting density boundaries and (b) host rock compression on residual gravity changes. Gravity contributions from (a) and (b) in our study are 0.2 μGal and −0.4 μGal, respectively, and hence almost cancel one another out. Therefore we deduce fluid density changes from hydrothermal perturbations as the main source of δg r changes. The inclusion of the aforementioned effects and the two-way coupling approach would be a next step of studying HTU at Ruapehu.
Neither of our simulations account for the temperature dependence of parameters such as permeability (Ikard & Revil, 2014), fluid properties (Arens et al., 2020) or elastic parameters (Head et al., 2021), all of which have an effect on geophysical anomalies modeled in our study; a dedicated analysis is required to assess this influence. Although we neglect thermoelectric processes caused by strong thermal gradients (Corwin & Hoover, 1979;Fitterman & Corwin, 1982) in the HTU simulations, we find that for a maximum temperature change of 0.18°C (unrest III) at the plateau, the thermoelectric potential (TEP) is ±0.3 and 0.1 mV using a thermoelectric coupling coefficient of ±1.5 mV/°C and ±0.5 mV/°C (Ikard & Revil, 2014;Revil & Mahardika, 2013), respectively. The TEP is only 5-15% of the maximum SP amplitude, so we conclude that electrokinetic processes dominate electrical potential field changes.

Conclusions
We have utilized multiphysics models to study volcanic unrest and concurrent geophysical anomalies at the active volcano Mt. Ruapehu. Our study was able to discriminate spatio-temporal anomalies that might help identify the nature of unrest (hydrothermal vs. magmatic). While gravitational and electrical potential field anomalies are indicative of magmatic processes (e.g., source pressurization and density changes) in the sub-volcanic mush zone, ground displacements (vertical and horizontal) in the proximity of the deformation source remain below detection limits for reference and parametric simulations. However, horizontal displacements become resolvable in the far-field and could provide additional insights into MU. In contrast, ground displacements, residual gravity changes and SP anomalies from HTU are detectable in the near-field.
Parameter space testing show the major control of some key model parameters (e.g., α BW , C SP , V) on the detectability of geophysical anomalies. For instance, magnitudes of SP and residual gravity changes correlate with key parameters C SP and Δ m , respectively. While the superposition of magmatic and hydrothermal perturbations need to be taken into account when interpreting observed precursors, we have identified unique sets of resolvable magnitudes of geophysical anomalies from either subsurface perturbation. We conclude that joint and simultaneously collected multi-parameter time series should provide valuable insights into unrest source mechanisms, especially when corrected for non-volcanic background processes. In order to distinguish between frequent HTU and less-frequent but potentially more violent MU at Ruapehu, we propose the implementation of routine SP and gravity monitoring to support ongoing monitoring efforts. The findings reported here may have implications for assessing unrest dynamics at other crater lake volcanoes.