Poroelastic responses of confined aquifers to subsurface strain and their use for volcano monitoring.

. Well water level changes associated with magmatic unrest can be interpreted as a result of pore pressure changes in the aquifer due to crustal deformation, and so could provide constraints on the subsurface processes causing this strain. We use ﬁnite element analysis to demonstrate the response of aquifers to volumetric strain induced by pressurized magma reservoirs. Two different aquifers are invoked – an unconsolidated pyroclastic deposit and a vesicular lava ﬂow – and embedded in an impermeable crust, overlying a magma chamber. The time-dependent, fully coupled models simulate crustal deformation accompanying chamber pressurization and the resulting hydraulic head changes as well as ﬂow through the porous aquifer, i.e. porous ﬂow. The simulated strain leads to centimetres (pyroclastic aquifer) to metres (lava ﬂow aquifer) of hydraulic head changes; both strain and hydraulic head change with time due to substantial porous ﬂow in the hydrological system. Well level changes are particularly sensitive to chamber volume, shape and pressurization strength, followed by aquifer permeability and the phase of the pore ﬂuid. The depths of chamber and aquifer, as well as the aquifer’s Young’s modulus also have signiﬁcant inﬂuence


Introduction
Pre-, syn-and post-eruptive changes in water levels have been reported for several volcanoes (Newhall et al., 2001).Many processes can lead to water level changes, including meteorological influences (e.g.rainfall or barometric pressure), the injection of magmatic fluids (e.g.suggested for Campi Flegrei by Chiodini et al., 2012), the opening and closing of fractures (e.g. proposed for Kilauea by Hurwitz and Johnston, 2003) or heating of the aquifer (suggested for Campi Flegrei by Gaeta et al., 1998).
Here, we focus on the very commonly suggested mechanism of strain-induced water level changes.Examples include well level changes of more than 9 m preceding the 2000 eruption of Usu volcano, Japan (Matsumoto et al., 2002) and the water level rise of more than 85 m in a geothermal well at Krafla volcano, Iceland, associated with a dyke intrusion in 1977 (Stefansson, 1981).The observations can be explained by poroelasticity (Wang, 2000).Compression or dilatation of an elastic porous medium leads to a decrease or an increase in pore space, respectively, which in turn influences the pore pressure and thereby the water level.Hence, measured well level changes can be interpreted as strain-induced changes in pore pressure in the aquifer due to crustal deformation.This is in line with observations of water level changes accompanying seismic events (e.g.Roeloffs, 1996;Jonsson et al., 2003;Shibata et al., 2010) or crustal spreading, as observed at the Juan de Fuca Ridge (Davis et al., 2001).
In volcanic environments, many processes can lead to substantial strain, including pressure changes in magma reservoirs and intruding dykes.Information about the local strain field is therefore highly valuable for volcano monitoring and eruption forecasting, as it could allow derivation of these subsurface magmatic processes (e.g.Voight et al., 2006;Linde et al., 2010;Bonaccorso et al., 2012).However, strain data are difficult to interpret and strainmeters are complex and expensive installations.The described poroelastic relations raise the question whether wells in aquifers can provide additional information on the subsurface strain field and if we could even use them as cheaper and somewhat simpler strainmeters.
Previous studies have indeed utilized the poroelastic behaviour of aquifers to infer magmatic processes from observed water level changes at volcanoes (e.g.Shibata and Akita, 2001;Takahashi et al., 2012).A method of assessing the strain sensitivity of an aquifer is to track water level changes as a result of predictable excitations such as Earth tides or measured barometric variations.The known strain sensitivity is then used to derive volumetric strain from observed water level changes during unrest, and combining this with analytical deformation models such as the Mogi model (Mogi, 1958), inferences can be made on magmatic drivers behind the level changes.
However, oversimplification of the coupling between solid and fluid mechanics may make these models inadequate.An example is the 2000 Usu eruption (Matsumoto et al., 2002), where pre-eruptive water table changes of several metres were observed in two wells at different locations simultaneously with a radial ground deformation of about 2 cm recorded about 8 km from the summit over the course of 2 weeks prior to the eruption.The water level changes were interpreted as a result of this crustal deformation.However, the two different wells apparently give inconsistent information about the source of strain: only one of the two well level changes agrees with the model proposed by Matsumoto et al. (2002).In order to make reasonable monitoring interpretations based on well level data, we therefore need to improve our understanding of how these hydrological signals are generated and identify the relative importance of the parameters that affect them.Changes in the hydrological conditions in volcanic areas are usually interpreted as a result of changes in the magmatic system, but the effects of non-magmatic parameters on the pressure-response in the aquifer should also be considered.
Analytical solutions exist for only a few comparatively simple poroelastic problems (e.g.Rice and Cleary, 1976).Numerical modelling of pressure changes in hydrological systems has focused on pressure and temperature transients in hydrothermal systems and resulting ground deformation due to the injection of hot magmatic fluids, using one-way coupling of solid deformation and porous flow (e.g.Todesco et al., 2004;Hurwitz et al., 2007;Hutnak et al., 2009;Rinaldi et al., 2010;Fournier and Chardot, 2012).Rutqvist et al. (2002) have developed a two-way coupled code and applied it to problems related to carbon dioxide injection in aquifers and the disposal of nuclear waste in porous media.Whether and how the pure deformational strain caused by magma bodies would induce water level changes has not yet been explored numerically.The full, two-way coupling of fluid and solid mechanics required has so far been avoided in volcanological applications, and so the effect of solid deformation on pore pressure and porous flow has been neglected.
We investigate the phenomenon of poroelastic responses to magmatic strain to better understand the hydrological signals one might observe in wells on a volcano before and during eruptions.We assess to what extent confined aquifers can serve as indicators of stress/strain partitioning in the shallow crust due to reservoir pressure changes and therefore if they could provide a tool to scrutinise pre-eruption processes.

Methods
Table 1 gives a list of all symbols used in this study.

Theory
We present a set of generic models using finite element analysis to perform parametric studies on several volcanic settings with an inflating magma chamber affecting overlying rock layers and hydrology.The models solve a series of constitutive equations that result from the full coupling of continuum mechanics equations for stress-strain relations of a linear elastic material with Darcy's law and mass conservation within the porous flow theory (for details see Biot, 1962;Wang, 2000;COMSOL, 2013).The calculations are based on the Navier equation for a solid: with σ being the stress tensor and F V a body force.Inertia terms in the Navier equation are neglected as the solid deformation is treated as quasi-static.The solid mechanics equations assume linear elasticity and do not allow for material failure, hence only work for sufficiently small strains.
The stress tensor σ is related to the strain tensor and the pore pressure p f by a generalized Hooke's law: Here, C is the drained elasticity tensor and α is the Biot-Willis coefficient.Strain is given through the displacement vector (u): Fluid flow is described by mass conservation and Darcy's law: Here, ρ f is the fluid density, S is the specific storage, v is fluid flow velocity, Q is a source/sink term, vol is the volumetric strain, κ is aquifer permeability, µ is water viscosity, ρ f is water density, g is acceleration of gravity, and D is elevation.The equations for fluid flow only consider single-phase, single-component flow.The aquifer is considered to be fully saturated and perfectly confined at all times.
In both Eqs. ( 2) and ( 4), the terms including the Biot-Willis coefficient describe the coupling between solid deformation and fluid flow, which manifests in stress absorption by the fluid and pore pressure changes due to the increase/decrease of pore space resulting from volumetric changes of the porous medium.The coupling parameter α is a measure of the strength of this coupling (with values between the porosity of the medium and 1), and is defined by the volume of fluid expelled from/sucked into a porous medium when subject to volumetric change.Fluid mechanical properties remain unchanged, including permeability and porosity.The coupling is achieved solely by the pore pressure and stress effects, as well as the expression for the specific storage, which includes elastic properties of both pore with φ being the porosity of the medium, χ f the fluid compressibility, and K the drained bulk modulus of the solid matrix.This set of equations is solved for solid deformation (u) and fluid pressure (p f ) using the structural mechanics and Darcy's law modules of the finite element analysis software package COMSOL Multiphysics, version 5.0.These modules provide the equations for solid deformation and fluid flow, respectively, and have been coupled by manually including the coupling terms in Eqs. ( 2) and (4) as described above and appropriately defining the specific storage.It should be noted that the readily provided poroelasticity module (which applies the same equations) has been avoided, as this causes problems when gravity is not neglected.In the provided coupling, COMSOL treats solid and fluid as one unit, meaning they compute the average density of the porous domain and gravity then acts on this.This causes downward displacement of the solid matrix and the gravity part of Darcy's law causes downward porous flow.While it is possible to equilibrate this flow by initializing hydrostatic pressure (using fluid density), this does not prevent the solid displacement, as this needs to be balanced using initial lithostatic pressure (using average density).As one needs to use different densities it is not possible to equilibrate flow and displacement at the same time in this module when gravity is turned on.In our solution, we switch on gravity only in the Darcy's law module, hence the solid matrix is not affected by the gravitational stresses, and the fluid pressure is initialized as hydrostatic.This way we can take flow following topographic gradients into account.

Model set-up
As a starting point to investigate hydrological responses to magma chamber inflation, we build a 2-D-axisymmetric model geometry in COMSOL Multiphysics following Hickey and Gottsmann (2014), who provide guidelines for volcano deformation modelling using finite element analysis.The initial model consists of a linear elastic solid block with an embedded spherical cavity, representing a magma chamber at depth.This cavity is pressurized by applying a boundary load, which is stepped up over 10 −8 s.We assume almost instantaneous pressurization for simplicity and to more easily recognise the different influences of parameters.Magma chamber pressurization can be generated by the injection of fresh magma, vesiculation, thermal expansion of the magma, melting of country rocks or volume changes during crystallization (Fagents et al., 2013).Using the relation for temperature-independent volume changes and assuming a magma compressibility of β = 10 −11 Pa −1 , a pressurization of 10 MPa could correspond to a volume change of V = 100 000 m 3 .V is not simulated in the presented models, as the magma chamber is represented by a pressurized cavity.The volume change calculated with Eq. ( 7) serves as a first order estimation and a guide to corresponding magmatic processes.Note however, that this equation does not account for several additional processes and does not deliver a perfectly accurate volume change of the source.The resulting deformation of the surrounding material is calculated by discretizing the model domain to solve the constitutive equations for continuum mechanics for stress-strain relations of a linear elastic material.Boundary conditions are also taken from Hickey and Gottsmann (2014): the Earth's surface is treated as a free surface, the bottom boundary is fixed and the lateral boundary has a roller condition (free lateral, but no vertical displacement).We then adapt this model set-up for our purposes by adding a shallow, rectangular, poroelastic aquifer, which is saturated with water.The internal boundary conditions bordering the aquifer domain are (a) no flow and (b) continuous stress and displacement.Note that changing the lateral aquifer boundary condition to a fixed pressure instead of "no flow", hence allowing water to leave or enter the domain, does not change the results of this study.The initial pore pressure is set as hydrostatic: p f = ρ f gz, with z being the depth coordinate.The duration of the time-dependent simulation is 1000 days.We solve the full set of coupled equations, giving solid displacement u and fluid pore pressure p f .To demonstrate its meaning for water table changes that could be observed during volcanic unrest, we present all results as hydraulic head h, which is proportional to pore pressure: It represents the maximum water level change in a small diameter well (ideally a piezometer) in a confined aquifer; its initial value is 0 m.As a well moves with the ground, we subtract the vertical ground displacement from this hydraulic head change prior to using the data in order to obtain the relative water level change that would be measured in a well.Note here that, when "initial" responses are shown, the time referred to is 10 −6 s.The final model set-up is shown in Fig. 1; reference values of geometric parameters can be found in Table 2.
The linear elastic material surrounding the magma chamber, from here on called "host rock", has elastic properties of a general granitic crust.Depth-and temperature-dependent changes of Young's modulus of the crust are ignored for simplicity.We chose a Young's modulus of 30 GPa, representing an average value for crustal rocks in arcs for depths up to 8 km as derived from seismic velocity data (Gottsmann and Odbert, 2014).We test for two typical aquifer types found in volcanic regions: unconsolidated pyroclastic deposits, commonly composed of coarse ash to fine lapilli sized clasts, and vesicular basaltic lava flows.Vesicular here means a sufficiently connected porosity of the lava to serve as an aquifer; note that while the petrological porosity might be higher, only the connected pores matter for the fluid flow.These two types differ substantially in their elastic and fluid flow properties, which have significant influence on the observed signals.The layer above the aquifer, from here on called "cap rock", has elastic properties of a soft, impermeable clay.Input material properties for the reference simulation are given in Table 2; we used medians of parameter ranges found in the literature (Freeze and Cherry, 1979;Fetter, 1994;Wang, 2000;Gercek, 2007;Gudmundsson, 2011;Adam and Otheim, 2013;Geotechdata.info, 2013).Note that elastic properties of poroelastic layers are always required to be the drained parameters (i.e.measured under constant pore pressure).However, very few data exist on poroelastic parameters so we used the dry Young's moduli and Poisson's ratios instead and increased respective ranges in parametric sweeps to account for this unknown error.Some experimental data on the elastic properties of porous volcanic rocks are provided by Heap et al. (2014) for permeable tuff in the Neapolitan area; they fall within the here explored range of parameters.Within the different layers, material properties are considered isotropic and homogeneous.Standard water parameters are also given in Table 2.

Parametric studies and sensitivity analysis
In the parametric studies we investigated the effects of magmatic source properties as well as poroelastic and geometric properties of the aquifer (Table 2).Ranges for material properties of the aquifer were taken from literature (Freeze and Cherry, 1979;Fetter, 1994;Wang, 2000;Gercek, 2007;Gudmundsson, 2011;Adam and Otheim, 2013;Geotechdata.info, 2013).Ranges for geometric parameters can of course never cover the whole natural variation, we attempted to cover a reasonable range to be able to make general statements about the influence of certain parameters.A  This strains the surrounding linear elastic host rock (granitic crust), the poroelastic aquifer and the overlying linear elastic cap rock (clay).The water-saturated aquifer is modelled as either a vesicular lava flow or unconsolidated pyroclasts.An aquifer not covering the chamber but starting at some lateral distance L is realized by setting the darker grey region impermeable.The bottom boundary is fixed, the upper boundary is treated as a free surface, the lateral boundaries have a roller condition.There is no flow outside the aquifer; stress and displacement at the internal boundaries are continuous.An extract of the finite element mesh is shown only for illustration.The mesh density is finer around the cavity, at aquifer boundaries and the free surface.
non-geometric source property explored is its pressurization strength, for which we investigated a range of 3 orders of magnitude in an attempt to cover a sufficient range to recognise general patterns.As mentioned above, instantaneous pressurization is assumed in the simulations for simplicity, however chambers in reality will more likely pressurise over longer time periods.We therefore compare hydraulic head changes produced by the reference simulation with the poroelastic response to a chamber that inflates over 100 days.
When sweeping over one parameter, all others are kept constant.This entails that in all geometric sweeps, the distance between magma chamber top and aquifer was fixed, except for the sweep over magma chamber depth because this distance is such an important parameter it would have otherwise overwhelmed the pure effects of, for example, aquifer thickness.When investigating the effects of magma chamber shape, we changed the vertical semi-axis b of an ellipsoidal chamber, which then defines the horizontal semi-axis via the constant chamber volume.Pore fluid (H 2 O) temperature was effectively changed by varying its density, viscosity and compressibility.We used the program provided by Verma (2003) to calculate these parameters for varying temperatures and a pressure of 4.5 MPa, which represents average lithostatic pressure in the aquifer that was kept at 200 m depth (Table 3).In a subset of simulations, the central portion of the aquifer is replaced with an area of zero permeability out to a radial distance L but with the same poroelastic properties as the aquifer, to avoid numerical errors at the inner boundary of the aquifer.Aquifer density and porosity have a negligible influence and have not been included in the parametric study results.
To investigate the importance of parameters on hydraulic head change, we performed a sensitivity analysis.The influence of lateral distance L between magma chamber and aquifer onset has not been included in this analysis due to the lack of comparable signals (i.e.comparing the central hydraulic head change is not possible as the aquifer only starts at some radial distance) -it will be discussed in detail later.To assess the sensitivity of head changes to changes in individual parameter values, we compare the range of relative hydraulic head change ( h h ref ) produced by varying the respective parameter.To account for the change in time and space, this has been done for different locations in the model domain and at different times during the simulation.

Reference simulation
The described model was run for each aquifer type, using reference values of parameters given in Table 2, with a magma chamber pressurization of 10 MPa.In both aquifer types, the pressurization of the magma chamber induces a fall in hydraulic head, which is strongest directly above the magma chamber and decreases with radial distance from the chamber (Fig. 2).At distances larger than 5 km from the axis of symmetry -hereafter termed "centre" -the initial head change is of opposite sign compared to the central areas, but the amplitude of the head rise here is small in comparison with the central signal.There is also a comparatively small vertical gradient in the hydraulic head values.Whilst the pattern of the head change is the same in both aquifers, the absolute value of the signal differs substantially.In the pyroclastic aquifer, the maximum head fall is about 1.4 cm, while the hydraulic head in the lava flow aquifer falls by a maximum of 6 m.The initial hydraulic head change profile perfectly mirrors the strain curves (Fig. 2), illustrating that strain is the driver for the head changes.The aquifer is subject to dilation (positive strain), with a maximum value centrally above the chamber, which changes to compression (negative strain) with radial distance.Like for the hydraulic head, the two aquifers show similar patterns in strain, but different absolute values.Maximum volumetric strain in the pyroclastic aquifer is about 3 microstrain, while it is 13 microstrain in the lava flow aquifer.
Figure 3 illustrates the fluid flow pattern in the simulations, which are very different for the two different aquifer types.In the pyroclastic aquifer, fluid flow is away from the centre, while flow in the lava flow aquifer is towards the cen- tre.Despite its lower permeability, flow speeds are higher by more than a factor of two in the lava flow aquifer -as a result of the larger pressure gradient.Fluid flow is important because it equilibrates pressure in the aquifer and is responsible for the changes of strain and hydraulic head signals with time.Figure 3 shows the change with time of hydraulic head and volumetric strain, respectively, in a point in the aquifer centrally above the chamber.As water flows away from the centre, hydraulic head continues to fall in the pyroclastic aquifer until it reaches an equilibrium value of about −4 cm.Volumetric strain decreases and changes sign to compression after about 10 days; this compression increases and reaches an equilibrium value of about 10 microstrain.In the lava flow aquifer, hydraulic head increases, also tending towards an equilibrium value of about −4 cm with time; volumetric strain increases and also evolves to an equilibrium value.Whilst time-dependent changes take place almost until the end of simulation duration (1000 days) in the pyroclastic aquifer, the values in the lava flow aquifer reach equilibrium after less than 10 days.4).For specification of symbols please refer to Table 1.

Sensitivity analysis
draulic head change to the different parameters.For both aquifers, plots as shown in Fig. 4 were produced for four different locations in the aquifer, respectively.As the influence of many parameters varies in time and space, the ranking of parameters according to their significance is a two-step procedure.First, three parameter groups A, B, and C are defined based on the influence of a parameter on hydraulic head change in one location at one of the three tested simulation times (Table 4 and Fig. 4).Group A has the strongest influence on hydraulic head change, C the least.Parameters are then ranked into four priority groups based on the number of tested occasions (locations and simulation times) in which they belong to a certain parameter group: -priority 1 -parameters that belong to group A in ≥ 85 % of tests; -priority 2 -parameters that belong to group A in > 50 % of tests; -priority 3 -parameters that belong to groups B or C in ≥ 50 % of tests, but belong to group C in < 85 % of tests; -priority 4 -parameters that belong to group C in ≥ 85 % of tests.
Following the definitions of the priority groups we can rank the investigated parameters as follows: -priority 1 -pressurization value, volume and aspect ratio of the chamber; -priority 3 -chamber depth, aquifer depth and Young's modulus of the aquifer; -priority 4 -Poisson's ratio and thickness of the aquifer.
This ranking is, however, only a relative one -even those parameters of the last priority group have a non-negligible influence on the resulting hydraulic head change.Furthermore, the ranking of a parameter depends partly on the range of values tested for that parameter.This is particularly important in interpreting the sensitivity to the Biot-Willis coefficient (α).Due to the scarcity of experimental data for this parameter, the sweeps in both aquifer cases were performed almost over the whole mathematical range of α between the porosity and 1.However, the true value of α for natural soft rocks should be close to 1, while it is close to the porosity for hard rocks.Therefore, although ranked here as priority 2, in reality the Biot-Willis coefficient might belong in a lower priority group.More information on the individual influence of α can be found in Appendix A.

Results of parametric studies
The parametric sweeps provided a number of interesting insights; we are focusing here on describing the most important ones.

Influence of material properties
Of the aquifer's elastic properties, namely the Poisson's ratio ν and the Young's modulus E aq , only the latter is significant for the poroelastic response to applied strain.The most important hydraulic property is the permeability κ. Figure 5 shows the influence of these two important material properties on the hydraulic head response and its change with time.The initial hydraulic head change in an aquifer is identical for different permeabilities, but can be changed by orders of magnitude by changing the aquifer stiffness.For values between 1 and 10 000 MPa, a higher Young's modulus, i.e. a stiffer aquifer, leads to a larger hydraulic head response.However, this relationship is not monotonous as the head response decreases again when increasing E aq from 10 to 100 GPa (Fig. 5b).Note that, in the pyroclastic aquifer, a change in flow behaviour can be seen when increasing the stiffness: as long as E aq is smaller or equal to 10 MPa, hydraulic head falls with time as it does in the reference simulation.When the stiffness is increased to 100 MPa, the hydraulic head behaviour with time is comparable to that in the lava flow aquifer, where flow is towards the centre of the domain; hence hydraulic head increases with time (Fig. 5a).The permeability determines porous flow velocity and hence how fast the hydraulic head signal changes with time.A larger κ leads to a quicker change in head, visible for example in Fig. 5a for E aq = 100 MPa: while the hydraulic head signal after 1 day is still at the initial value for a permeability of 10 −14 m 2 , it is decreased by half for a permeability of 10 −10 m 2 .
The non-monotonous influence of the Young's modulus stems from the fact that not only its absolute, but also its value relative to the surrounding lithology is important.We therefore also performed parametric sweeps over the Young's Moduli of the host and cap rock, E h and E c , respectively, while E aq is kept constant.The ratios of elastic properties of the three rock layers are defined as and Figure 6a and b show the hydraulic head change depending on ER c for different ER h values.A larger ER h for a fixed E aq indicates a stiffer host rock and results in a smaller strain in the aquifer and hence smaller hydraulic head change.The relative cap rock stiffness ER c has negligible influence when it is small (generally: ER c < 0.1; for small ER h : ER c < 0.01).However, it becomes increasingly important when the cap rock stiffness is close to or larger than that of the aquifer (ER c > 1): a stiff cap rock can decrease the hydraulic head change and even change its sign to a head rise.This "signflipped" signal increases with larger ER c , as can be seen for the pyroclastic aquifer for ER c values larger than 100 (Fig. 6a).For small ER h values, this effect of ER c already takes place at ER c values larger than 0.1, visible for the lava flow aquifer (Fig. 6b).When sweeping over the Young's modulus of the aquifer, these elastic ratios are changed as well and contribute to the resulting head change which leads to the non-monotonous influence of E aq .
By sweeping ER c of the pyroclastic aquifer together with sweeping other parameters, we found that the ER c value at which the strain sign is flipped ("ER c flip ") is determined by the geometry of the system, in particular by the distance between aquifer and magma chamber and the thickness of the cap rock (Fig. 6c and d).The shorter the chamber-aquifer distance and the thicker the cap rock layer, the smaller is ER c flip .
The elastic stratigraphy determines the strain distribution in the domain, visible in Fig. 7a that shows vertical strain profiles from 2 km depth to the surface for different settings.Strain in the host rock is larger in the pyroclastic reference case -despite the fact that the same host rock is used in the reference simulations, the Young's modulus of the aquifer hence also influences the strain in the underlying rock.The graph also illustrates how strain changes at the boundary between different elastic mediums: strain increases when hitting a stiffer medium and vice versa.ER h determines this change at the host rock-aquifer boundary.In the reference settings, ER h is larger for the pyroclastic aquifer, leading to a smaller strain compared to the lava flow aquifer.The relative cap rock stiffness ER c then determines the strain change at the aquifer-cap rock boundary, but it also influences the strain change at the aquifer-host rock boundary.Figure 7a shows the strain profile for a simulation with a sufficiently large ER c to flip the sign of the signal.This change in sign is due to the strain jump at the host rock-aquifer boundary, where the dilatational strain in the host rock is turned into compression in the aquifer in the sign-flipped case.
Figure 7b shows the hydraulic head change along a horizontal profile in an aquifer with a sign-flipped signal.In contrast to the reference case, the central head change here is positive and changes sign twice: at about 3 km radial distance to a fall, and again at about 6 km to a head rise -mirroring sign-flipped volumetric strain.
It is common that aquifers are heated in volcanic settings.Figure 8 shows the substantial influence of changing the pore fluid temperature on the initial hydraulic head change and its evolution with time, especially when temperatures are above the pressure-dependent boiling point and the aquifer pores are no longer filled with liquid water, but steam.With increasing temperature of liquid water, the initial hydraulic head change is reduced in the lava flow aquifer, whereas it is very slightly increased in the pyroclastic aquifer (Fig. 8).For steam-filled pores (above 300 • C), the initial central hydraulic head change in the pyroclastic aquifer is 1 order of magnitude larger than for liquid water -in the lava flow Solid Earth, 6, 1207-1229, 2015 aquifer the opposite relation is true.Interestingly, the order of magnitude of hydraulic head change is the same in the two different aquifer types when steam saturated, while there is a 2 order of magnitude difference in the signals for the water saturated aquifers.In both aquifers, hydraulic head change increases with increasing temperature of the steam.Flow velocities are up to 1 order of magnitude faster in steam saturated pyroclastic aquifers when compared to water aquifers -lava flow aquifers show 1 order of magnitude higher flow velocities when saturated with water.Additionally, a change in flow behaviour is visible when conditions in the pyroclastic aquifer are changed: in a steam aquifer, hydraulic head increases with time, while it falls in the water-aquifer.

Influence of the geometry
To demonstrate combined geometric effects we plot the central initial hydraulic head change vs. the distance between magma chamber and aquifer for different chamber radii and absolute chamber depths in Fig. 9.The larger the chamber radius the larger is the resulting hydraulic head change in the aquifers.For the lava flow aquifer, the absolute magma chamber depth has no influence on hydraulic head change as long as the distance between aquifer and magma chamber is constant.The smaller this distance, the larger is the corresponding hydraulic head fall.This relation is somewhat more complicated for the pyroclastic aquifers, where hydraulic head depends on the distance between chamber and aquifer but also on the chamber depth.The central hydraulic head in the pyroclastic aquifers (Fig. 9a) is positive (hence sign-flipped) for a sufficiently small distance between aquifer and chamber, then switches sign to a head fall at a value dist flip , which increases with further increasing distance.The value dist flip is larger, the deeper the chamber.For a constant distance, hydraulic head change is larger for deeper chambers, if hydraulic head change is positive.Extrapolating the curves shows that if the hydraulic head change is negative, it is larger for shallower chambers.
Normally, the maximum hydraulic head fall is directly above the chamber.However, when considering a hydraulic head profile through the pyroclastic aquifer for a shallow magma chamber without a sign-flipped strain (e.g.z CH = 1 km and z aq = 200 m, see Appendix B), the maximum head change is no longer central but laterally offset by up to 1 km.
We also evaluated the influence of the shape of the magma chamber by incorporating tests for a prolate and oblate spheroid.Although chamber volumes are constant, the shape can change the hydraulic head signal by 1 order of magnitude.Figure 10 shows that the amplitude is highest for oblate chambers, intermediate for a sphere and smallest for prolate chambers.
Instead of having an "infinite" aquifer covering the whole volcano, we also varied the lateral distance L between the centre of the model and the onset of the aquifer, realized by a zero permeability zone in the centre of the domain (compare Fig. 1).The initial hydraulic head in these shorter aquifers equals the respective value at the same location in the reference aquifer.After some time, however, the head signal in the shorter aquifer differs from the reference case.Figure 11 shows the head changes after 10 days of simulation in the pyroclastic aquifer and after 1 day in the lava flow aquifer, respectively (the different timescales were used to account for the faster processes in the latter case).Compared to the profile of hydraulic head in the reference simulation, the maximum hydraulic head fall after this time (at the respective locations) in aquifers starting at 2 km radial distance is about 50 % larger in the pyroclastic aquifer and about 50 % smaller in the lava flow aquifer.This difference is strongest close to the lateral aquifer boundary facing the domain centre, where head falls are largest -with radial distance, the head change profile of the shorter aquifer approximates the reference profile.In the pyroclastic aquifer the difference between reference and shorter aquifer is negligible after the first kilometre, while in the lava flow aquifers head values differ considerably from each other over longer distances.
For larger values of L, the pyroclastic aquifers also differ in a comparable manner from the reference case, but become indistinguishable on the centimetre scale at distances larger than 6 km.Lava flow aquifers starting at 4 km or further radial distance are all significantly different from the reference case as they show positive head changes (up to 10 cm), while the reference aquifer at this time shows negative values everywhere less than about 10 km radial distance from the centre.
Figure 12 shows the different flow patterns in the first 8 km of three different aquifers after 0.1 days of simulation.Radial distance (km) Radial distance (km)

Influence of source pressurization
The initial hydraulic head response linearly depends on pressurization strength of the source.We assumed instantaneous pressurization for simplicity; however, real magma chambers will more likely inflate over longer time periods.In Fig. 13 we compare central hydraulic head evolution due to instantaneous pressurization with hydraulic head evolution due to a pressurization that is stepped up over 100 days (reaching the same maximum pressurization value).In this simulation, hydraulic head in the pyroclastic aquifer decreases more or less in parallel to the increase of chamber pressurization, until it reaches its maximum fall after 100 days.The maximum fall and its further development with time is not distinguishable from the hydraulic head evolution in the reference simulation at the respective time.Hydraulic head change in the lava flow aquifer at first also falls more or less in parallel with the increase of chamber pressurization, but reaches its maximum earlier at about 50 days, when pressurization slows down (due to the definition of the step function; compare inlet in Fig. 13a).This maximum fall of about 8 cm is significantly smaller than the maximum (initial) fall in the reference simulation of about 6 m.Note though that the hydraulic head in the reference simulation is already equilibrated to −4 cm at 50 days.After reaching its maximum, the hydraulic head fall decreases and reaches approximately the same equilibrium value as the reference simulation when stepping up of pressurization is complete.

Model limitations
In order to investigate poroelastic aquifer responses to crustal deformation, we made some simplifying assumptions.For one, the presented models only consider singlephase, single-component flow under constant temperature conditions.However, our parametric studies have shown that the pore fluid properties significantly influence the resulting head changes.Hydrothermal systems can contain steam, water and a number of solutes, and temperatures can change substantially.This can also affect the solid matrix, as its mechanical behaviour may deviate from elastic when it is sufficiently heated.Additionally, the injection of hydrothermal fluids into the aquifer can lead to a pore pressure increase, heating and further deformation (see e.g.Fournier and Chardot (2012) for a one-way-coupled model).We focused on the pure poroelastic response, but the poroelastic, heating and phase change processes superimpose.
Secondly, the aquifer was fully saturated and confined.To keep this study feasible, we did not investigate unconfined aquifers as this would imply a non-saturated permeable zone, and the coupling of linear elastic behaviour with non-saturated porous flow is associated with a high computational effort and often the solvers fail to converge due to the high nonlinearity of the problem.The model also ignores any hydrological sources and sinks, such as meteoric recharge, which can significantly influence well level observations in reality.
The discussed models are most applicable to confined aquifers that do not undergo extensive heating during the observation period (e.g.aquifers at some distance of the volcanic centre).They present a good opportunity to better understand poroelastic aquifer responses that have been used for monitoring.Their advantage over previous models is the full two-way coupling of flow and linear elastic behaviour and that we are able to simulate various geometries.The comparatively short computation time (on the order of 10 to 15 min per simulation depending on geometric complexity and number of necessary time steps) allows the study of a large number of parameters and their influence on hydraulic head changes and flow pattern.

General aspects
Our simulations show that neither injections of fluids nor flow within a fracture -hereafter termed "fracture flow"is needed to induce hydraulic head changes of several metres in an aquifer.Volumetric strain induced by a quasiinstantaneous magma chamber pressurization causes immediate hydraulic head changes in local aquifers.Dilation above the chamber, due to ground uplift, leads to a fall in pore pressure, while the accompanying compression at more than 5 km distance from the centre of the uplift induces a head rise.Poroelastic processes are therefore a reasonable approach to interpret rapid and large water level changes observed at volcanoes and they should not be ignored when studying hydrological systems in volcanic areas.For the same source and model geometry we could observe large differences between the two typical aquifer types.These differences are mainly due to the different elastic properties of the aquifers: the pyroclastic aquifer is much softer than the lava flow aquifer and therefore strain attenuation is stronger, hence the resulting hydraulic head change is smaller.
Both the strain-induced pressure gradient in the aquifers and the topographic gradient due to the ground uplift induce porous flow; groundwater flows from larger to smaller pressure/hydraulic head and from higher to lower elevation (Eq.5).The chamber pressurization in our reference simulation leads to a central ground uplift of about 4 cm, leading to a topographic gradient that is opposing the pressure gradient induced by the strain: the topographic gradient points away from the centre, while the pressure gradient points towards the centre.As the hydraulic head change in the pyroclastic aquifer is smaller than this uplift, gravitational flow dominates over strain-induced flow, and hence fluid flow is away from the centre.In the lava flow aquifer, strain-induced hydraulic head changes are two orders of magnitude larger than the topographic change; therefore flow in this aquifer is towards the centre, following the dominant pressure gradient.Despite its lower permeability, the lava flow aquifer shows higher flow velocities, as the strain-induced pressure gradient in this stiffer aquifer is large enough to overcome not only the topographic change but also the difference in permeability.
Fluid flow leads to the changes of strain and hydraulic head signals with time.Hydraulic head continues to fall in the pyroclastic aquifer as water flows away from the centre, while the opposite flow direction in the lava flow aquifer leads to a decrease of the initial hydraulic head fall with time (i.e. head increases).In both aquifers the equilibrium hydraulic head is approximately balancing the change in elevation (about 4 cm in the centre), as this is the equilibrated water level signal that would be recorded in the well that moved upwards with the ground.When the Young's modulus of the pyroclastic aquifer is sufficiently high (Fig. 5), the poroelastic response of this aquifer also becomes large enough to overcome the topographic gradient; hence the flow pattern changes to flow towards the volcano, causing hydraulic head in the centre to increase with time.These considerations show how the entire process is mostly governed by the elastic properties of the aquifer, not its hydraulic properties.
Strain changes simultaneously with hydraulic head due to the poroelastic nature of the aquifers.As water flows away, the pyroclastic aquifer responds to the removal of pore fluid with compaction -explaining the change of strain from dilation to compaction.The volumetric strain increase in the lava flow aquifer stems from the initial stress absorption by the pore fluid (final term in Eq. 2), which manifests as the pore pressure change.With equilibration of the pressure in the aquifer this stress absorption effect vanishes and strain approximates an equilibrium value that represents the strain value in an elastically equivalent, but dry material.Here, stress absorption of the fluid leads to an initial strain reduction by about 15 %; this value can be increased by increasing the Biot-Willis coefficient.
As pointed out by, for example, Rice and Cleary (1976), the initial elastic response (i.e.instantaneous deformation) of a porous medium can be calculated with the pure elastic solution by using undrained elastic parameters in the stress-strain relations.Indeed, the initial strain results of our poroelastic simulations agree with solutions calculated with just the structural mechanics module applying undrained parameters, which can be derived from the drained parameters with the Biot-Willis coefficient.However, the equilibrium strain signal in the aquifers does not correspond to the elastic solution calculated with drained parameters, except when one ignores gravitational flow, as flow down the slopes of the uplifted volcano adds strain changes that are not taken into account by purely elastic formulations.
The above findings highlight the necessity of a full coupling of fluid and solid mechanics.Both the effect of ground deformation on the pore fluid, as well as the influence of a pore fluid on strain in the solid matrix need to be considered to fully understand well level and/or strain signals.
Parametric studies have shown that poroelastic aquifer responses are complex processes that are strongly influenced by source, geometrical and aquifer parameters as well as the elastic stratigraphy.Chamber radius and pressurization determine the strength of the deformation source and the subsurface strain it causes.Strain partitioning in the crust is regulated by the elastic properties of the different layers; both the absolute and relative elastic properties of the aquifer and its surrounding lithology have a complex influence on the strain and head signals.A special case occurs when the cap rock is sufficiently stronger than the aquifer.A stiff cap rock prevents the dilation of the aquifer and turns the strain into compression, hence causing sign-flipped signals.In the reference set-up, the cap rock needs to be 2 orders of magnitude stiffer than the aquifer; this could be fulfilled if an unconsolidated, permeable pyroclastic layer is overlain by a lava flow.But other, perhaps more common, geological settings exist in which a sign-flipped response can be expected, as the geometry plays an important role as well: the thicker the cap rock, the smaller is the necessary ratio of cap rock to aquifer stiffness to change the sign of strain.For example, a cap rock that is only 3 times stiffer than the aquifer can already lead to a sign-flip if the aquifer is about 1 km deep.For our reference pyroclastic aquifer, this means a cap layer with a stiffness of more than 30 MPa, such as another pyroclastic layer that is slightly stiffer (e.g.due to a different lithology, rate of consolidation or grain size distribution) -a situation particularly feasible at stratovolcanoes.
The subsurface stress and strain fields are also substantially dependent on the shape of the chamber.For oblate chambers, the aquifer area that is exposed to vertical stress is larger than for prolate chambers and it is therefore subject to stronger strain.Additionally, the centre of the oblate chambers is shallower than the centre of the prolate chambers (as the depth of the chamber top is fixed in the simulations).The distance between aquifer and magma chamber is another factor contributing to the strength of the strain field affecting the aquifer.Generally, the closer the aquifer to the source, the stronger the strain and hence its pressure response.However, if elastic properties are close to values causing a sign-flipped signal, i.e. if the aquifer is rather soft, a sufficiently close aquifer-source distance can lead to a sign-flipped strain (because ER c flip is changed) and hence sign-flipped hydraulic head change in the aquifer.
The elastic properties of the solid matrix as well as the pore fluid together with the Biot-Willis coefficient of the aquifer determine the initial pressure response of the aquifer to the strain.Permeability then determines the velocity of pressure equilibration and gravitational flow and therefore the development of head and strain signals with time.Of particular interest is the influence of pore fluid temperature.It can change the hydraulic head response by 1 order of magnitude as well as influence the flow behaviour in the aquifer.This is especially important in volcanic environments, where heat flow is high and therefore temperature changes are likely.Changing the temperature means changing compressibility, density and viscosity of the water.We attempted to distinguish their individual influence with simulations in which only one of the three parameters was changed to a value corresponding to steam, while the others were kept at values corresponding to water.The viscosity does not influence the initial head fall, but affects the speed of equilibration, which is slower for a liquid water viscosity than for a lower steam viscosity (compare Eq. 5).Changing the phase of the pore fluid does not have this straightforward effect however, as flow velocities are also determined by the initial pressure gradients.These are influenced by fluid density and compressibility: decreasing the density of the pore fluid increases the initial hydraulic head change, while increasing the compressibility decreases it.At higher temperatures high enough for a phase change from liquid to steam, fluid density is reduced, while its compressibility rises.We therefore see a complex combination of these two effects.In the pyroclastic aquifer, the density effect dominates (hence hydraulic head change is larger), while in the lava flow aquifer the compressibility effect is more important (hence hydraulic head change is smaller).The stronger initial hydraulic head fall in the steam saturated pyroclastic aquifer is then large enough to overcome the topographic gradient, such that flow is towards the volcano.Therefore, as opposed to the reference simulation, the initial fall in hydraulic head diminishes with time in the soft aquifer just as it does in the stiffer aquifers.Hydraulic head falls in the lava flow and pyroclastic aquifers are of the same order of magnitude when saturated with steam, suggesting that the elastic properties of the solid matrix are less important and the processes are now governed by the fluid properties.While we investigated the temperature effect, other processes could also change pore fluid properties, such as dissolved minerals, and thereby play a role in determining the hydraulic head change.
Porous flow in the lava flow aquifer and therefore evolution of signals with time is also significantly influenced by the lateral distance between the magma chamber and the aquifer, even though initial hydraulic head values at respective locations are the same.For L = 2 km, i.e. in an aquifer that starts at 2 km radial distance from the centre, hydraulic head is up to 50 % smaller than the reference value (where L = 0 km).This effect is due to locally faster flow in the shorter aquifer.The positive hydraulic head changes in the lava flow aquifers starting at distances larger than 4 km are due to very different flow processes caused by changed initial pressure gradients.Generally, the initial hydraulic head is negative in the first 5 km of lateral distance and is positive at locations further from the centre; this pressure profile causes flow towards the centre in the lava flow aquifer.However, if the aquifer onset falls in the "positive head area", the driving pressure gradient and hence the flow direction are reversed.For lava flow aquifers that onset near the transition zone, two flow directions can be observed -one towards and one away from the centre of the volcano.These two directions are due to the maximum initial head change being not directly at, but lateral offset from the lateral aquifer boundary.This can be seen in the profile of the initial hydraulic head change in the reference simulations: beyond 5 km from the centre, hydraulic head change first increases with distance before decreasing again.While this comparatively small gradient is negligible in the reference simulation, as flow is dominated by the much stronger gradient towards the centre, it matters for flow in aquifers that start close to the transition zone from positive to negative strain -at least early in the simulation.
Horizontal flow directions in the softer pyroclastic aquifers are not affected by the changed initial strain and head gradient for different L, as flow is dominated by the topographic gradient, which still points away from the volcano.However, the vertical component of flow changes.In the reference simulation, the vertical gradient in hydraulic head change is overwhelmed by the horizontal gradient; hence lateral flow is dominant.In the shorter aquifers (i.e. with larger L values), the vertical gradients become more important and lead to up-or downward components of flow, respectively.That hydraulic head changes are up to 50 % larger close to the lateral boundary of the shorter pyroclastic aquifers than in the reference aquifer is due to the difference in mechanical properties of the central, impermeable portion and the outer permeable portion of the aquifer.While initially elastically equivalent, the change of head and strain with time due to porous flow in the outer aquifer leads to a mechanical boundary at the lateral aquifer onset.Especially in the pyroclastic aquifer, where strain undergoes significant flowinduced changes, this discontinuity in strain, head and flow is hard to simulate (compare Fig. C1a) and hence numerical signals close to this boundary should be interpreted with caution.The strain discontinuity is negligible in the lava flow aquifer, where strain does not change much with flow (compare Fig. C1b).
We only briefly studied the effect of long-term inflation, but results show that the time scale of pressurization is non-negligible as flow processes act simultaneously with the response to increased pressurization and can significantly change the signals.While hydraulic head responses in soft aquifers, where flow mostly follows the topographic changes, are comparable to instantaneous pressurization, the hydraulic head signals in stiffer, strain-dominated aquifers are reduced as the flow quickly equilibrates strain-induced pressure changes.Flow works against the increased pressurization and the rate of change of inflation determines which effect dominates, i.e. whether hydraulic head continues to fall (pressurization dominates) or starts to reach its equilibrium value (flow dominates).

General considerations
We have shown that wells can reflect the deformation at volcanoes, suggesting that their implementation in volcano monitoring systems could provide insights into subsurface processes causing the strain.However, prior to the interpretation of well signals, one needs to carefully consider that water levels can also be changed by several other processes, e.g.meteorological influences (rainfall), hydrothermal fluid injection, heat transferred conductively through the crust or changes in flow conditions due to the opening or closure of fractures.These processes can also act simultaneously and overcome hydraulic head changes caused by a poroelastic response.Under certain circumstances the different processes can be distinguished.First, the general hydrological behaviour -i.e. the meteorological responses -and up-to-date meteorological information should be tracked and therefore be reasonably well known if wells are to be included in a monitoring system.Then, a water level response to strain will be a transient signal on top of the background behaviour.Well level monitoring can form an important component for volcano monitoring in conjunction with geophysical or geochemical observations to track magma reservoir evolution.For example, ground deformation data will be useful for identifying hydrothermal injections.When hydrothermal fluids from a magma reservoir are injected into surrounding rocks, hydraulic head in the hydrothermal systems will rise and the ground will be uplifted (e.g.Todesco et al., 2004;Hurwitz et al., 2007;Rinaldi et al., 2010;Fournier and Chardot, 2012), while strain due to chamber inflation generally leads to a water level fall together with ground uplift.Temperature sensors in monitoring wells provide additional valuable information, as hydrothermal fluid injections should lead to relatively fast temperature effects compared to heat transferred from an inflating chamber or a dyke intruding through the crust, which would reach the aquifer much later -in fact, after the poroelastic response is already equilibrated.More complicated is the recognition of fracture flow.A hydraulic head change caused by rapid volumetric straining may be distinguished from flow-induced changes on the basis of timescales of the changes: the pore pressure response to strain is instantaneous, flow processes are slower.However, a gradual evolution of volumetric strain (see Fig. 13) will result in a slow hydraulic head response.Again, the combination with other monitoring data will help constrain the strain rate to distinguish between reservoir induced level changes from those by fluid flow.
Our parametric studies show how poroelastic aquifer responses are influenced by a variety of source, geometrical and aquifer parameters, which each have the potential to significantly alter the signal amplitude and development with time and space making the poroelastic processes highly complex.Consequently, a change in any of these parameters could lead to a change in an observed hydraulic head.In addition, the porous flow alters the initial hydraulic head signal with time.Therefore, not all observed aquifer pressure transients are necessarily related to a change in the magmatic system, which needs to be carefully considered when interpreting observed water level changes.
Yet another limitation lies in the fact that chamber inflation generally is not instantaneous.As the focus of this study was to identify the different influences of model parameters on hydraulic head signals, we assumed instantaneous pressurization for simplicity.However, it is clear that the analysis of signals during long-term inflation is even more complicated as one needs to decipher flow and inflation effects -again emphasizing that well data should ideally only be used in conjunction with surface deformation data.

Strain sensitivity and fluid flow
If the level changes are thought to be caused by strain, our models suggest that volumetric strain in the aquifer can be directly inferred from measured water level changes, as the simulated initial hydraulic head change perfectly mirrors the strain.This requires a known strain sensitivity, the change of hydraulic head in the aquifer in metres per unit applied strain, which can be assessed by tracking water level changes as a result of predictable excitations.Figure 14 shows the theoretical strain sensitivity of the two aquifers used in the reference simulation, determined by dividing the simulated hydraulic head change by the volumetric strain.This has been done along profiles through the aquifers.The very small strain and head values close to the transition zone from dilatational to compressional strain lead to numerical errors in the determined strain sensitivities in these locations (which can be reduced by increasing the mesh density), but in general we calculate a consistent value.Strain sensitivity of the pyroclastic aquifer is about −5 × 10 3 m; the lava flow aquifer has a strain sensitivity about 2 orders of magnitude larger.
However, the influence of flow on strain sensitivity is problematic; Fig. 15 shows how the theoretically calculated strain sensitivity changes with time: in the pyroclastic aquifer it quickly increases to about −250 × 10 3 m before it changes sign to a value of about 100 × 10 3 m, followed by a decrease -due to the change of strain from dilation to compression associated with the removal of pore fluid.In the lava flow aquifer, strain sensitivity shows a decrease (approaching  zero), which is comparable to the decrease of the hydraulic head change.Hence, the strain sensitivity value determined from aquifer responses to known strains only provides accurate strains when applied to the initial hydraulic head change, as it does not take flow and resulting poroelastic processes into account, and dense time series of well data (and ideally simultaneously recorded ground deformation data) are necessary to catch this response.
The better the local hydrology is known, the more value lies in well monitoring.Our simulations show how different soft and stiff aquifers behave in a strain field.Hence, knowledge of the lithology and/or determination of the strain sensitivity is important to discriminate between aquifer types.A low sensitivity value would indicate an aquifer similar to the presented pyroclastic example, which entails crucial information: the topographic gradient due to ground deformation can easily dominate over a strain-induced pressure gradient, and if this is the case, dilatational strain will quickly change to compression.Additionally, softer aquifers are more prone to the sign-flip effect.
Information about flow in the aquifer is important and the acquisition of permeability data, e.g.via pumping tests, help decipher flow processes.While these tests usually provide only a local value, high-resolution time series of head data can in fact provide information on the effective permeability of the aquifers.If a poroelastic head response has been identified, one can observe the equilibration of the pressure change with time -giving information on groundwater flow velocities.In any case, observing the flow behaviour in local aquifers, by installing several observation wells, is a valuable addition to existing monitoring efforts as they can reveal flow patterns caused by head changes, be they strain-induced or caused by other (volcanic) processes.Finally, it is important to know aquifer geometry as the models show that the flow pattern can strongly depend on the lateral distance from the aquifer to the source.

Application of analytical volcano deformation models
Even if strain sensitivity has been accurately used to infer volumetric strain, we still face the problem of interpretation of this signal.To invert for the source of volumetric strain, analytical volcano deformation models can be applied.However, these models commonly assume a source in a homogeneous half space and some only consider spherical or pointlike chambers (e.g.Dzurisin, 2007, and references therein).
Additionally, all previous approaches to model the deformation due to reservoir inflation treated their data as a result of dry deformation.Our results underline that any model using these simplifications will likely be inadequate when used for interpretation of poroelastic processes.Firstly, the assumption of a spherical chamber is a likely source for substantial mistakes and several chamber shapes have to be tested, as chamber shape is one of the two most important parameters influencing the signal.But even if chamber shape is taken into account, the assumptions of a homogeneous half space and dry deformation can lead to further misinterpretation.The stress absorption of a pore fluid leads to a reduction of initial strain in the aquifer when compared to an elastically equivalent dry layer.If the initial strain is used to infer the magmatic source based on a model for dry deformation, its strength can therefore be underestimated."Dry" strain is reached in the lava flow aquifer after porous flow has equilibrated the strain-induced pressure gradient.So, this problem could be solved when sufficiently dense time series of hydraulic head data are available: strain sensitivity can be combined with the evolution of signals with time to infer initial as well as equilibrium "dry" strain.In aquifers that are dominated by the topographic gradient, the presence of the pore fluid even leads to a reversal of dilatational to compressional strain and the application of dry deformation models is not possible.
The third assumption of a homogeneous half space is precarious as volcanoes are strongly heterogeneous -several previous studies have already shown that mechanical heterogeneities in the subsurface affect the ground deformation at volcanoes (e.g.Folch and Gottsmann, 2006;Manconi et al., 2007;Geyer and Gottsmann, 2010).Our investigation of the elastic stratigraphy has shown that hydraulic head change and consequently derived strain can also significantly deviate from signals in a homogeneous crust.Especially in settings with a sign-flipped signal, i.e.where the dilatational strain in the aquifer is turned into compression by a stiff cap rock, this influence becomes crucial.The hydraulic head rises and hence interpretation of the hydraulic head data alone would suggest a deflating chamber, while it is really inflating.We simulated ground deformation signals with the aim of investigating whether they -if available -could aid with this problem.Figure 16 shows the central hydraulic head and surface displacement signals with time for the pyroclastic aquifer with two different cap rocks: a soft one (E c = 70 MPa) and a stiff one (E c = 10 GPa).The latter leads to a change of sign of the volumetric strain.The surface deformation however does not change sign and shows inflation of the ground in both cases and can hence be used to indicate that the strain in the aquifer is sign-flipped.
The above considerations hint that the apparent inconsistency of observed well data and model predictions in the 2000 Usu case (Matsumoto et al., 2002) probably stems from the simplifications used in the applied Mogi model, which assumes a spherical pressure source in a homogeneous, dry half space.
In summary, while water level data can be a valuable addition to monitoring systems and give indications on subsurface strain, one needs to be careful when interpreting the head as well as strain data.We need to take into account that many parameters influence water level changes and that most of the commonly used analytical dry deformation models might fail to explain them.

Conclusions
In this study we presented fully coupled numerical models to investigate the interaction between solid mechanics and fluid flow in porous media.We have shown that strain due to the inflation of a magma chamber leads to significant hydraulic head changes and porous flow in the local hydrology.The flexibility of the finite element analysis method allowed us to perform extensive parametric studies providing detailed insights in these poroelastic processes.Parameters controlling aquifer behaviour are in order of importance (i) the shape, volume and pressurization strength of the magma chamber (ii) the phase of the pore fluid and the permeability of the aquifer (iii) chamber and aquifer depths and the aquifer's Young's modulus.Magmatic source properties and the distance between chamber and aquifer determine the strain field; strain partitioning is defined by the elastic stratigraphy of the crust.Elastic and flow parameters of the aquifer define its response to this strain and how head and strain signals change with time due to porous flow.
One aim of this study was to investigate the accuracy of the method to combine strain sensitivities with deformation models to interpret observed hydraulic head changes.Our models show that volumetric strain in the aquifer can indeed be inferred from hydraulic head changes using strain sensitivities, under certain conditions.Firstly, other causes for hydraulic head change have to be excluded, ideally by consulting other monitoring systems.Dense time series of well level data need to be acquired in order to account for flow processes and to measure the initial hydraulic head change.Additionally, we need to ensure that strain sensitivities have been accurately determined and have not changed with time due to changes in the hydrology.
However, using common analytical deformation models for the interpretation of this strain information is problematic, as several assumptions of these models can lead to substantial misinterpretation.They are only applicable for a comparatively homogeneous crust (i.e.E aq ≈ E c ≈ E h ), when one either accounts for fluid-induced strain reduction or considers an aquifer with very little strain reduction.The shape of the chamber needs to be taken into account as well.
The hydraulic head signal is very sensitive to source volume, shape and pressurization value.This suggests that if we have a detailed knowledge on the hydrology, some information about the source can be gained from hydraulic head changes -although solutions will always be non-unique.Our analysis has shown the necessity of numerical models to account for the large number of parameters that significantly influence the results.Nevertheless, well water levels and groundwater flow reflect subsurface strain and therefore are a valuable complement to other monitoring systems.So, in the pyroclastic aquifer for α ≤ 0.9, φχ f is 1 order of magnitude smaller than the right summand, which therefore dominates the definition of S. For α approaching 1, both terms become important.For the lava flow aquifer, E ≈ 10 10 and therefore the right summand has the order of magnitude 10 −11 × (1 − α) and φχ f is the dominating term in the definition of S for all α.Therefore, in the pyroclastic aquifer, changing α changes the coupling terms and the specific storage, while in the lava flow aquifer changing α has almost no effect on S.

Appendix B: Extra information on the influence of chamber depth
For a shallow magma chamber, in a situation with no signflipped strain, the maximum head change in the pyroclastic aquifer is no longer central, but laterally offset by up to 1 km (shown for z CH = 1 km and z aq = 200 m in Fig. B1).

Figure 1
Figure1.2-D axisymmetric model set-up: a boundary load P is applied on a cavity at depth, with the radius r for the spherical case or vertical semi-axis b for the ellipsoidal case, respectively.This strains the surrounding linear elastic host rock (granitic crust), the poroelastic aquifer and the overlying linear elastic cap rock (clay).The water-saturated aquifer is modelled as either a vesicular lava flow or unconsolidated pyroclasts.An aquifer not covering the chamber but starting at some lateral distance L is realized by setting the darker grey region impermeable.The bottom boundary is fixed, the upper boundary is treated as a free surface, the lateral boundaries have a roller condition.There is no flow outside the aquifer; stress and displacement at the internal boundaries are continuous.An extract of the finite element mesh is shown only for illustration.The mesh density is finer around the cavity, at aquifer boundaries and the free surface.

Figure 2 .
Figure 2. Results of the reference simulation, shown as the initial (i.e. 10 −6 s) hydraulic head change (blue, solid line) and volumetric strain (red, dashed line) along profiles through the two aquifer types due to a magma chamber pressurization of 10 MPa.(a) pyroclastic aquifer, (b) lava flow aquifer.Both aquifers show a fall in hydraulic head mirroring the dilatational strain curves.

Figure 4 Figure 3 .Figure 4 .
Figure4shows the range of relative hydraulic head change ( h h ref ) produced by varying individual parameters, demonstrating the significant differences in the sensitivity of hy-

Figure 5 .
Figure 5. Influence of Young's modulus and permeability of the aquifer on the central hydraulic head change and its evolution with time.(a) Pyroclastic aquifer, (b) lava flow aquifer.

Figure 6 .
Figure 6.Influence of the elastic stratigraphy on the central, initial hydraulic head change, shown using the ratios of Young's Moduli ER c = E c E aq , and ER h = E h E aq .(a) Pyroclastic aquifer, (b) lava flow aquifer.Also shown is the behaviour of central initial head change in the pyroclastic aquifer with ER c for different distances between chamber and aquifer (c) and the behaviour of central initial head change in the pyroclastic aquifer with ER c for different cap rock thicknesses (d).

Figure 7 .Figure 8 .
Figure 7. (a)shows vertical strain profiles through the centre of the domain from 2 km depth to the surface, for different elastic stratigraphies: the two reference simulations and a setting in which the sign of strain (and consequently hydraulic head) is flipped from dilation to compression due to a sufficiently stiff cap rock (ER c = 1000).(b) shows hydraulic head change and volumetric strain along a horizontal profile through the aquifer in the reference pyroclastic aquifer in comparison to the aquifer, where sign of strain and head is flipped.

Figure 9 .Figure 10 .
Figure 9. Influence of the geometry, i.e. chamber radius, chamber depth z CH , and distance (dist) between aquifer and chamber, on the initial central hydraulic head change.(a) pyroclastic aquifer, (b) lava flow aquifer.

Figure 11 .
Figure 11.Influence of lateral distance L between chamber and aquifer on the hydraulic head change profile through the aquifer after t = 10 d (pyroclastic aquifer) and t = 1 d (lava flow aquifer), respectively.(a) pyroclastic aquifer, (b) lava flow aquifer.

Figure 12 .
Figure 12.Influence of lateral distance L between chamber and aquifer on the flow pattern in the aquifers, here shown for t = 0.1 d for the first 8 km of each aquifer.(a) pyroclastic aquifer, (b) lava flow aquifer.Arrows indicate flow direction at the point where the arrow is attached, their length is proportional to flow velocity (different scales for a and b), colours show velocity magnitudes.Note that arrows are upscaled for better visibility; therefore, they cross the upper aquifer boundary.However, this is only a plotting effect, there is no flow leaving the aquifer.

Figure 13 .
Figure 13.Comparison of the development of central hydraulic change with time due to an instantaneous chamber pressurization and a pressurization over 100 d (with the same maximum pressurization value).(a) Pyroclastic aquifer, (b) lava flow aquifer.Inlet in (a) shows how the pressurization is stepped up over time.

Figure 14 .
Figure 14.Strain sensitivity in the aquifers, determined by dividing simulated hydraulic head change by the volumetric strain, along a profile through the aquifers.(a) Pyroclastic aquifer, (b) lava flow aquifer.Very small strains close to the transition zone from dilatational to compressional strain lead to numerical errors (can be reduced with increasing mesh density).

Figure 15 .
Figure 15.Strain sensitivity in the aquifers, determined at a point centrally above the chamber for different simulation times.(a) pyroclastic aquifer, (b) lava flow aquifer.The value significantly changes with time, depending on flow processes in the aquifers.

Figure 16 .
Figure 16.Central vertical surface deformation and hydraulic head change with time for the pyroclastic aquifer overlain by different cap rocks, showing the effect of a sign-flipped strain in comparison to the reference case.

Figure B1 .Figure C1 .
Figure B1.Hydraulic head change profile in the pyroclastic aquifer for z CH = 1 km.

Table 2 .
Input parameters: reference values and ranges for parametric studies (where performed).

Table 4 .
Parameter groups definition for the ranking resulting from sensitivity analysis.