Numerical models for ground deformation and gravity changes during volcanic unrest: simulating the hydrothermal system dynamics of a restless caldera.

. Ground deformation and gravity changes in restless calderas during periods of unrest can signal an impending eruption and thus must be correctly interpreted for hazard evaluation. It is critical to differentiate variation of geophysical observables related to volume and pressure changes induced by magma migration from shallow hydrothermal activity associated with hot ﬂuids of magmatic origin rising from depth. In this paper we present a numerical model to evaluate the thermo-poroelastic response of the hydrothermal system in a caldera setting by simulating pore pressure and thermal expansion associated with deep injection of hot ﬂuids (water and carbon dioxide). Hydrothermal ﬂuid circulation is simulated using TOUGH2, a multicomponent multiphase simulator of ﬂuid ﬂows in porous media. Changes in pore pressure and temperature are then evaluated and fed into a thermo-poroelastic model (one-way coupling), which is based on a ﬁnite-difference numerical method designed for axi-symmetric

Abstract. Ground deformation and gravity changes in restless calderas during periods of unrest can signal an impending eruption and thus must be correctly interpreted for hazard evaluation. It is critical to differentiate variation of geophysical observables related to volume and pressure changes induced by magma migration from shallow hydrothermal activity associated with hot fluids of magmatic origin rising from depth. In this paper we present a numerical model to evaluate the thermo-poroelastic response of the hydrothermal system in a caldera setting by simulating pore pressure and thermal expansion associated with deep injection of hot fluids (water and carbon dioxide). Hydrothermal fluid circulation is simulated using TOUGH2, a multicomponent multiphase simulator of fluid flows in porous media. Changes in pore pressure and temperature are then evaluated and fed into a thermo-poroelastic model (one-way coupling), which is based on a finite-difference numerical method designed for axi-symmetric problems in unbounded domains.
Informed by constraints available for the Campi Flegrei caldera (Italy), a series of simulations assess the influence of fluid injection rates and mechanical properties on the hydrothermal system, uplift and gravity. Heterogeneities in hydrological and mechanical properties associated with the presence of ring faults are a key determinant of the fluid flow pattern and consequently the geophysical observables. Peaks (in absolute value) of uplift and gravity change profiles computed at the ground surface are located close to injection points (namely at the centre of the model and fault areas).
Temporal evolution of the ground deformation indicates that the contribution of thermal effects to the total uplift is almost negligible with respect to the pore pressure contribution during the first years of the unrest, but increases in time and becomes dominant after a long period of the simulation. After a transient increase over the first years of unrest, gravity changes become negative and decrease monotonically towards a steady-state value.
Since the physics of the investigated hydrothermal system is similar to any fluid-filled reservoir, such as oil fields or CO 2 reservoirs produced by sequestration, the generic formulation of the model will allow it to be employed in monitoring and interpretation of deformation and gravity data associated with other geophysical hazards that pose a risk to human activity.

Introduction
Variations in geophysical observables, such as ground deformation at active volcanoes, are useful indicators of subsurface mass and density changes and can be evaluated as precursory signals to an impending eruption via data modelling. For caldera volcanoes in particular, earlier models focused on explaining ground deformation by magma emplacement (Anderson, 1937;Mogi, 1958;Bonafede et al., 1986;Bianchi et al., 1987;De Natale et al., 1991). Beside this interpretation, more recently models also consider the perturbation of hydrothermal systems (by pore pressure changes, variations in gas saturation and thermal expansions) as a possible (additional) source of spatio-temporal variations in deformation and gravity signals (Casertano, 1976;Gottsmann et al., 2003Gottsmann et al., , 2006aTodesco et al., 2003Todesco et al., , 2010Chiodini et al., 2007;Hurwitz et al., 2007;Hutnak et al., 2009;Ingebritsen et al., 2010;Rinaldi et al., 2010Rinaldi et al., , 2011Troiano et al., 2011).
The origin of unrest activities is still under debate in many restless calderas (such as at the Campi Flegrei, Italy), although for pre-eruptive hazard assessment it is fundamental to disentangle the signals generated by hydrothermal perturbations (e.g. Todesco and Berrino, 2005;Hurwitz et al., 2007;Hutnak et al., 2009;Todesco et al., 2010;Rouwet et al., 2014) from those related to magma movement towards the surface (e.g. Amoruso et al., 2008Amoruso et al., , 2014bWoo and Kilburn, 2010;Trasatti et al., 2011). Few deformation models account for significant complexities such as heterogeneities in key hydrological and mechanical properties of matrix and faults, which might influence both the path of ascending magma and the sub-surface circulation of hydrothermal fluids. Here we present a numerical model to evaluate ground deformation and gravity changes caused by the hydrothermal fluid circulation in restless calderas, taking into account the abovementioned complexities.

Background and motivation
Although the model is applicable to any caldera system, the model parametrisation in this paper is based on data available from the Campi Flegrei (CF) caldera in Italy. The CF, situated to the west of Naples, formed as a result of two structural collapses associated with the eruptions of the Campanian Ignimbrite (39 ka) and the Neapolitan Yellow Tuff (14 ka) (Orsi et al., 1996;Rolandi et al., 2003;Deino et al., 2004). The CF has received growing attention from the scientific community due to its reawakening in the last 50 years after a period of quiescence since the last eruption in 1538 with background slow subsidence at a rate of ∼ 1.5 cm yr −1 (Parascondola, 1947). Renewed unrest was associated with two periods of bradyseism (1969-1972 and 1982-1984), with a total vertical deformation of about 3.5 m (Troise et al., 2008). To date these uplifts have not culminated in an eruption. After 1984 a period of more than 20 years of general subsidence followed, interrupted sporadically by a series of minor uplift events. Since 2006 the caldera started uplifting again with an increased rate from 2011 (De Martino et al., 2014). Maximum ground deformation is recorded near the town of Pozzuoli, while the main fumarolic activities occur ∼ 800 m away at La Solfatara.
Significant gravity changes associated with unrests are usually observed in caldera systems (Berrino et al., 1984;Battaglia et al., 2003;Gottsmann et al., 2003;Todesco and Berrino, 2005), either at the centre of maximum deformation or at the structural boundaries of the caldera complex, which are likely associated with caldera ring faults (e.g. Gottsmann et al., 2006a).
Ring faults significantly alter strain partitioning and fluid propagation and hence must be considered for the interpretation of geophysical signals (De Natale and Pingue, 1993;De Natale et al., 1997;Beauducel et al., 2004;Folch and Gottsmann, 2006;Troiano et al., 2011;Jasim et al., 2015). In this paper we explore the impact of vertical and lateral mechanical heterogeneities in the shallow crust beneath the CF, including ring faults, on monitoring signals at the surface (ground deformation and gravity changes) as a consequence of unrest caused by a perturbation of the shallow hydrothermal system. Unrest is modelled by the injection of a mixture of hot water and carbon dioxide at the centre of the caldera system, which is associated with the main fumarolic activity at La Solfatara, and at the base of the ring faults, which simulates fluid release from a deeper pressurised reservoir (Jasim et al., 2015). We investigate the separated contribution of pore pressure and thermal effects to total ground deformation through a series of generic test cases which compare the single (central) injection model with the simulation of multiple injection points. We then show that different injection rates alter the timescales and amplitudes of deformation and gravity changes during periods of unrest. A sensitivity analysis of fault mechanical properties is also provided.
It is important to note that, while models are informed by data on the solid and fluid mechanics of the CF, we do not attempt to replicate or fit observations made during the ongoing unrest at CF.

Model parametrisation
In order to account for the complex mechanical structure of the shallow crust and the caldera infill at a restless caldera (such as CF caldera), the modelling domain is subdivided into several regions with different hydrological and mechanical properties. The model is 2-D axi-symmetric and defined by the coordinates (r, z), with r the radial distance and z the vertical position. The hydrological model is 1.5 km deep and is closed to heat and fluid flow in the radial direction and to fluid flow across much of the basal boundary (Fig. 2, detailed in Sect. 3.1), whereas the mechanical model is unbounded in the radial and downward vertical direction (Fig. 3, detailed in Sect. 3.2). Both models are based on information available for the CF and designed such that a central fumarolic field is situated on its rotational axis.
Two high-angle faults (Faults A and B) are implemented with parameters informed by data on the ring faults of the Neapolitan Yellow Tuff (14 ka) and Campanian Ignimbrite (39 ka) eruptions, respectively (De Natale and Pingue, 1993;Orsi et al., 1996;Folch and Gottsmann, 2006;Piochi et al., 2014). The fault geometry is represented in Fig. 1 Orsi et al., 1996;Jasim et al., 2015). The fault zone is divided into two sub-zones with different hydrological and mechanical characteristics: a central narrow (25 m wide) core zone is bordered on both sides by a wider (100 m wide) damage zone, the latter having properties intermediate between those of the core and the rock surrounding the fault zone (Tables 1 and 2).  Table 1. Inclination and radial placement of faults are not in scale. The domain extends toward infinity in the radial and vertical (downward) directions. Free-stress boundary conditions are ascribed at the top boundary, while vanishing displacements are assigned at infinite distances.

Hydrothermal model
Simulation of the hydrothermal circulation is performed by the well-known TOUGH2 software, a fluid flow and heat transport simulator of multiphase multicomponent fluids in porous media accounting for phase changes, relative permeability of each phase and capillarity pressure (Pruess et al., 1999). TOUGH2 solves a system of mass and energy balance equations that can be summarised as follows (for a general case of a fluid with k components): where Q is the accumulation term, F the flux and q the source (or sink) term, while the subscript α = M i or E refers to the mass balance equation for the ith component or the energy balance equation, respectively. The accumulation terms and fluid fluxes (based on the extended Darcy law) for mass balance equations are where the subscript β = l or g refers to the liquid or gas phase respectively, φ is the porosity, ρ β the density, S β the saturation, χ i β the mass fraction of the ith component in the β phase, K and K rβ are the absolute and relative permeability,  Fig. 2: rock density ρ r (kg m −3 ), porosity φ, permeability K (m 2 ), thermal conductivity λ (W (m × K) −1 ), specific heat capacity C r (J (kg × K) −1 ). Matrix permeability is isotropic, but enhanced in the vertical direction k z by almost 2 orders of magnitude in the fault damage zone and by 3 orders of magnitude in the core of the faults. In other respects the fault zones have the same hydrological characteristics as the matrix (star symbol ).

Rock Porosity Permeability
Thermal Specific heat density conductivity capacity  respectively, µ β the viscosity, P β the fluid pressure andĝ the gravitational acceleration. For the energy balance equation, the accumulation term (Q E ) and the heat flux (F E ) are where e β and h β are the specific internal energy and enthalpy of the phase β, T is the temperature, and ρ r , C r and λ are the density, specific heat and the thermal conductivity of the rock respectively.
In this paper we simulate fluids of magmatic origin entering the domain as a mixture of two components (k = 2): hot water and carbon dioxide. This mixture is simulated by the EOS2 module of TOUGH2. The depth of the domain for the hydrological model is 1.5 km, since the focus is the shallow hydrothermal activity, maintaining temperature and pore pressure of the entire simulation within the range considered by TOUGH2-EOS2 equation of state modules (which does not extend to super-critical fluids).
Atmospheric boundary conditions (P = 0.101325 MPa and T = 20 • C) are prescribed on the top of the domain z = 0; lateral boundaries are assumed to be impervious and adiabatic. A heat flux of 0.195 W m −2 is assigned at the impervious bottom boundary during the entire simulation, specified in order to sustain a temperature gradient comparable to that estimated for CF -∼ 130 • C km −1 (Rosi and Sbrana, 1987;De Siena et al., 2010;Piochi et al., 2014).
Cell centres of the finite-volume mesh used in TOUGH2 for the 1.5 km depth domain are shown in Fig. 4a. Hydrological parameters (permeability, density and porosity) are obtained from averaging drilling data for AGIP's report (Piochi et al., 2014), while the thermal properties of the rocks (thermal conductivity and specific heat) are derived from Rosi and Sbrana (1987) and Todesco et al. (2010) (see Table 1).
Although all parameter values are specified according to measured data at CF, the rock permeability may vary over several orders of magnitude, and this variation may substantially influence the fluid flow and heat transport in all the simulations. Jasim et al. (2015) explore the sensitivity of the hydrological system to matrix (caldera fill) and fracture hydrological properties. However, exploration of a wide range of possible hydrological values goes behind the scope of this paper. Fault zones are assigned the hydrothermal properties of the surrounding rock, except for the permeability, which is represented by an anisotropic tensor K of Eq. (2): where k r and k z are the radial and vertical permeabilities, respectively. While k r equals the isotropic permeability of the surrounding rock (set at 5 × 10 −15 and 10 −15 m 2 for layers A and B), a higher value of k z is chosen for those cells of the TOUGH2 finite-volume mesh whose centre falls into the core (k z = 10 −12 m 2 ) and damage (k z = 10 −13 m −2 ) zones of the faults. In order to simulate the fumarole activities at the centre of the domain, a central conduit with a higher permeability is placed at the centre of the domain and represented by a vertical cylinder with a radius of 200 m. A transition zone is specified between this conduit and the bulk of the caldera fill which has intermediate hydrothermal properties, as in previous simulations of Todesco et al. (2010) and Jasim et al. (2015) (Table 1).

Geomechanical and gravity models
The elastic response of a porous medium to pore pressure and temperature changes associated with the circulation of hot fluids is modelled by linear thermo-poroelasticity theory. The thermo-poroelastic effects are taken into account by including the pore pressure and temperature terms in Hooke's law (Rice and Cleary, 1976;McTigue, 1986): where and σ are the strain and stress tensors, respectively, µ is the rigidity modulus, ν the Poisson's ratio, tr(σ ) = σ xx + σ yy + σ zz the trace of σ , I the identity tensor, P and T are pore pressure and temperature changes, respectively, K d is the bulk modulus in drained conditions, K s is the bulk modulus of the solid constituent (Rice and Cleary, 1976;Rinaldi et al., 2010), and β is the volumetric thermal expansion coefficient of the solid matrix. Since we assume that deformations occur slowly, the governing equations are represented by the equations of equilibrium ∇ × σ = 0 with σ obtained by the inversion of Eq. (4), leading to the following set of Cauchy-Navier equations (Fung, 1965): where α = 1 − K d /K s is the Biot-Willis coefficient and u is the deformation vector, and where we have used the relation K d = 2µ(1+ν) 3(1−2ν) . The third Eq. (5) represents the linear approximation of the strain-deformation relation for small deformations.

A. Coco et al.: Numerical models for ground deformation and gravity changes
Free-stress boundary conditions σ ·n = 0 are prescribed on the surface, where n is the outward unit vector orthogonal to the surface. Unlike the domain for the hydrothermal model ( Fig. 2), the computational domain of the problem defined by Eq. (5) is unbounded in the radial r and vertical z downward directions, and a vanishing displacement is assigned at infinite distance: lim r→∞ u = lim z→−∞ u = 0. Since we assume that the problem is axi-symmetric, we solve the 2-D axi-symmetric version of Eq. (5) in the unknown u = (u, v), where u and v are the radial and vertical deformation, respectively.
The unbounded domain is discretised by a quasi-uniform grid (Fig. 4b), whose resolution is finest close to the axis of symmetry and smoothly decreases toward infinity (Grosch and Orszag, 1977;Fazio and Jannelli, 2014). In this way artefacts introduced by artificial truncation of the domain are avoided. Equation (5) is discretised and solved by extending the finite-difference numerical method proposed by Coco et al. (2014) for Cauchy-Navier equations to thermo-poroelasticity equations.
Heterogeneities in mechanical properties (µ and ν) are taken into account. In particular, the rigidity modulus µ for each layer of Fig. 3 is derived from seismic p wave velocity v p data (Orsi et al., 1996;Zollo et al., 2008;Piochi et al., 2014) by the application of the formula of Mavko et al. (2009): Density values of the porous medium ρ are derived from V p by the Brocher equation (Brocher, 2005): An appropriate value of the Poisson ratio for volcanic regions of ν = 0.25 is specified for the whole domain, except in the damage and core zones of the fault areas, where higher values (0.30 and 0.40, respectively) are specified on the basis of the nature of the rock (Gercek, 2007). Rigidity values are reduced in the fault zones (µ = 0.385 and 0.0357 GPa for the fault core and damage zone, respectively, which correspond to E = 10 9 and 10 8 , where E is the Young modulus). The volumetric thermal expansion coefficient is β = 10 −5 K −1 after Rinaldi et al. (2010) and Todesco et al. (2010). All values are reported in Table 2.
In order to separate the contribution of pore pressure to the total ground deformation from thermal effects, we solve two different sets of differential equations for the mechanical simulation: where σ = λ tr( )I + 2µ is the elastic stress tensor (i.e. without taking into account pore pressure and temperature contributions). Let u T and u P be the solutions of the two problems Eq. (6), respectively. As a result of the linearity of the stress-strain relationship σ (u) and the divergence operator, the total ground deformation u can be expressed as the sum of the solutions to the two problems (namely u = u T + u P ). In practice, it is sufficient to solve only one of the problems Eq. (6) and obtain the other solution by difference.
Gravity changes g are computed by solving the following boundary value problem for the gravitational potential φ g (Currenti, 2014): where G is the gravitational constant and ρ is the density distribution change. The finite-difference method presented by Coco and Russo (2013) is applied to solve the problem Eq. (7) on an infinite domain, using the coordinate transformation method (Coco et al., 2014).

Numerical simulation scenarios and results
The background hydrothermal fluid circulation is driven by the injection of a mixture of hot water and carbon dioxide at a temperature of about 350 • C from the base of the central high-permeable conduit, simulating the input of fluids of magmatic origin. A heat flux is assigned at the bottom impervious boundary at a rate of 0.195 W m −2 . The steadystate solution, obtained after a long-lasting injection period (c. 4000 years), is used as the initial condition for the unrest simulations (run-up to a final time of 100 years), which are divided into the three scenarios described below.

Modelling scenarios
Scenario I: central injection at the base of the conduit (radius of 200 m) at the same temperature but at an increased rate with respect to that used for the steady-state quiescent solution (see Table 4); Scenario II -constant mass rate: Scenario I plus injection at the bases of each fault core zone of a total mass flow rate equal to that of the central injection (see Table 3); Scenario III -constant flux rate: Scenario I plus injection at the bases of each fault core zone at a specific (per square metre) mass flow rate equal to that of the central injection (see Table 3).
Injection at the base of the faults (core zone of Fig. 1, 25 m wide) for Scenarios II and III simulates the possible release of gas from a deeper reservoir ascending along preferential pathways of the fault zone during unrest periods.

Injection rates
Once the rates of the central injection are established, the corresponding injection rates at the base of the faults are de- Table 3. Different injection values (mass and flux rate) for the central conduit and faults, normalised to the injection of a mass of 1 kg of fluids. Base area of the central conduit is π × 200 2 = 125 664 m 2 , of the Fault A core zone is 2 π × 3000 × 25 = 471 239 m 2 , of the Fault B core zone is 2 π × 6500 × 25 = 1021 018 m 2 .
0 9.79 × 10 −7 7.96 × 10 −6 termined by Table 3. Rates of hot water and carbon dioxide central injection for both the steady-state and unrest simulations are selected in order to match observed data at CF, following other models present in the literature for simulating the unrest activity associated with the perturbation of the hydrothermal system (e.g. Chiodini et al., 2003;Todesco and Berrino, 2005;Rinaldi et al., 2010;Todesco et al., 2010). In particular, the injection rates for the steady-state simulation are chosen so that the total flux (3400 tons day −1 ) and the molar ratio CO 2 /H 2 O of 0.17 (equating to 1000 tons day −1 of CO 2 and 2400 tons day −1 of H 2 O) are based on average degassing measured prior to the 1982-1984 bradyseismic crisis, while an increased molar ratio of 0.40 is used for the unrest simulation. Regarding the magnitude of the injection rates, several values have been adopted in the literature in different contexts, albeit the rates used in Rinaldi et al. (2010) and Todesco et al. (2010) (6000 tons day −1 of CO 2 and 6100 tons day −1 of H 2 O) provide a good match to observed data. Recently, a constraint on the magnitude of the injection rates has been discussed by Afanasyev et al. (2015). Although there are many other parameters that can influence the mechanical response (including depth of injection and temperature of the injected fluid), in this paper we focus on the influence of the injection rates on the timescale and amplitude of the deformation (Table 4). Where not specified, the injection rates of the unrest × 1 column of Table 4 are used.

Initial conditions
Initial conditions for the unrest simulation are the same for all scenarios and represented in Fig. 5. Due to the injection of hot fluids, the central conduit shows a higher temperature with respect to the rest of the domain, while the pressure approaches hydrostatic, indicative of a steady-state condition. A slight temperature variation is observed at the fault zones, where the locally increased permeability focuses convective fluid flow, with downward flow of cold water via the fault (Jasim et al., 2015). A two-phase plume forms close to the central conduit, according to the results of previous fluid flow simulations Rinaldi et al., 2010;Todesco et al., 2010).

Pore pressure, temperature and density changes during unrest
At each time step of the unrest simulations we evaluate changes in pressure ( P = P − P 0 ), temperature ( T = T − T 0 ) and density ( ρ = ρ − ρ 0 ), relative to initial conditions (subscript 0) and use these to compute ground deformation and gravity changes at the surface by Eqs. (5) and (7). Density change is in practice computed as ρ = φ β (ρ β S β − ρ β0 S β0 ), where subscript β = l or g refers to the liquid or gas phase, respectively. We observe that ρ is mainly driven by the gas saturation change, since densities of liquid and gas do not significantly change during the simulation. For this reason we plot the gas saturation change S g = S g − S g0 rather than ρ (Figs. 5 and 6). Analysing Scenario I (Fig. 5), we observe that after 6 months of simulated unrest the zone of perturbed pore pressure has already approached the surface (z = −500 m) at the central conduit, with a maximum P of about 4 MPa observed at the injection point. Temperature and gas saturation changes remain small and confined to the areas surrounding the injection point. The maximum P of the entire simulation (about 5 Mpa) is observed at 3 years. At the same time, gas saturation changes reach the shallower layer (z = −400 m), while no changes in temperature are apparent. After 3 years P decreases; at 10 years hot fluid (warmed by up to T ∼ 100 • C) rises up to about z = −1000 m and the gas region extends up to the surface. At 100 years, which is the end of the simulation, P continues to decrease towards a new steady state, while T keeps increasing (with a maximum T ∼ 130 • C), extending the central plume laterally by up to 250 m. Gas saturation changes approach the steadystate solution, and a single-phase gas region is forming close to the surface. We do not observe any significant variation in pore pressure, temperature or density close to the faults, where the values remain the same as the initial condition.
The location of regions where significant changes in pore pressure, temperature and density are observed depends on the background simulation. During the steady-state simulation, fluids are injected only at the centre of the model, and thus a two-phase plume develops only in the central conduit,     Table 4) and additionally for Scenarios II and III injecting at the base of the core zone of the two faults according to Table 3. Note that the colour scale of initial conditions is different from the respective colour scale of unrest plots.
with complete liquid saturation within the fault zones. During the unrest the increased rate of injection at the conduit leads to an increase in pore pressure most markedly at depth within the conduit, but increases in temperature and gas saturation occur at the border of the expanding two-phase plume.
If we vary injection rates in Scenario I (Fig. 6), the amplitudes of P , T and ρ are strongly (nonlinearly) affected. Regardless of the injection rate, T continues to increase for the entire simulation (100 years), while P peaks at ∼ 3 years. Therefore, the timescale for pore pressure changes to reach the maximum value does not significantly depend on the injection rate. In particular, the maximum P is 2.15 for Unrest ×0.5, 9.85 for Unrest ×2 and 14.1 MPa for Unrest ×3 (all at t = 3 years). The maximum T is observed at the final simulation time (t = 100 years) and is 92.1 for Unrest ×0.5, 171 for Unrest ×2 and 181 • C for Unrest ×3. The extent of the central plume increases for the entire simulation: after t = 100 years the plume has extended laterally by up to 200 m for Unrest ×0.5, 450 for Unrest ×2 and 550 m for Unrest ×3.
In contrast to Scenario I, in Scenarios II and III injection at the base of the faults induces a perturbation in pore pressure, temperature and density at the fault zones (mainly located on the hanging wall), while the behaviour at the central conduit is similar in all three scenarios (Fig. 5). Due to the higher injection rates at the base of the faults, Scenario III shows more pronounced perturbations than Scenario II. Both faults behave similarly in Scenario III: the region with significant pore pressure change approaches the surface after 6 months (with a maximum P of about 2.5 MPa), while temperature

Ground deformation
At each time step of the unrest simulations, changes in pore pressure and temperature are interpolated from the finite-volume mesh of the hydrological model to the finitedifference mesh of the mechanical model (the two meshes are represented in Fig. 4) and fed into Eq. (5). This is known as one-way coupling between hydrological and mechanical models, as used previously by a number of studies (Hurwitz et al., 2007;Hutnak et al., 2009;Rinaldi et al., 2010;Tode-sco et al., 2010). It is a simplified approach compared with a fully coupled model that also takes into account the influence of stress and strain on permeability and porosity during the simulation (Neuzil, 2003;Rutqvist, 2011). In Scenario I (Fig. 7), for the first 10 years of unrest the uplift is maximum at the centre of the domain and decays radially. Vertical and horizontal displacements reflect the Mogi solution for a small spherical source (Mogi, 1958). The profile obtained at t = 100 years does not reflect a Mogi solution and presents a maximum total uplift of 21 cm at r = 300 m, decaying rapidly as radial distance increases. Temporal evolution of the ground deformation at the centre of the domain throughout 100 years of unrest (Fig. 8) indicates that the contribution of thermal effects (v T ) to the total ground deformation is almost negligible with respect to the pore pressure contribution (v P ) during the first years of the unrest, but increases in time and eventually becomes dominant. In particular, for lower injection rates (unrest ×0.5 of Table 4) the vertical deformation due to thermal effects only exceeds the pore pressure contribution after more than 100 years, while for higher injection rates (unrest ×2 and ×3 of Table 4) it takes less than 50 years. The amplitude of the deformation is nonlinearly dependent on the injection rate, while the timescale of the first local maximum is largely independent   Table 4). The solid line is the total vertical displacement v = v P + v T , while the dashed and dotted lines are the vertical displacement due to pore pressure v P and thermal effects v T , respectively. of injection rate, occurring after ∼ 3 years of unrest in all simulations. Vertical displacement due to pore pressure ef-fects (v P ) increases very rapidly with the onset of unrest. After this strong initial pressurisation (lasting about 3 years), vertical deformation starts decreasing towards a new steadystate value. Thermal effects strongly affect the long-term behaviour and their importance increases with increasing injection rates. Consequently, the timing of the local minimum, prior to the thermally induced later monotonic increase, occurs earlier for higher injection rates. Although we show only the temporal variation of the vertical deformation at the centre of the model for Scenario I, a similar pattern is observed localised around both faults for Scenarios II and III. In Scenario II (Fig. 9) the deformation profile reflects the injection of fluids at the fault zones. Maximum vertical deformation is observed at the centre of the model and two local maxima correspond to the faults (Fig. 9a). Magnitude of peak displacements both horizontal and vertical reduces from centre to Fault A and from Fault A to Fault B, reflecting the different injection rates.
After about t = 3 years the vertical deformation at the centre of the model reaches a temporary maximum (see solid line in Fig. 8 for Scenario I), then decreases toward a lower value (at about t = 10 years) while deformation on faults continues to increase. At t = 100 years the vertical displacement at the centre of the model increases again toward a steady-state solution (solid line of Fig. 8), while deformation on faults decreases toward a lower value. We observe in Fig. 9c Figure 9. Ground deformation computed for Scenario II at the surface after t = 0.5, 3, 10 and 100 years of unrest: total vertical deformation (a), total horizontal deformation (b), vertical deformation due to pore pressure (c), and vertical deformation due to thermal effects (d).
Vertical lines refer to the boundary of the central conduit and to the injection and shallowest points of faults.
vertical deformation profile at t = 100 years is almost exclusively attributable to thermal effects, which are negligible in the first years of the unrest simulation. Horizontal deformation (Fig. 9b) shows a Mogi-like pattern close to the central conduit (Mogi, 1958), while two peaks are observed close to the fault zones. For both peaks the deformation profile is steeper on the side towards the centre of the domain due to the fault inclination (dip-angle smaller than 90 • , Fig. 1), since the steeper deformation profile is always observed in the hanging wall. We finally observe for all the plots that the deformation profile is relatively smooth above Fault A, while there is a sharp kink above Fault B, because such fault reaches the surface (Fig. 2). Vertical deformation at the centre of the domain throughout the entire simulation (100 years) is practically the same as for Scenario I (Fig. 8), indicating that pore pressure and temperature changes along the faults do not significantly affect the mechanical behaviour of the fumarole.
In Scenario III (Fig. 10) vertical deformation on faults is greater than at the centre of the model (up to t = 10 years). Although pore pressure change at the faults shows a lower value compared with that close to the injection point, it is more vertically extensive (Fig. 5) due to the lower vertical permeability of the central conduit compared to the faults, causing a larger uplift. Vertical deformation at the axis of symmetry is also slightly amplified (by the mechanical influ-ence of faults) with respect to the one observed in Scenarios I and II.
Except for faults, the mechanical heterogeneities described so far depend only on depth, resulting in a 1-D heterogeneity structure. A complex mechanical structure for CF could be used, taking into account the lateral variation in mechanical properties to reflect differences between the two caldera infills, as proposed in the models of Trasatti et al. (2005), based on tomographic studies of Aster and Meyer (1988). Some simulations (not shown) have been performed with different matrix properties around faults, maintaining the same mechanical properties for fault core and damage zones. No significant differences were obtained close to fault areas, highlighting that the amount of deformation is mainly driven by the values of µ and ν assigned to the fault core and damage zones, especially when these values are much smaller than those assigned to the surrounding area (Table 2). A sensitivity analysis of the rigidity modulus on faults is provided below.

Sensitivity analysis on fault rigidity modulus
In this section we analyse the influence of rigidity of fault core and damage zone on ground deformation. For simplicity we restrict our analysis to the vertical component of deformation. In detail, µ c , µ d andμ are the rigidity values of the core zone, damage zone and the surrounding rock, respec- tively. We reduce the rigidity on the fault core (and damage) zone with respect to the surrounding rock by s 1 (and s 1 /2) orders of magnitude, i.e. µ c =μ 10 s 1 , µ d =μ 10 s 1 /2 .
In the baseline simulation the rigidity of the faults is the same as the surrounding area (i.e. s 1 = 0). For each value of s 1 in the range 0 ≤ s 1 ≤ 3 we obtain a variation in the ground deformation of s 2 orders of magnitude, i.e.
where v 0 is the uplift observed for the baseline simulation (i.e. s 1 = 0). Figure 11 shows the values of s 1 and s 2 computed at the centre of the model and at faults for simulation times t = 3 and t = 100 years. Reducing the rigidity values (i.e. increasing s 1 ), the deformation increases for the simulations at t = 3 years and decreases for t = 100 years. At t = 3 years the deformation is mainly driven by pore pressure changes (Figs. 8 and 10) close to injection points (therefore at a depth of ∼ 1.5 km), while at t = 100 years deformation is mainly driven by temperature changes, which constitute a shallow source of deformation (thermal effects reach the surface at t = 100 years, see Fig. 5). In the latter case, the region where the rigidity is reduced (fault core and damage zones) is below the source of deformation, causing less uplift than that observed for the baseline simulation. After t = 3 years sensitivity of deformation to fault rigidity is greater for Fault B than for Fault A, whilst the reverse is true at t = 100 years. Changes in deformation at the centre of the domain are minimal throughout all simulations, showing the limited lateral influence of the mechanical properties at the faults.

Gravity changes
The solution of Eq. (7) is the gravity change g = g − g 0 , where g 0 and g are the gravity distributions observed at the initial condition and at a fixed time of unrest, respectively. Evaluating g at a particular point of the surface (r, z = 0) means that also g and g 0 refer to the same geometric location (r, z = 0). Gravity change measured in the field g = g − g 0 is actually influenced by ground deformation, since g is measured at the same material point of g 0 , but at a different geometric (translated) point (r, z = 0) + u(r, z = 0), which takes into account the absolute movement of the gravimetry associated with the ground displacement. The value g = g −g 0 is often referred in literature as residual gravity (Bonafede and Mazzanti, 1998;Fernández et al., 2005;Gottsmann et al., 2006a), since it does not include the gravity change associated with the ground deformation (Telford and Sheriff, 1990). Gravity changes computed at the centre of the model (r = 0, z = 0) for different injection rates (Table 4)  in Fig. 12. After a transient increase (maximum 16.1 µGal for the Unrest ×1 model) over the first months of unrest (Fig. 12b), gravity changes become negative and decrease monotonically towards a steady-state value, although this is not reached within 100 years (Fig. 12a). The modulus of the gravity changes is more pronounced for higher injection rates, with a maximum increase after 0.5 years of 33.9 µGal for the Unrest ×3 model and a much larger negative value. The behaviour is, however, nonlinear at a fixed time with respect to injection rates, due to both the change of molar ratio from the steady state to the unrest phase and the nonlinearity of the hydrothermal model. The increase in injection rate causes only a minor increase in the time needed to change sign (from positive to negative, Fig. 12b). Figure 13 compares the gravity changes computed at the surface for different simulation times and three injection scenarios (D-I) and the associated vertical gravity gradient (A-C), computed as g/v, where v is the vertical deformation computed in Sect. 4.3. Again, this is usually referred as the residual gravity gradient, since it does not take into account the free-air correction (Gottsmann et al., 2006a). Data are plotted for up to 20 years of unrest, since after a long period of unrest the gravity gradient becomes unstable in most of the domain due to very small vertical deformation far from the faults and central conduit. Maximum values in modulus are observed at a radial distance of ∼ 570 m at the boundary of the two-phase plume, and are almost equal for the three sce-narios. However, local maxima of the modulus of the signals are present at the faults for Scenarios II and III. The absolute value is significantly higher for Scenario III, reflecting the higher mass flux.
The sign of the vertical gravity gradient is the same as that of the gravity changes, since the sign of ground deformation is almost always positive (i.e. uplift) in all the simulations. The pattern observed close to the axis of symmetry is similar to that for the gravity changes, presenting a local extreme at the border of the plume. In Scenarios II and III, the gravity gradient presents local extremes on the faults (most evident for t > 10 years) because of local extremes in both gravity changes and vertical deformation (see Appendix). In Scenario II the local extreme on Fault A is a minimum, since the wavelength of the gravity change profile on Fault A is lower than that of the vertical displacement, after a proper normalisation (see Appendix for more details). Local extreme on Fault B is a maximum, since g has a greater wavelength than v (see, for instance, the g and v profiles at t = 20 years in Fig. 13h). In Scenario III both extremes are minima, since the wavelength of the gravity change profile on the faults is lower than that of the vertical displacement, after a proper normalisation (see Appendix for more details). The value observed at the faults is much greater (due to greater gravity changes associated with greater injection rates). Heterogeneities in hydrological and mechanical properties as well as the presence of faults within caldera forming volcanoes in the model substantially affect the hydrothermal circulation of hot fluids and the consequent variation in geophysical signals. Models of the CF caldera suggest that the higher permeability of a central conduit at La Solfatara favours the uprising of hot fluids from the deep portion of the reservoir to the surface. Steady-state simulations show formation of a two-phase plume Rinaldi et al., 2010;Todesco et al., 2010), with radius and gas composition that depend on the permeability structure of the caldera fill . According to our simulations the two-phase plume occupies the entire central conduit and part of the transition zone, leading to a 300 m radius plume at 1.5 km depth. The radius of the plume reaches 500 m in a shallow region close to the surface. Two gas regions form at the edges of the plume: one surrounding the injection point and a shallower region which extends to the surface, simulating the gas discharging observed during the fumarolic activities at La Solfatara. The transition zone of intermediate hydrological properties favours pressurisation of the system during the first 3 years of the unrest and then allows depressurisation as injected fluids ascend and discharge at the surface (Fig. 5).
This behaviour is reflected by the fast initial vertical deformation at the centre of the domain, which is followed first by a rapid and then by a slower subsidence period (Fig. 8). This pattern would not be observed if the permeability contrast between the central conduit and the rest of the domain was stronger. The close relationship between deformation and fluid flow is highlighted in this simulation. If we lower the permeability of the caldera fill, subsidence after uplift does not occur. In fact lower-permeability caldera fill would inhibit the recharge of cold water to the base of the domain and the plume would be confined to a considerably narrower area, resulting in a hotter gas-saturated region, as shown in Todesco et al. (2010). Pressure release after the initial uplifting would not be present and the following period of subsidence would not be observed.
Although the deformation profile observed in Scenario I reflects the solution of a Mogi-type source (Mogi, 1958) in the first years of the unrest, over time it develops into a more complex pattern that cannot be explained by a simple deformation source (Fig. 7). In the long timescale the ground deformation is therefore mainly driven by the thermo-poroelastic response of the hydrothermal system.
Usually deformation observed at the centre of the model and associated with a central source located at the axis of symmetry is amplified by the mechanical heterogeneities of lateral fault zones (Folch and Gottsmann, 2006). This behaviour is not observed in the simulations of this paper, due to the small ratio between the central source depth (injection depth of 1.5 km) and radial distance of the closest fault (∼ 3 km), although in Scenario III the higher injection rate at the base of the faults gives a small amplification of deformation.
Rock expansion due to temperature changes is slower than that due to changes in pore pressure. Temperature changes are confined to the areas surrounding the injection points during the first 10 years of the unrest and take more than 50 years to reach the surface. Thermal contribution to the total ground deformation is therefore almost negligible within the first 10 years but becomes dominant after some decades of unrest (Fig. 8). Fournier and Chardot (2012) suggest that the relative contribution of temperature and pore pressure is directly proportional to the injection depth. Rinaldi et al. (2010) modelled the effect of a short unrest period (20 months) of high injection rate, and showed that the pore pressure declines immediately after cessation of fluid injection, while the temperature continues to increase until hot fluids discharge at the surface. Most recently, Chiodini et al. (2015) examined the accelerating rate of ground deformation affecting CF between 2005 and 2014, and suggested that the observed deformation pattern requires both an extended period of heating of the rock and short pulses of injection of magmatic fluids into the hydrothermal system.
In our simulations, maximum temperature change is located close to the edge of the plume (Fig. 5). Consequently, the maximum uplift observed at t = 100 years is slightly displaced from the centre. The shape of this temperature change is elongated in the vertical direction, resembling a prolate source, and causes the rapid decay of the vertical deformation. The same behaviour is observed for the gravity changes at the centre of the domain. Density changes are localised at the boundary of the plume, where replacement of water by gas over an increasingly large area occurs and gravity changes present a local extreme. Gas saturation changes are small during the first years of unrest and restricted to an area close to the injection point (Fig. 5). As a consequence, gravity changes take about 2 years to exceed 50 µGal in absolute value (for the Unrest ×1 case). Indeed, the initial period of the unrest is characterised by an increase in density, since a substantial amount of water is rapidly introduced to regions with positive gas saturation, following the increase in injection relative to the background rate. This perturbation is amplified for Unrest ×2 and ×3 models since a larger mass of water is injected, as inferred by the positive sign of gravity changes at the beginning of the unrest in Fig. 12. After a transient period this pattern is inverted, since the higher molar ratio of CO 2 /H 2 O of the fluid injected during unrest pushes the system toward a steady-state solution in which a substantial amount of gas will replace fluid-saturated regions, causing a negative change in density and consequently in gravity changes. In contrast, gravity changes over the fault zones are negative for the whole simulation time, since the base of the faults are liquid saturated at the beginning of the unrest (no background injection is performed at the base of the faults).

A. Coco et al.: Numerical models for ground deformation and gravity changes
The inclusion of faults in the model fundamentally alters the dynamics of fluid flow and heat transfer in the surrounding of fault areas. Jasim et al. (2015) show that the permeability contrast between the fault zone and surrounding rock affects local temperature gradients, causing faults to act as preferential pathways for either recharge or discharge of groundwater, depending both on fault/matrix permeability ratio and on the vertical extension of the fault. Temperature changes are more pronounced around the faults than at the central conduit, since the background hydrothermal circulation in the fault zones is not driven by any fluid injection, locally enhancing the contrast between the steady-state and unrest simulations. Gravity change and deformation associated with thermal effects are thus larger on the faults than close to the axis of symmetry.
Fault mechanical properties strongly influence the deformation profile in the vicinity of faults. In particular, a lower rigidity for the fault core and damage zones is associated with increased uplift on the fault where the source of deformation is deep (as in the case of pore pressure change during the first years of unrest, mainly localised around injection points) but with reduced uplift where the source of deformation is adjacent to the surface (as in the case of temperature changes after a long period of unrest). There is only minor perturbation of uplift observed at the centre of the domain, showing that mechanical properties of faults have a limited lateral influence. Such influence would be amplified if a deeper domain was considered (Folch and Gottsmann, 2006).
Fault geometry (inclination, vertical extension, penetration depth, radial distance, etc.) also influences the amplitude and pattern of deformation and gravity changes. Profiles of vertical deformation vary smoothly on Fault A, while a sharper contrast is present at the Fault B, likely because Fault B extends up to the surface z = 0. This sharp behaviour is mainly associated with the mechanical heterogeneities of fault core and damage zones rather than with hydrological causes (it would not be observed if the mechanical heterogeneities do not reach the surface).
Although the simulations performed in this paper provide a qualitative assessment of the contribution of hydrothermal fluid circulation at restless calderas, a more quantitative study and comparison with observed data from a particular caldera (such as the CF) is beyond the scope of this study.
It is important to consider limitations of the approach adopted in this paper. First, the shallow fluid injection (only 1.5 km deep) is constrained by the range allowed by TOUGH2 (which does not account for supercritical fluids), while several studies at CF have speculated that there is a deeper source, between 2.7 and 5 km (Gottsmann et al., 2006a, c;Amoruso et al., 2014a, b). Afanasyev et al. (2015) recently investigated deep supercritical regions of the hydrothermal system at CF using MUFITS, a multiphase multicomponent fluid flow in porous media simulator that accounts for high-temperature processes (Afanasyev, 2013a, b), more realistic for representing restless calderas.
In addition, whilst assuming that simple layering of rock properties is appropriate in the absence of detailed subsurface data, in reality it is probable that the stratigraphy of the caldera fill is more complex. Representing the effects of such heterogeneity, and in particular the strong local contrasts in the vicinity of the faults, is difficult using standard gridding approaches (Geiger and Matthäi, 2014). Small-scale geological heterogeneities observed in nature, usually modelled by geostatistical methods (Journel et al., 1998;Strebelle, 2002), cannot be correctly represented by a coarse cell blocks and identifying appropriate upscaling methods is challenging (Gerritsen and Durlofsky, 2005;King et al., 2006). On the other hand, using an extremely fine grid would radically increase the computational cost, making the model unusable for practical purposes where a number of simulation runs is required, such as optimisation and uncertainty reduction (Oliver and Chen, 2011).
The 2-D axi-symmetric representation of ring faults is obviously not able to describe the complex fault networks which characterise collapse calderas. For example, circulation along fault planes is a purely 3-D phenomenon that cannot be represented by a 2-D model. However, this study provides a first approximation of the influence of fluid flow mechanics around faults on relevant geophysical observations and indicates the importance of this area for future research.
Last but not least, the one-way coupling adopted in this paper, although provides a reasonable simplification for short period unrests, is not appropriate for the simulation of prolonged processes, since a significant variation in key hydrological parameters (permeability, porosity) can be induced by a change in stress and strain (Neuzil, 2003;Rutqvist et al., 2002), altering the long-term behaviour of fluid flow in the porous medium and the consequent evaluation of geophysical signals. For example, since an increase in the effective stress may cause a permeability and porosity reduction (Davies et al., 1999;Rutqvist et al., 2002), a drop in these hydrological parameters is expected where higher deformation are observed, namely at the centre of the domain and close to the fault zones. This may reduce the deformation and gravity change profiles over time. In addition, since these changes in permeability and porosity would be less pronounced where deformation is lower, the permeability contrast between the central conduit and the transition zones would be attenuated, modifying the dynamics of the rapid uplift and subsequent deflation observed in Fig. 8. However, a qualitatively analysis is difficult to perform at this stage for a number of uncertainties, such as the sensitivity to parameters regulating the relationship between effective stress and permability/porosity.

Conclusions
The model proposed in this paper is targeted at evaluating the variations in geophysical parameters associated with the perturbation of the hydrothermal system in a restless caldera. A correct evaluation is fundamental to discriminate between magmatic and hydrothermal unrest. Although the model can refer to a generic system, parameters have been chosen on the basis of the CF caldera, to simulate a behaviour proposed to explain the periodic unrests at the CF caldera since 1969. This periodic behaviour can be explained by a series of brief injections of hot fluids into the hydrothermal system Todesco et al., 2010) or after a long thermal process following an increase in rock heating, as highlighted by Chiodini et al. (2015). Similarly, Jasim et al. (2015) show that periodic behaviour of gas composition can be associated with sharp increase of the heat flux, with periodicity comparable to the decennial cycle observed at CF. Simulations performed in this paper evaluate the ground deformation and gravity changes caused by a long period of unrest associated with a prolonged injection of fluid of magmatic origin into the shallow hydrothermal system at a higher rate compared to that of the background simulation. To represent the inherent complexities at collapse calderas, we considered the effects of heterogeneities in the vertical and lateral distribution of hydrological and mechanical parameters and the effect of faults. Permeability contrasts considerably affect the fluid flow pattern Jasim et al., 2015) and consequently ground deformation and gravity changes at the surface.
The presence of the ring faults formed as a consequence of the episodes of collapse can significantly alter the behaviour of the system in the surrounding of the fault zones. Higher permeability parallel to the plane of the fault favours convection and upward discharge of hot fluids from depth, perturbing the hydrothermal system by changing pore pressure, temperature and fluid density, dependent on injection rate (compare Scenarios II and III). These perturbations, together with weaker mechanical properties of fault core and damage zones, substantially alter geophysical signals (ground deformation, gravity changes) at the surface close to the faults; furthermore, in Scenario III, a minor influence on the centre of the model is observed.
Investigation of different scenarios shows that the qualitative and quantitative perturbations of the fluid dynamics are sensitive to fluid injection rates, whose correct evaluation is one of the key challenges to improve the understanding of restless caldera systems.