A hydro-geophysical simulator for fluid and mechanical processes in volcanic areas.

Eﬃcient and accurate hydrothermal and mechanical mathematical models in porous media constitute a fundamental tool for improving the understanding of the subsurface dynamics in volcanic areas. We propose a ﬁnite-diﬀerence ghost-point method for the numerical solution of thermo-poroelastic and gravity change equations. The main aim of this work is to study how the thermo-poroelastic solutions vary in a realistic description of a speciﬁc volcanic region, focusing on the topography and the heterogeneous structure of Campi Flegrei (CF) caldera (Italy). Our numerical approach provides the opportunity to explore diﬀerent model conﬁgurations that cannot be taken into account using standard analytical models. Since the physics of the investigated hydrothermal system is similar to any saturated reservoir, such as oil ﬁelds or CO 2 reservoirs produced by sequestration, the model is generally applicable to the monitoring and interpretation of both deformation and gravity changes induced by other geophysical hazards that pose a risk to human activity.


Introduction
Volcano geophysics focuses on integrating data from different monitoring techniques and developing numerical and physical models to explain the observations. The geophysical observations collected in volcanic areas are the surface expressions of processes that occur deeply within the volcanic edifice. Integration of different geophysical observations and mathematical models enables to identify renewed volcanic activities, forecast eruptions, and assess related hazards []. Unrest periods are defined as variations in the geophysical and geochemical state of the volcanic system with respect to the background behaviour. Usually, geophysical changes observed during unrest periods are modelled in terms of volume and pressure changes in a magma chamber embedded within an elastic medium [-]. Indeed, this approach often appears at odds with the complex processes accompanying volcanic unrest. Particularly, the interplay between magma chambers and hydrothermal systems may result in the heating and pressurization of hydrothermal fluids, which in turn induce ground deformation and variations in the rock and fluid properties [-]. A quantitative evaluation of this interaction is fundamental for a correct hazard assessment in volcanic areas.
A thermo-poroelastic numerical model is here proposed to jointly evaluate ground deformation and gravity changes caused by hydrothermal fluid circulation in complex media with surface topography and mechanical heterogeneities in a D axis symmetric formulation.
Although the model is applicable to a generic volcanic system, in this paper we set up the model parametrisation on the Campi Flegrei caldera (CFc), a volcanic area situated to the west of Naples. After a long period of quiescence, two main uplift episodes (- and -) highlighted a reawakening phase of the CFc without culminating in an eruption []. The total amount of uplift of . m induced the evacuation of about , people from the town of Pozzuoli and surrounding areas. A slow subsidence followed the second uplift episode, periodically interrupted by mini-uplifts. Since  CFc started uplifting again, with a particularly increased rate from  [, ].
To reduce the uncertainty related to hydrothermal activities, significant complexities must be accounted when a mathematical model is adopted, such as real topography and mechanical heterogeneities, which are indeed neglected by most of existing models [, , , ]. In addition, highly efficient and accurate numerical solvers are required when a large number of simulation runs must be executed for a single simulation scenario, such as for optimization purposes related to geophysical data inversion or when the thermo-poroelastic response must be calculated at each time step in a coupled hydrothermal/mechanical model [].
The thermo-poroelastic numerical method proposed in this paper is second order accurate and based on a finite-difference ghost-point discretization for complex geometries successfully adopted in [] to solve elliptic equations and in [] for elasto-static D problems with plane-strain assumptions. Here, the methods proposed in [, ] are extended to account for thermal expansion and pore pressure effects (caused by the perturbation of the hydrothermal system) in a D axis symmetric framework. Although a D model would be suitable to represent a more accurate scenario, the axis symmetric assumption is a reasonable approximation for representing caldera systems and volcano edifices in general, which are usually characterized by radial structures.
In order to avoid artifacts introduced by finite truncation of the domain, a coordinate transformation method is adopted in order to prescribe vanishing solutions at infinite distances maintaining the second order accuracy []. An extension of the multigrid solver described in [, ] is adopted as well.
Hydrothermal activity is simulated by TOUGH, a well known multi-phase multicomponent software for fluid flow and heat transfer in porous media []. The hydrothermal system is perturbed by injecting fluids of magmatic origin at the base of a central conduit for a prolonged period (until steady-state conditions), simulating the fumarolic activity at La Solfatara []. Following unrest periods are modelled by an increased injection rate. Variations in relevant geophysical signals (pore pressure, temperature and density) between the initial condition (steady-state) and unrest period are computed at each time step and fed into the thermo-poroelastic model (one-way coupling) to compute associated ground deformations and gravity changes. A similar one-way approach has been implemented at CF in previous works (e.g., [, , ]).
The model proposed in this paper is widely applicable to a large number of relevant scientific and engineering problems for addressing subsurface flow and transport problems, such as geological carbon sequestration, nuclear waste disposal, energy production from geothermal, oil and gas reservoirs as well as methane hydrate deposits, environmental remediation, vadose zone hydrology. The coupling with the suite of TOUGH simulators makes the model suitable for simulating and interpreting geophysical changes induced by a wide range of alterations in the subsurface flows, such as CO  sequestration and geothermal exploitation. However, the main scope of this paper is related to the methodology rather than the application to case studies and comparison with real data, and then the focus is restricted to the simulation of ground deformation and gravity changes in a generic volcanic area whose model parametrisation is set up on the Campi Flegrei caldera.

Hydro-geophysical model
The mathematical model is based on the governing equations of the thermo-poroelasticity theory (Mechanical model, Section ..), which describes the elastic response of a porous medium to the propagation of hot fluid through pores (Hydrothermal model, Section .) by the effects of variations in pore pressure and temperature. Alterations in the gravity field associated to fluid density variations are evaluated by the model described in Section ...

.. Mechanical model
Assuming that the deformation induced by pore pressure and temperature changes occurs slowly compared to the time scales of elastic waves, and that the system behaves as elastic (neglecting then the elasto-plastic and viscoelastic effects, which, to a first approximation, is a reasonable assumption for shallow crusts), the rock can be considered in static equilibrium and the displacement is found by solving the equations of equilibrium coupled with thermo-poroelastic extension of Hooke's law [], giving the following set of equations []: where σ and ε are the stress and strain tensors, respectively, u is the deformation vector, λ and μ are the Lame's elastic medium parameters and I is the identity tensor. Eqs. () are solved in the domain M (the mechanical domain) for the unknown u. In our simulations M is an infinite domain in the radial and downward vertical directions, i.e. Figure ). A suitable set of boundary conditions is posed to close the mathematical problem, namely zero displacements at infinity and stress-free boundary conditions σ · n s =  on the ground surface M = {z = f TOP (x, y)}.
Analyzing the stress tensor (second equation of ()), we observe that two terms are added to the elastic stress tensor of the general Hooke's law: the P pore-pressure contribution from poroelasticity theory through the α BW Biot-Willis coefficient and the T temperature contribution from thermo-elasticity theory through the linear thermal expansion coefficient α. Both P and T refer to perturbations with respect to a previous equilibrium state.

.. Gravity model
Fluid circulation, temperature and pore-pressure changes necessarily alter the fluid density distribution, which, in turn, affects the gravity field. Gravity variations in active volcanic areas range in amplitude between a few μGal ( μGal =  - m/s  ) and several hundred μGal with a spectrum varying from seconds to years depending on the involved subsurface mass-redistribution processes. The gravity change g, arising from fluid density redistribution, can be calculated by solving the following boundary value problem (Poisson equation) for the gravitational potential φ g []: where G is the gravitational constant and ρ is the density distribution change. The problem is closed imposing the condition of vanishing gravitational potential at infinity. Density changes in the medium are computed with respect to an initial density distribution with φ being the porosity, S β and ρ β the saturation and density of the phase β, respectively.

Axis symmetric formulation
In order to reduce the computational burden associated with D models we reformulate problems () and () in cylindrical coordinates and assume that they are axis symmetric, reducing then to a D formulation ( Figure ). We recall the coordinate transformation expressing the cylindrical coordinates (r, θ , z) in terms of Cartesian coordinates (x, y, z): where r, θ and z are the radial, angular and vertical coordinates, respectively. For simplicity, we maintain the same nomenclature for all quantities, even if they are expressed in cylindrical coordinates, {σ , ε, u, λ, μ, P, T, ρ}(r, θ , z), f TOP (r, θ ).
The components of tensors σ , ε, u are named in such a way the cylindrical coordinate formulation is highlighted by a proper subscript, as follows.
Strain tensor ε components expressed in terms of u components read [, ]: Stress-strain relationship presents the same functional formulation as in Cartesian coordinates, namely the second Eq. of (): Using the cylindrical coordinates formulation of the divergence operator acting on a second order tensor field, Eq. () (∇ · σ = ) becomes [, ]:  r The axis symmetric formulation corresponds to assuming that: (i) all quantities are independent of the angular coordinate θ , namely ∂q/∂θ =  for q ∈ {σ ij , ε ij , u i , λ, μ, P, T, f TOP }, (then they can be expressed only in terms of the radial and vertical coordinates by q(r, z)); and (ii) the displacement along the angular coordinates is zero (i.e. u θ = ). By these assumptions, Eqs. () simplify to: Using () and (), the stress components become (observe that σ θθ is not necessarily zero): and finally Eqs. () reduce to a PDE system of two equations (since the first Eq. of () vanishes) in the  unknowns u r and u z : The D normal vector to an axis-symmetric surface M = {z = f TOP (r)} is n = (n r , , n z ), leading to the following axis symmetric formulation of the stress-free boundary conditions σ · n = : The D domain M and the surface M reduce to a D domain and a curve, respectively (that we continue to call M and M for simplicity), expressed as (see Figure ): In summary, Eqs. () (with σ ij expressed by Eqs. ()) in the domain () with stress-free boundary conditions () on () and zero displacement at infinity constitute the D axissymmetric thermo-poroelastic model. By a similar argument, the axis-symmetric formulation of Eq. () becomes:

Numerical method
Pressure, temperature and density changes induced by fluid circulation are computed using the TOUGH numerical code []. Starting from these quantities, ground deformation and gravity changes are solved through Equations () and () using a finite-difference ghost-point method developed in the context of elliptic problems [] and recently extended to solve elastostatic equations in D plane strain configuration for unbounded domains []. Here, we extend the methods developed in [, ] to solve the axis symmetric formulations of Section .. The method consists of three stages: firstly, the unbounded domain problem is transformed in a bounded one by the coordinate transformation method. Secondly, the new set of equations and boundary conditions in the bounded domain is discretized by a finitedifference ghost-point approach. Finally, the discrete equations are solved by a proper geometric multigrid technique.

.. Coordinate transformation method
By the coordinate transformation method, the unbounded problem with variables (r, z) is mapped into a bounded one with new variables (ξ , η) as follows: where χ : [-, ] → [-∞, +∞] is a differentiable strictly increasing function in ] -, [ such that χ() =  and lim s→± χ(s) = ±∞. Equations () and () (as well as the associated boundary conditions) are transformed according to (), namely following the transformation of differential operators: The computational domain C = [, ] × [-, ] is therefore discretized by a uniform Cartesian grid, which automatically results in a quasi-uniform grid for [, +∞]×[-∞, +∞] (see, for example, Figure ). Imposing vanishing solution at infinity is equivalent to require zero Dirichlet boundary condition on ∂C for the transformed problem. The choice of the mapping function χ is crucial to define the mesh distribution in such a way a loss of accuracy of the numerical method is avoided. This is accomplished by two requirements: (i) the finest region of the quasi-uniform grid must cover the areas where the sources are concentrated, and (ii) the asymptotic behavior of the solution must be correctly captured. To satisfy the requirement (i), we observe that in the case of Eqs. () and () sources are concentrated where P, T and ρ are different from zero. Since these quantities are evaluated by the hydrothermal model (detailed in Section .), this area is certainly contained in H . Therefore, an acceptable mesh distribution should guarantee to cover a large A level-set function φ is used to discriminate between internal points (φ < 0, represented by black dots) and external points (see (15)). Among the latter, grid points close to the surface (ghost points, represented by red circles) are included to the discretized problem and a suitable value is assigned to them by enforcing boundary conditions (by the ghost-point method described in Section 2.3.2). The coordinate transformation method is used also to solve Eq. (12) in the domain R 2 (then without a ghost-point level-set method).
region containing H by a fine grid. The resolution of the grid then decays going to infinity at the same order as the solution, to satisfy (ii). In practice, knowing the a priori asymptotic decay of solutions of Eqs. () and (), it is possible to calibrate the mapping function in order to satisfy (ii) (namely, the gradient of the solution of the bounded problems does not have to develop singularities at the boundary of [, ] × [-, ]). A good choice of mapping function is χ(s) = cs (-s  ) . We refer the reader to [] for an accurate discussion on the choice of a proper mapping function of elasto-static problems. The parameter c regulates the length scale of the computational grid (and then the maximum resolution). In practice, to ensure a fine resolution in H , we impose that the spatial resolution at the boundary of H is half the resolution at the origin (where there is the minimum grid spacing). We choose c so that this condition is satisfied. In detail, we impose χ (χ - (α)) =  · χ (), where α is the radius of the hydrothermal system (in our case α =  km). If m = , after some algebra we obtain the condition: By requirements (i) and (ii) the numerical method will achieve second order accuracy in a natural way, unlike the case when artificial truncation of the domain is used to assign zero displacement and zero potential boundary conditions.

.. Finite-difference ghost-point method
Observe that, if the topography function f TOP is known, a level-set function can defined by φ = zf TOP (r). In the applications, the topography is often defined by a set of grid points (ξ k , η k ) with k = , . . . , N T . A reasonably accurate level-set function can be obtained by φ(P) = ±d(P, t), namely the signed distance between the point P and the line t passing through the two topography points (ξ k , η k ) closest to P. The sign is determined by the mutual position between P and the two closest boundary points in such a way φ(P) <  if and only if P ∈ b M . Level-set methods are a powerful tool to implicitly describe moving geometries as well as complex topological changes []. Although in this paper we do not face moving geometries, the level-set function comes useful for two main reasons: (i) it provides accurate information, such as normal direction and curvature, that can be adopted in the discretization technique to improve the accuracy of the method, and (ii) it allows us to provide a general numerical method that can be embedded in a more generic framework of time-dependent problems with (complex) moving geometries. Reason (i) is strictly related to the discretization technique adopted, which is described below.
The level-set function is able to discriminate between internal (φ < ) and external (φ ≥ ) grid points. For each internal grid point Eqs. () are discretized using the standard central finite-difference discretization. In details, observe that the differential terms of Eqs. () (with σ components given by ()) are a linear combination of: where γ is a smooth coefficient and w = u r or u z . The finite-difference discretization of the previous terms is: where h = ξ = η is the (uniform) spatial step and γ i+/,j = (γ i,j + γ i+,j )/. Finitedifference discretization for an internal grid point involves then the surrounding grid nodes. Some of these surrounding points may lie outside of b M , defining the so-called set of ghost points. Therefore the list of unknowns of the discrete linear system consists of the discrete values {u P r , u P z }, where the subscript P varies within the set of internal and ghost nodes. To close the linear system (i.e. to have a number of equations as much as the number of unknowns) two discrete equations must be generated for each ghost point G (to equate the two unknowns u G r and u G z ). This discrete equations are recovered by prescribing the stress-free boundary condition on b M with second order accuracy using the level-set function. In practice, referring to Figure , the discretization technique associated to a ghost point G consists of the following three stages.
• First we compute the projection B of G to the boundary by the level-set function. In details, we compute the outward unit normal vector n G by discretizing Eq. () on G using a standard finite-difference scheme; then we apply the bisection method to solve the D nonlinear equation φ(B) =  on the segment between points G and G  = G -√  h n G , with h being the spatial step. • Then, we identify a nine-point stencil made by internal and ghost points, containing G and whose convex hull (the smallest convex set containing the stencil) contains B. whereσ rr ,σ rz ,σ zz andφ are the biquadratic interpolant on the nine-point stencil of σ rr , σ rz , σ zz and φ, respectively.

.. Multigrid solver
The discrete linear systems involving the unknowns {u P r , u P z }, with P varying within the sets of internal and ghost points, is solved by an efficient geometric multigrid solver, which is an extension of the multigrid approach proposed in [] (to which the reader is referred for more details).

Hydrothermal model
The circulation of hot fluids in the hydrothermal system is modelled by TOUGH [], a well-known software to simulate multiphase multicomponent fluid flow and heat transfer in porous media. For a generic multiphase fluid with k components a set of k+ equations (k mass equations and one energy equation) is solved. The k +  equations can be resumed as follows []:   T e m p e r a t u r e U n k n o w n Dependence on x and/or t is also highlighted in the first column. All quantities that do not depend on time are a priori assigned. The subscript β refers to the liquid (β = l) or gas (β = g) phase.
where the subscript β = l or g refers to the liquid or gas phase, respectively, and F β = In the simulations considered in this paper, we use k =  or k =  to simulate pure water (EOS module of TOUGH) or water and carbon dioxide (EOS module of TOUGH) injection, respectively.
We call H ⊂ R  the spatial domain in which we solve Eq. () (the domain of the hydrothermal model). Due to the axis-symmetric configuration of the problem, the domain can be expressed in terms of the radial distance r = x  + y  and the vertical variable z, as H = { ≤ r ≤ R, z min ≤ z < f TOP (r)}, where f TOP (r) represents the topography function (see Figure ). Since TOUGH does not account for infinite domains, H differs from M and is truncated at a radial distance of r = R and a depth of z = z min . On the portion of the ground surface within a radial distance r ≤ R, H , atmospheric (Dirichlet) boundary conditions for pressure and temperature are prescribed.
In multi-phase system, each phase may be at a different pressure P β due to interfacial curvature and capillary forces. The difference between the gas and liquid pressures is referred as capillary pressure P c . Capillary pressure and relative permeabilities k β , β = l, g, are usually posed as a function of the liquid saturation S l . In this paper we use the capillary pressure and relative permeabilities set according to the following Brooks-Corey functions []: where S e is the effective saturation, S lr = . the residual liquid saturation and S gr = . the residual gas saturation.

Simulation of unrest at Campi Flegrei
The H for the TOUGH simulations extends up to a radial distance of R =  km and a depth of z min = . km computed from the intersection between the ground surface and the axis of symmetry. The domain H was discretized in the radial direction by a set of logarithmically spaced nodes with a horizontal resolution starting from . m along the axis of symmetry and decreasing to  m at the external boundary. Vertically, the domain was divided into  equally spaced layers, which corresponds to a resolution of  m/layer. This discretization leads to , cells.
In the thermo-poroelastic and gravity models the domain extends toward infinity (Figure ), warranting zero displacement and potential. The computational grid is vertexcentered and has a maximum resolution of r = z = . m at r = z = .
Firstly, TOUGH is run to evaluate temperature, pressure and density variations with respect to their initial distributions, which, then, are fed into the thermo-poroelastic solver to compute the deformation and gravity changes. Due to the different grids adopted, these quantities are interpolated from the TOUGH grid to the thermo-poroelastic solver nodes.
In all the TOUGH simulations, atmospheric boundary conditions P = . MPa and T =  • C are prescribed on the surface H = {z = f TOP (r)}, while adiabatic and impervious     Table ). The elastic medium properties are defined on the basis of tomographic studies [], which depict the heterogeneity of the shallow structure of Campi Flegrei (Table ). Following [], the shallow area of the medium is divided into three horizontal layers having a thickness of  km. The inner caldera is modelled as a  km wide cylinder, coaxial with the axis of symmetry, with an internal variation of the elastic parameters. A further layer extends from  km depth to the bottom of the (infinite) computational domain ( Figure ). The thermal expansion coefficient α is  - K - and the Biot-Willis coefficient α BW is ( -K/K s ), where K is the isothermal drained bulk modulus ( GPa) and K s is the bulk modulus of the solid constituent ( GPa) [].
Although hydrothermal fluids are not pure substances, but generally mixtures of several mass components or chemical species, the dominant fluid component is usually water, and it is often reasonable to ignore other components. However, in volcanic regions the aqueous phase generally contains also dissolved incondensable gases, such as CO  . To investigate the effect of a CO  component in water, two different scenarios are investigated, simulating: () a pure water injection using the TOUGH/EOS module, which models a water system in its liquid, vapor, and two-phase states; () a mixture of water and carbon dioxide injection using the TOUGH/EOS module, which accounts for the non-ideal behavior of gaseous CO  and dissolution of CO  in the aqueous phase.
We observe that the depth of H is chosen in such a way the focus of the simulations is on the shallower hydrothermal activity, in order to maintain temperature and pore pressure within the range considered by TOUGH, which does not take into account super-critical fluids. Therefore, in all the simulations of this paper we investigate only sub-critical fluids.
Observe that we do not aim to simulate the subsidence period following the main unrest phase, and therefore the simulation is confined to the first years of unrest.

Injection rates
In the following sections two scenarios will be studies: pure water injection and water and carbon dioxide injection. The injection rates adopted in the scenario of water and carbon dioxide injection are obtained from previous works according to geochemical data collected at the CFc and providing a good matching to observed data [, , ]. The injection rates adopted for the scenario of pure water injection are chosen in such a way the total mass of fluid injected is the same for both scenarios.

Pure water injection
In this first scenario the fluid is composed of pure water, whose properties and phase transitions are calculated based on the thermodynamic conditions, according to the steam table equation []. Initial conditions are obtained reaching a steady-state solution by simulating a  thousand year long phase with a deep injection of hot water with a constant flux rate of . kg/s at temperature of °C. The fluid is injected in a  m wide inlet located at the bottom of the domain around the symmetry axis. This extent has been estimated from the area affected by fumaroles activity [].
The unrest is simulated increasing the flux rate to . kg/s for  years. The distributions of saturation and pressure, temperature and density changes at the end of the unrest with respect to the steady-state initial conditions are displayed in Figure . Temperature changes are restricted near the inlet, whereas pressure changes of a few MPa are distributed in a larger area. The saturation shows that the fluid is gas-saturated at the base and a two-phase fluid reaches the ground surface. The saturation distribution controls the fluid density, which shows negative variations at the edge of the inlet and positive variations at a mid-depth of the domain. Using these solutions the evolution of ground deformation and gravity changes are then evaluated.
The deformation field ( Figure ) computed solving Eq. () resembles the patterns obtained in [-], although there are some discrepancies in the amplitudes, likely arising from the medium heterogeneity rock properties and the single component injection. It is similar to deformation field arising from spherical point pressure sources. During the unrest the deformation monotonously increases reaching an amplitude of . cm and . cm in the horizontal and vertical component, respectively.
A comparison with a homogeneous model, in which the rigidity is  GPa and the Poisson's ratio is . is performed. An enhancement in the deformation field is observed since the rigidity is lower with respect to the heterogeneous model ( Figure ).

Water and carbone dioxide injection
In this second scenario the fluid is a mixture of water and carbon dioxide. The initial conditions are achived by injecting for  thousand years a fluid at temperature of °C composed of a water component with a flux rate of . kg/s and a CO  component with a flux rate of . kg/s. Then, the unrest is simulated increasing the flux rates to . kg/s for the water and to . kg/s for the CO  . The distributions of saturation and pressure, temperature, and density variations after  years from the unrest are displayed in Figure .

Conclusion
A finite-difference ghost-point method for solving the elasto-static and Poisson equations on an arbitrary unbounded domain has been presented. The proposed strategy, which adopts the coordinate transformation method, has been applied to thermo-poroelastic models to evaluate deformation and gravity changes in volcanic areas. The method has been coupled with the TOUGH code to investigate the role of pressure, temperature and density changes on the geophysical observables. The results from two scenarios have been compared in order to explore the effect of the presence of carbon dioxide in the fluid mixture. Although the total flux rate is similar in both scenarios, the presence of CO  strongly alter the solutions. Particularly, the injection of the CO  engenders an enhancement in the deformation field and perturbs the temporal evolution of the gravity changes. This last feature could be a useful signature to constrain the relative ratio between water and carbon dioxide content. It turns out that the ratio between gravity and vertical deformation de-  creases for higher ratios of CO  /H  O. The comparisons between gravity and deformation may help to discriminate between injection with or without carbon dioxide content and to provide inferences on the nature of the source. Since unrest periods are accompanied by increases of the CO  /H  O ratio, a decrease in the ratio between gravity and vertical deformation could help in detecting the onset of unrest phases. Moreover, our findings seem able to justify what observed in some volcanic regions, where major gravity changes appear without any significant deformation. In many volcanoes worldwide, there are evidences that the ratio between gravity and height changes is far beyond what could be predicted by simple models, in which volume and pressure changes in a magma chamber are considered. When significant gravity changes occur without any significant deformation, or vice versa, it is often difficult, if not impossible, to jointly explain the observations using the popular Mogi model. Our results may provide an alternative explanation to the observations and help in resolving the controversy between geodetic and gravity observations as a volcano moves from rest to unrest state.
However, there are some limitations of the model presented in this paper that must be considered. The shallow hydrothermal system is only . km deep, while a deeper model is more realistic to represent volcanic areas [, ]. The shallow injection depth considered in this paper is constrained by the range of applicability of TOUGH, which does not account for supercritical fluids. Some recent models investigated the influence of super-critical fluids to the results, allowing to consider deeper domains []. The D axi-symmetric representation may not be able to describe the complex D network that characterises caldera systems. However, to a first approximation calderas present a radial structure and then can be approximated by the axis symmetric assumption. The one-way coupling between the hydrothermal model and the thermo-poroelastic model represents a reasonable simplification when a short period of unrest is simulated, but it is not adequate for simulations of longer processes, where stress and strain alterations may induce a significant variation in the relevant hydrological parameters (permeability, porosity), modifying the long term processes of fluid flows in the porous medium and then the associated geophysical signals [, ].
The effects of these limitations may be reduced by considering a more realistic fully coupled D model that accounts for super-critical fluids. The implementation of a more effective simulator is part of our future investigations, with the aim to extend the simulation to more realistic cases of multi-component fluids.