Influence of permeability on the hydrothermal system at Vulcano Island (Italy): inferences from numerical simulations

Volcano-hydrothermal systems are governed by complex interactions between fluid transport, and geochemical and mechanical processes. Evidence of this close interplay has been testified by distinct spatial and temporal correlations in geochemical and geophysical observations at Vulcano Island (Italy). To understand the interaction between fluid circulation and the geochemical and geophysical manifestations, we perform a parametric study to explore different scenarios by implementing a hydro-geophysical model based on the equations for heat and mass transfer in a porous medium and thermo-poroelastic theory. Numerical simulations allow us to define the controlling role of permeability distribution on the different modeled parameters as well as on the geophysical observables. Changes in the permeability within the highly fractured crater area could be responsible for the fluctuations in gas emission and temperature recorded during the crisis periods, which are accompanied by shallow volcano-seismicity in the absence of significant deformation and gravity variations. Despite the general medium permeability of the volcanic edifice, the presence of more highly permeable pathways, which allow the gas to rapidly escape, as testified by the presence of a well-developed fumarolic field, prevents the pressure buildup at shallow depths.


Introduction
In recent decades, a number of efforts have been made to improve the comprehension of the interplay between fluid flow and mechanical rock response (Chiodini et al. 2016;Ingebritsen et al. 2010;Rinaldi et al. 2010;Tanaka et al. 2018;Troiano et al. 2011). Thermo-poro-elastic approaches, coupling fluid flow and rock strain, have been explored in several volcano geothermal areas, which have demonstrated the important role of elastic and transport properties in leading to changes in geochemical and geophysical observables at the surface related to hydrothermal activity (Chiodini et al. 2016;Peiffer et al. 2018;Rinaldi et al. 2010;Tanaka et al. 2018). Vulcano Island, whose hydrothermal activity has been characterized by different unrest phases since 1890, when the last eruption occurred (Paonita et al. 2013;Selva et al. 2020), may well be considered a representative case. On the island, during the last 30 years, the main unrest episodes, lasting no more than a few months, develop with the increase in magmatic components and in temperature of the fluids emitted through the crater fumaroles. In particular, fumarole temperatures, 350 °C on average, show pulses reaching values up to 700 °C (Nuccio et al. 1999;Diliberto 2011Diliberto , 2013. The emission of CO 2 , the most abundant anhydrous component in the fumarolic gases of the La Fossa cone (Inguaggiato et al. 2012), increases up to 1000 t/d, while in quiescent periods, it is 482 t/d on average. These increases in the fluid emission and in the temperature have been accompanied by the occurrence of shallow micro-seismicity at a depth of about 500-1500 m below La Fossa crater ). On

Open Access
*Correspondence: santina.stissi@ingv.it Stissi et al. Earth, Planets and Space (2021) 73:179 the contrary, no significant ground deformation or deep volcano-tectonic seismicity has been observed (Alparone et al. 2019). Concurrently, gravity changes (Di Maio and Berrino 2016) did not show any temporal correlation with these cyclic pulses in the fumarolic gas contents and temperature. These episodic and abrupt geochemical changes seem to involve the shallow hydrothermal system, which can be provoked by the input of deep magmatic fluids with no contribution of direct magma movement upward. In this framework, massive escape of fluid from the deep magma reservoir or the activation of new degassing levels (Paonita et al. 2013) represents the most likely unrest scenario (Selva et al. 2020), but we cannot rule out the possibility of an enhanced (sea) water inflow at depth. Preliminary investigations on the unrest phases of the hydrothermal system have already been attempted by  and , who evaluated the imprint of the hydrothermal activity on the geophysical signals. In those previous studies, with the aim of simulating multiphase fluid flow in porous media undergoing deformation, two well-known numerical codes were coupled: (i) TOUGH2-EOS2 for nonisothermal CO 2 -H 2 O multiphase flow in porous media (Pruess et al. 1999) and (ii) COMSOL Multiphysics for computing thermo-poro-elastic deformations and gravity changes (COMSOL 2012). Numerical thermo-poroelastic models made it possible to investigate several scenarios and assess their suitability by comparing the simulation results with the available observations. However, limitations on the investigations of fluid temperature were posed by the use of the TOUGH2-EOS2 module, which can only work at temperatures in the range between 0 and 350 °C, below the water criticalpoint temperature. This low-temperature range implies that the behavior of super-critical fluids, present in the hydrothermal areas, cannot be properly simulated, hampering a realistic evaluation of the fluid dynamics and rock response. Indeed, the temperature values measured at the crater fumaroles (Diliberto 2013) and the geochemical model (Federico et al. 2010) indicate an average temperature of 400 °C in the shallow hydrothermal system, and hence, much higher values for the deep input fluid could be expected. The fluid temperature not only affects the thermo-elastic response of the rock but also the associated fluid density variations and, hence, it influences both the deformation and gravity field at the surface.
To overcome this limitation and investigate more realistic models, here, we explore the use of MUFITS with the BINMIXT module (Afanasyev 2012;Afanasyev 2013a, b). MUFITS determines phase equilibria of the mixture in sub-and super-critical regions of the hydrothermal system and simulates both the sub-and super-critical condition of CO 2 , modeling a non-isothermal three-phase flow of the mixture, in which two different phases of CO 2 are present: liquid like and gas like phases. By means of these numerical simulations, a set of elastic and transport properties are examined to reconcile geochemical, deformation, and gravity observations.

Geological setting
Vulcano Island is a composite volcanic edifice located approximately 20 km north of the Sicily coast (Fig. 1a). The island, belonging to the Aeolian Archipelago, rises 500 m a.s.l., but represents the emerged part of an NW-SE elongated submarine volcanic complex over 1000 m high (Romagnoli et al. 2013). The subaerial eruptive activity was characterized by long periods of quiescence interrupted by explosive eruptions that emit lava flows and small volumes of pyroclastic rocks (Nicotra et al. 2018). It has progressively migrated from SE to NW, strongly controlled by regional tectonic and influenced by the Tindari-Letojanni fault system (Ventura et al. 2013;Ruch et al. 2016), a right-lateral strike-slip fault that trends northward from the Sicily coast to the center of the Aeolian Archipelago (Selva et al. 2020). Since the last magmatic eruption in 1890, Vulcano has undergone frequent fumarolic activity, located in two main exhalative areas, the Baia di Levante area, where fumarole temperatures are lower than 100 °C, and La Fossa crater which is characterized by the presence of a high-temperature fumarolic field where gas emissions currently reach temperatures of around 400 °C (Madonia et al. 2015). The distribution and appearance of fumaroles at Vulcano Island are closely related to the tectonic, geologic, and volcanic history of the volcanic complex (De Astis et al. 1997;Diliberto 2017). In particular, at the La Fossa cone, multi-scale observations revealed that hydrothermal vents are preferentially located on surface fractures and cracks controlled by the main structural features (Schöpa et al. 2011). These fumarolic emissions bear witness to the existence of a very active hydrothermal system below the La Fossa cone. There is extensive geophysical literature (Chiarabba et al. 2004;Blanco-Montenegro et al. 2007;Alparone et al. 2010;Napoli and Currenti 2016;Ruch et al. 2016) on its presence at a depth between 0.5 and 1.5 km below sea level (b.s.l.). Moreover, two geothermal drillings (Gioncada and Sbrana 1991) and several geochemical studies enabled imaging the inner structure of La Fossa cone, identifying the characteristics of the hydrothermal system and modeling its hydrothermal circulation (Revil et al. 2008(Revil et al. , 2010Barde-Cabusson et al. 2009). Several anomalous degassing periods, with significant temperature variations, changes in chemical composition and increases of exhaling areas at the summit area of the La Fossa cone, have been observed in the last decades, such as during the crises of 1996, 1998, and 2004-2005(Chiodini et al. 1996Granieri 2006;Diliberto 2017). These periods, that may also be accompanied by local shallow seismicity below La Fossa crater ), indicate periodic fluctuations in the level of volcanic activity, due to interaction of the magmatic and the hydrothermal systems.

Numerical modeling
With the aim of comprehending the interplay between fluid circulation and the geochemical and geophysical manifestations, we implement a hydro-geophysical model, based on hydrothermal and thermo-poroelasticity theory. The reservoir simulator MUFITS (Multiphase Filtration Transport Simulator) is used for hydrodynamic simulations of non-isothermal multiphase flows of the CO 2 -H 2 O two-component fluid in a porous medium. Assuming that the fluid is in local thermal equilibrium with the rock, MUFITS solves the mass balance equation (Eq. 1) (one for each component j = 1, 2 of the mixture), the energy balance equation (Eq. 2), and Darcy's law (Eq. 3), discretized by finite volume method and the upwind scheme for fluxes, to compute the distribution of gas saturation, composition of the fluid phases, temperature, pore pressure, and density variations In Eqs. (1-3), the given parameters are the porosity φ , the rock grain density ρ r , the thermal conductivity of saturated rock ∼ , the absolute permeability of the rock K, the gravity acceleration g, and the source terms (the influx of the jth component Q j or the energy Q E ). The primary variables for which the equations are solved are the pressure, the total fluid enthalpy, and its total composition. The temperature T, the internal energy of the rock e r , and for each phase i, the density ρ i , the saturation S i , the jth component mass fraction c ij , the mass enthalpy h i , the internal energy e i , the relative permeability k i , and the viscosity η i are nonlinear functions of the primary variables, which are defined by an iterative numerical procedure for predicting the fluid phase equilibria (Afanasyev 2012, b). The mass flow rate of phase i, F i = ρ i w i = −K k i η i ρ i gradP − ρ i g , follows Darcy's law (3) and w i is the Darcy velocity. Substituting (Eq. 3) into (Eqs. 1 and 2) allows to exclude w i from the system of governing equations and consider (Eqs. 1 and 2) as a closed system on the primary variables, i.e., the pressure, the total fluid enthalpy, and its total composition. Using these primary variables, MUFITS avoids singularities in the vicinity of the critical thermodynamic conditions, which provides a greater numerical stability, and avoids the reduction of the simulation time step under near-critical conditions ensuring a fast convergence (Afanasyev 2013a, b).
The circulation of a three-phase (i = 1, 2, 3) hot fluid composed of a mixture of H 2 O and CO 2 can be simulated using the BINMIXT module of MUFITS that employs the described governing equations. It works in the temperature and pressure range of 0-1000 °C and 0-350 MPa, respectively. It enables simulating properties of a binary CO 2 -H 2 O mixture under both subcritical and super-critical conditions (with CO 2 -rich and H 2 O-rich phases). The CO 2 -H 2 O fluid can be in 1-phase, 2-phase, and 3-phase equilibria at temperatures above ~ 0 °C (Afanasyev 2013b). The three phases are the liquid phase (mainly consists of H 2 O), the gas phase (H 2 O vapor and gaseous CO 2 ), and the liquefied CO 2 phase (a liquid phase mainly consisting of CO 2 ). This third phase and, thus, the 3-phase equilibria are possible only at low temperatures (generally below critical point of CO 2 , T < 31 °C) and high pressures. Such a condition is not encountered in our simulations, in which only 1-phase and 2-phase equilibria occur.
The thermo-poroelasticity model  is implemented using the finite-element COM-SOL Multiphysics software (COMSOL 2012). We solve the elastostatic equation and the Poisson equation for the gravitational potential, to evaluate ground deformation and gravity changes induced by pore pressure, temperature, and density changes obtained by MUFITS simulations.
The displacement is calculated by solving the stress equilibrium Eqs. (4-6) coupled with the thermo-poroelastic extension of Hooke's law (5), obtained by adding the pore pressure changes P through the Biot-Willis coefficient β = 1 − K d K s ( K d and K s drained and solid bulk moduli, respectively) and the temperature changes T through the volumetric thermal expansion coefficient α T (Fung 1965;Jaeger et al. 2007) In the Eqs. (4-6), σ and ε are the stress and strain tensors, respectively, u , the unknown of the problem, is the deformation vector, and I the identity matrix. The Lame's elastic medium parameters and µ are related to the Poisson's ratio and Young's modulus. The problem is closed by imposing zero displacements at infinity and stress-free boundary condition σ · n s = 0 on the ground surface, n s being the normal vector at the ground surface. Pore pressure and temperature changes also alter the fluid density distribution, which in turn affects the gravity field. The gravity changes g induced by the density variations are evaluated by solving the gravitational potential ϕ g in the following Poisson equation (Coco et al. 2014): in which G is the universal gravitational constant. The density variation �ρ f of the fluid, which occupies the porous volume fraction φ (rock porosity), is computed as follows: �ρ f = φ ρ f − ρ f 0 , where ρ f = i ρ i S i is the density at a generic time, obtained by adding over the phases i their densities multiplied by the phase saturation S i , and ρ f 0 is the fluid density at a reference time. The problem is closed by imposing that the gravitational potential is zero at infinity.
Specific subroutines have been implemented to automate the data transfer between MUFITS and COMSOL and to routinely evaluate the spatio-temporal evolution of ground deformation and gravity changes ).

Modeling set-up
The model has been implemented in an axis-symmetric formulation, suitable for describing the radial structure of the main volcanic edifice, whose details can be found in . Along the NE-SW profile (AA') selected to represent the volcanic edifice ( Fig. 1), the computational domain, extending from the topographic surface down to 1500 m b.s.l., was discretized with a cell-centered grid of nodes, with a constant vertical resolution of 20 m, while the horizontal resolution decreases from 20 m, near the axis of symmetry, to 216 m at the external boundary at 10 km away.
To define the heterogeneous structure and hydrological properties of the model, we started from a simplified layered structure (Fig. 1). Along the AA' section, we identified different regions (Fig. 1b) on the basis of the available geological and geophysical data (Gioncada and Sbrana 1991;Blanco-Montenegro et al. 2007;Napoli and Currenti 2016). The first one (FC) extends from the symmetry axis for about 100 m and represents the main conduit, namely the preferential central pathway for hydrothermal fluids within La Fossa crater. FC is bounded by a transition zone, consisting of strongly altered Primordial Vulcano products (AP) extending from the bottom to a depth of about 500 m b.s.l. At that depth, AP is replaced by La Fossa caldera filling products (CP) topped by tephritic lavas (TL) and pyroclastics (PY). At about 500 m from the axis of symmetry lavas of the Primordial volcano extend for the whole domain (PL).
The rock hydrological properties (absolute permeability, density, porosity, thermal conductivity, and specific heat capacity) were attributed to the different regions (Tables 1, 2). Owing to the lack of detailed values, the lithological units are grouped into different classes whose permeability values were chosen based on the values of similar rock units and related modeling studies Okubo and Kanda 2010;Rinaldi et al. 2011;Troiano et al. 2011;Coco et al. 2016). In particular, a permeability of 5 × 10 −15 m 2 was assigned to TL and PL, and the intermediate value of 10 −14 m 2 was chosen for PY and AP, while the highest value of 5 × 10 -14 m 2 was defined for FC and CP ( Table 1). The porosity values, ranging between 0.1 and 0.4, are instead obtained from laboratory tests (Tommasi et al. 2016) carried out on the rock samples collected from the BL1 borehole, which is 100 m deep and located at the base of the northern flank of La Fossa cone, very near to VP1 well ( Fig. 1; Gioncada and Sbrana 1991).
The relative permeability follows Corey's law (Brooks and Corey 1964): where k l and k g represent the relative permeabilities of the liquid and gas phases, respectively, and S = (S l − S lr )/(1 − S lr − S gr ) , with S l the liquid saturation, S lr = 0.33 the irreducible liquid saturation, and S gr = 0.05 the irreducible gas saturation.
The fluid viscosity changes with pressure, temperature, and composition using a corresponding state model (Afanasyev 2013a).
The elastic parameters were defined in the thermoporoelastic simulator to evaluate the rock mechanical response to the circulation of the hydrothermal fluid. Since at Vulcano Island, there is no detailed information on the volumetric thermal expansion coefficient α T and the solid bulk modulus K s , we assume for these parameters the average values of literature data acquired in a similar volcanic environment, which are summarized in Table 2 Troiano et al. 2011;Coco et al. 2016). The values of the drained bulk modulus K d and of the Young's modulus E , that express the rock elastic properties, are derived from seismic tomography investigations (Chiarabba et al. 2004) and gravity studies (Budetta et al. 1983). In particular, the Young's modulus E , depending on the average density ρ of the heterogeneous medium and on P-wave seismic velocity V p , is defined according to the following formula (Kearey and Brooks 1991): where ρ and V p are assumed to increase linearly with depth and the Poisson's ratio ν is equal to 0.3. In this way,  The initial and boundary conditions are defined as follows. Atmospheric pressure and ambient temperature are set along the upper ground surface boundary. In the submerged boundary, the hydrostatic pressure loading is assigned according to the bathymetry. Hydrostatic pressure and a geothermal gradient of 130 °C/ km (AGIP 1987) are assigned in the whole domain. The molar composition in the domain and at the upper boundary is obtained from Dalton's law, assuming an atmospheric value of CO 2 partial pressure.
The initial conditions are computed as steady-state solutions by simulating for about 4500 years the injection of a mixture of water and carbon dioxide at a temperature of 600 °C through a 100 m-wide inlet placed at the bottom of the domain near the axis of symmetry. At the inlet, Neumann boundary conditions are imposed, fixing the flow rate of the mixture. Also, the enthalpy and composition of the injected fluid are specified, closing the mathematical problem in Eqs. (1-3). Therefore, the temperature is not fixed during all the simulation time and could also change near the boundary due to the switching from quiet to unrest period when the flow rate is enhanced.
Starting from these initial conditions, the unrest phase, characterized by a greater contribution of the flow rate, was simulated for a period of 1 year. The injected CO 2 reflects the values from monitoring of CO 2 output flow in a 500 × 500 m crater target area (CTA, Fig. 1) of La Fossa cone in the period 1995-2006 (Granieri et al. 2006). In this time span, the minimum value of 68 t/d was recorded in correspondence of a quiet period (March 2002), while a maximum value of 330 t/d was observed in December 2004, when the increase in the quantity of CO 2 (Granieri et al. 2006) and in the vent temperature (Diliberto 2011) are recorded, marking the culmination of the crisis that began in November 2004. Considering that these outflow rates of CO 2 were measured for CTA (Granieri et al. 2006), we assume for the entire crater that is about 2.2 times larger than CTA, an outflow rate of CO 2 of 150 t/d during a quiet period. Assuming a balance between injection and outflow, the same value of the injected CO 2 was considered for the steady state and a value of 726 t/d was assumed for the unrest. Based on the measurements gathered by Chiodini et al. (1996) and Inguaggiato et al. (2012), a flow rate of H 2 O of 433 t/d was defined for simulating the steady state and of 1485 t/d for the unrest phase, corresponding to measurements performed during quiet and unrest periods.

Parametric study
To better understand the observed variations at the surface induced by the thermal and mechanical response of the porous medium subjected to changes in the hydrothermal system, we perform a parametric study on some crucial parameters, which regulate the spatial and temporal evolution of some monitored quantities at Vulcano Island. A full parametric study is hindered by many uncertain parameters and the expensive computational burden, and therefore, our investigations are limited to a range of more sensitive parameters. First, particular attention is focused on the rock permeability that plays a fundamental role in guiding the ascent of fluids from the hydrothermal source to the surface, and in regulating the vent temperature of the fumaroles. The permeability reduction in the zone through which fluid reaches the surface, indeed, may produce an increase in the pressure, a buildup of heat below the obstruction zone and a decrease in vent temperature (Harris et al. 2012;Tanaka et al. 2017). On the other hand, an increase in permeability can cause a faster rise of fluid and larger ground deformation ).
Consequently, a period of heating may be due to an increase in permeability in the most superficial part of the system due to rock fractures and not necessarily to a greater contribution of heat of the source and of magma components of the fluid (Cannata et al. 2012). The increases in the occurrence of seismic hybrid and high-frequency events observed during unrest phases, indeed, testify to the occurrence of fracturing processes , induced by pore pressure changes and thermal stresses, of the rocks surrounding the hydrothermal system, as formulated by the thermo-poro-elastic models (Eqs. 4-6). Thermal stress and pore pressure changes could be greatly enhanced by the inflow of cold water and boiling processes due to a successive heating up, that could induce relevant microseismic activity (Alparone 2010,2019) and fracturing. A further contribution for rock fracturing could arise from the rock weakening due to a geochemical reaction between the rock matrix and corrosive acidic fluids, due to the presence of dissolved gas such as sulfur dioxide, chloridric, and fluoridric acid that could alter the rock properties (Farquharson et al. 2019).
To evaluate the effects produced by changes in the permeability of the domain and of the central conduit FC, within La Fossa crater, we performed a set of simulations under different configurations and parameters (Additional file 1: Table S1).

Initial permeability configuration
Two different permeability configurations have been considered for the domain. The former configuration, HP (higher permeability), is characterized by values ranging between 5 × 10 -14 and 5 × 10 -15 m 2 according to , while, in the latter configuration LP (lower permeability), permeability values are one order of magnitude smaller than those in HP. In both permeability configurations, a 100 m-wide central conduit, characterized by higher permeability values (5 × 10 -14 m 2 for HP and 5 × 10 -15 m 2 for LP), is assumed near the axis of symmetry to simulate a preferential pathway for the flow of water and carbon dioxide toward the surface within La Fossa crater.
Temperature (T), pressure (P), and gas saturation ( S g ) distributions obtained in the steady state (Additional file 1: Fig. S1) and for an unrest phase (Fig. 2) are evaluated in all the model domain for both HP and LP configurations. At steady state, higher values of P, T, and S g are obtained in the LP configuration where they affect a wider area compared to the HP configuration. By contrast, at the end of the unrest, positive temperature changes are mostly limited to the inlet zone and are higher in the HP configuration, where variations of about 200 °C, with respect to the steady-state solution (Additional file 1: Fig. S1), are attained (Fig. 2c, d). Pressure changes are higher in the LP configuration, where the maximum value of about 9 MPa is reached close to the inlet area at the depth of 1500 m b.s.l. These pressure changes decrease toward the surface and become lower than 1 MPa at a depth less than 1000 m b.s.l. (Fig. 2f ). On the other hand, in the more permeable HP configuration, pressure variations are less than 1 MPa, but are noted, in a wider area, up to 1 km from the axis of symmetry and affect the superficial area up to the sea level (Fig. 2e). Gas saturation distributions show different patterns in the two configurations. Indeed, in the HP configuration (Fig. 2g), despite being more permeable, it is possible to observe a gas saturated zone only in the inlet zone, and in the most superficial area, but not within the central conduit. On the contrary, in the LP configuration (Fig. 2h), gas is present throughout the conduit, from the inlet area to the ground surface, laterally extending for about 500 m from the axis of symmetry. It is worth noting that at the sea level, a thin gas saturated area spreads laterally in the shallow subsurface in correspondence of the lithological discontinuity between tephritic lavas TL and pyroclastics PY, where a contrast in permeability is present.
Although a greater pressurization of the system occurs in the less permeable LP configuration, ground deformations are greater in the HP configuration (Fig. 2i-l). This is due to the different depths of the pressure buildup that is too deep in the LP configuration to generate significant ground deformation at the surface. The simulations reveal that the pore-pressure-induced strain contributes more than the thermo-elastic effect (Fournier and Chardot 2012). The horizontal deformations of a few millimeters are observed in both HP and LP configurations, but they reach their maximum amplitude at different distances from the axis of symmetry, about 1200 m in HP and about 2000 m in LP. The vertical deformations reach the maximum amplitudes of about 1 and 2 cm, in LP and HP configurations, respectively. In both configurations, the maximum is located above the source on the axis of symmetry and decays radially. In the HP configuration, the vertical deformation decays more rapidly away from the axis of symmetry. No significant variations in total density are detected, except for some areas where negative density changes (Additional file 1: Fig. S2) are noted for both permeability configurations as lower density gas-rich fluids replace the colder and higher density fluids. The gravity changes are greater in the HP configuration, where they reach about 13 µGal, while in the LP configuration, they do not exceed 2 µGal.
The gravity change is prominent right above the hydrothermal source, and, in the same way as the horizontal deformation, it vanishes at about 2000 m and 2500 m in HP and LP, respectively.
An increase of two orders of magnitude of the permeability in the central conduit with respect to the initial configuration during the unrest phase is successively simulated. The effects on the modeled parameter changes, T and P , with respect to the initial steady state in Additional file 1: Fig. S1, the distribution of S g and on the geophysical observables are shown in Fig. 3. A depressurization of the system is seen in both configurations regardless of the different permeability values. In particular, negative pressure variations are more significant in the LP configuration, where we can observe a value of about − 8.5 MPa in the inlet zone that gradually increases along the conduit reaching the value of about − 0.7 MPa at a depth of about 500 m b.s.l. (Fig. 3f ). In the HP configuration, the pressure change reached a minimum of about -0.2 MPa around the inlet and on the ground surface, while along the central conduit from the sea level up to a depth of about 1330 m b.s.l., positive variations up to about 1 MPa are prominent up to 1 km from the axis of symmetry (Fig. 3e). Temperature variations affected the entire conduit in the HP configuration. In particular, a cooling develops from the inlet area to the surface, where it reaches the values of − 180 °C. This cooling affects the distribution of both gas saturation (Fig. 3g) and the total density variations (Additional file 1: Fig. S3e). Indeed, a gas saturated zone is present only at the ground surface where temperatures reach the lowest values (Fig. 3g). Moreover, the maximum density variation (about 1000 kg/m 3 , Additional file 1: Fig. S3e) is reached in the shallowest part of the central conduit where a very dense fluid, rich in water, has infiltrated into the gas zone. Also, in the LP configuration, a cooling zone is present, but only in the proximity of the inlet where a temperature change of − 20 °C is reached. No significant temperature variations are noted along the conduit except in a small superficial area where a heating of about 80 °C is present. The density variations have negative values in a thin superficial area, indicating the rise of a fluid with a lower density compared to the steadystate condition (Additional file 1: Fig. S3f ), while no significant variations are obtained in the whole domain. The depressurization and the cooling of the system produce significant effects on the ground deformation and on the gravity changes ( Fig. 3i-n). Thermo-elastic effects greatly contribute to ground deformation in the HP configuration (Additional file 1: Fig. S4i, k, m, o) and show less and less variability after 0.6 year. On the contrary, the contribution of pore pressure to ground deformation continues to vary over time (Additional file 1: Fig. S4m, o). Negative pore pressure changes contribute to most of the deflation in the LP configuration. The heating in the surface zone of LP contributes to positive peaks of a few millimeters in horizontal deformation in correspondence of the edge of the two-phase plume (Additional file 1: Fig. S4j, l, n, p). Depressurization is accompanied by negative gravity changes in the LP configuration (Fig. 3n), while in the HP configuration, positive gravity changes with a maximum amplitude of about 500 μGal are prominent right above the source (Fig. 3m). This is attributable to the positive density variations in the shallow part of the main conduit (Additional file 1: Fig. S3e), reflecting the inflow of liquid water.

Size of the inlet and the central conduit
The sizes of the inlet and the central conduit can influence the distribution of temperature, pressure, and gas saturation as well as the geophysical observables. To evaluate their effects, the steady state and the unrest phases are simulated for the HP and LP configurations, by varying them in the range 40-200 m and keeping the same flow rate. The results are compared with those already obtained using a 100 m-wide inlet and central conduit.
The simulation of the steady state reveals that narrower sizes of both the inlet and the central conduit lead to higher surface temperatures (Additional file 1: Figs. S5 and S6a-c) and in the LP configuration produce an enlargement of the area affected by heating. No significant variations of the pressure distribution are obtained in the HP configuration, while a pressure buildup in the inlet zone appears in the LP configuration. The gas saturation is insensitive to the size changes in the LP configuration, while in the HP configuration with the narrower size, it spreads laterally in the shallow subsurface and through the central conduit. By narrowing the conduit, an increase of the average mass flux rate is attained. In particular, when the conduit width reduces from 200 to 40 m, the average mass flux rate within the conduit area increases non-linearly from 0.011 to 0.032 t/d/m 2 and from 0.0037 to 0.0076 t/d/m 2 , for the HP and LP configurations, respectively. In the unrest phase (Figs. 4,5), for both the wider (200 m) and the narrower (40 m) conduit, the temperature changes are restricted to the deepest part of the domain, close to the boundary, so they are not very relevant for the processes at the ground surface. Significant variations in the distributions of temperature changes are achieved at the edge of the plume in the HP configuration, where the temperature change reaches 100 °C at about 90 m and 190 m from the axis of symmetry if widths of 40 m and 200 m are used, respectively (Fig. 4d, f). Also in the LP configuration, the temperature change reaches 100 °C at about 190 m from the axis of symmetry if a width of 200 m is used (Fig. 5f). The temperature variations could be ascribed to the  changes in the fluid flow direction. In both HP and LP configurations (Figs. 4,, the pore pressure significantly increases if the inlet becomes narrower, and decreases for a wider inlet affecting a wider area. Different gas saturation, distributions are noted for different sizes of the inlet and the conduit for the HP configuration, where the narrowest size produces a wider saturated area at the ground surface. On the other hand, there is no significant variation in the gas saturation distribution in the LP configuration, except that the horizontal extension of the gas saturated area near the sea level is greater for a conduit width of 40 m. Ground deformation and gravity changes are gradually enhanced as the size of the inlet increases for both HP and LP configurations. The exception is the gravity changes after 1 year from the unrest phase for the HP configuration, when the inlet and the main conduit sizes are 200 m (Fig. 4u). In this case, from 0.8 to 1 year, the gravity changes decrease from about 13.4 μGal to 12.4 μGal. The different parametric simulations highlight the controlling role of permeability and the size of the central conduit on the different modeled parameters' distribution as well as on the geophysical observables. In particular, in the unrest phases, we observe that significant variations of temperature at the surface are reached only when the permeability of the conduit is enhanced (Fig. 3). Increasing the flow rate of the mixture and the CO 2 content is not enough to produce significant temperature variations at the ground surface. The flow rate increase induces temperature changes that are concentrated in the inlet zone and cannot propagate upwards toward the surface (Fig. 2).

Narrow inlet and conduit with higher permeability
To better estimate the rising of the temperature front, we simulate the combined effects of a very narrow inlet and conduit FC with an increase in permeability of FC in the unrest phase (Fig. 6).
To model a narrower conduit, the computational grid is refined with respect to the steady-state condition (Additional file 1: Fig. S1) by reducing the size of the cells within and nearby the conduit with a minimum value of 10 cm near the axis of symmetry. The values of temperature, pressure, and mixture composition computed in the steadystate phase for a 100 m-wide inlet and conduit (Additional file 1: Fig. S1) on the previous grid nodes have been interpolated to the new grid nodes to be used as initial conditions in the unrest phase. To obtain significant temperature variations along the conduit and at the surface, the width of the inlet and that of the conduit need to be decreased. Since the real size of the inlet is unknown, we performed different simulations by gradually decreasing its width from 50 m to 20, 10, and 5 m. The narrowest width of 5 m is useful to simulate the flow propagation in the permeable fractured zones, that prompts the development of the fumaroles field, characterized by permeability increases, affecting the crater area. Simultaneously, the permeability of the conduit is increased by one or two orders of magnitude with respect to the steady-state phase. The most significant results for both HP and LP configurations are obtained when the FC conduit, that is two orders of magnitude more permeable, is 5 or 10 m-wide. Figure 6a-d shows the temperature and pressure variation distributions when the inlet and the conduit are 5 m-wide. The temperature change distributions show different patterns in the two configurations. In the HP configuration, the most significant temperature changes, which reach values above 200 °C up to a depth of approximately 1 km, are attained along the conduit up to the sea level. In the LP configuration, on the contrary, the size of the conduit does not significantly affect the pattern of temperature variations, where large variations appear in the deepest and most superficial areas, as we have already found when only the permeability is increased in the larger conduit (Fig. 3b).
The distributions of the pressure variations reveal a behavior similar to that already noted in the previous simulations. The pressure variations, albeit higher in the inlet zone of the LP configuration, reach lower surface values than in the HP configuration, generating smaller surface ground deformations (Fig. 6g-j).
Horizontal deformations of a few millimeters are obtained, which reach their maximum amplitude at about 1.14 km and 2 km away from the axis of symmetry, in the HP and LP configurations, respectively. The vertical deformations reach their maxima of about 2 and 1 cm in the HP and LP configurations, respectively, close to the hydrothermal source and then vanish moving away from the axis of symmetry. Also, the gravity changes are greater in the HP configuration, in which values slightly higher than 6 µGal are reached right above the main conduit, and then vanish at a distance of about 2 km from the axis of symmetry. In the LP configuration, instead, we observe values slightly less than 2 µGal that vanish at a distance of 2 km from the axis of symmetry.
(See figure on next page.) Fig. 6 Narrowest inlet and higher permeability in the conduit. (a-f) Temperature (∆T) and pressure (∆P) changes and gas saturation ( S g ) during the unrest phase; (g-l) radial distributions of horizontal (U) and vertical (V) displacements and of gravity changes ( g ) on the ground surface over time (from 0.0 to 1.0 year) induced by the injection of 1485 t/d of H 2 O and 726 t/d of CO 2 , at a temperature of 600 °C from a 5 m-wide inlet. The unrest phase was simulated for the permeability configurations HP (left panels) and LP (right panels) after increasing the permeability of the central conduit FC by two orders of magnitude Stissi et al. Earth, Planets and Space (2021)  In the unrest phase, an increase of the average mass flux rate is attained as the conduit is narrowed for the LP configuration. For conduit widths of 5, 10, and 20 m, the values of the average mass flux rates within the conduit area are 1.4588, 0.8840, and 0.5236 t/d/m 2 . For the HP configuration, we have obtained 4.0290, 2.9589 and 3.4252 t/d/m 2 .

Simulating a cyclic crisis event
The main outcomes of the parametric study give us some indications on the main response of the volcano edifice and provide useful clues to simulate one of the main crisis episodes, from January 2001 to April 2005, that occurred at Vulcano Island in the last 30 years. This period was characterized by an increase in the vent temperature, that started in November 2004 and culminated with a temperature peak of 398 °C in January 2005 (Diliberto 2011). It was also accompanied by an increase in CO 2 output, in the flow content, and in the local seismicity Diliberto 2011). Considering the temperature measurements ( Fig. 7) recorded on La Fossa cone at the fumarole F5AT (Fig. 1), which is the most representative fumarole of deep processes (Diliberto 2011), and the CO 2 concentration (Granieri et al. 2006;Selva et al. 2020), we simulate three different periods: (i) a linear trend, from January 2001 to 5 November 2004, preceding the crisis period, characterized by a linear temperature increase that reached the maximum value of 351 °C; (ii) the crisis period, from 6 November 2004 to 13 January 2005, in which the temperature peak was 398 °C; (iii) the period from 14 January to April 2005, when temperature decreased to about 375 °C. The initial conditions in January 2001 were achieved by simulating the injection of a mixture composed of 528.5 t/d of water and 68 t/d of CO 2 (Granieri et al. 2006), which constitutes only 5% molar fraction of the entire mixture (Selva et al. 2020), injected from a 100 m-wide inlet. When the steady conditions are reached, considering the effects of the size  (Tables 3, 4). The temperature measurements are extrapolated from Diliberto (2011) of the inlet and permeability of the main conduit on the modeled parameters, we appropriately tune both parameters to mimic the temperature pattern observed by Diliberto (2011). The inlet and conduit are restricted to 5 m, to allow the hot fluid to reach the ground surface rapidly by simulating the rising of the flow along a highly fractured zone. A gradual increase of the CO 2 flow rate from 68 t/d to 300 t/d and a stepwise increase in the molar concentration of CO 2 from 5 to 15% and in the permeability of the central conduit for the LP configuration from 1.25 × 10 -12 m 2 to 3.87 × 10 -12 m 2 are necessary to simulate the temperature linear increase preceding the crisis period. We tested both the HP and LP configurations. A good fit is found only for the LP configuration. Using the HP configuration under several model parameters, high surface temperatures are not reached even in the steady state. Therefore, HP does not allow to mimic the temperatures measured by Diliberto (2011).
In the next 69 days, flow rate and molar concentration of CO 2 were increased to 330 t/d and 20%, respectively, as well as the conduit permeability to 8.2 × 10 -12 m 2 to obtain a gradual increase in temperature until the peak of about 398 °C, as measured on 13 January 2005. Successively, a decrease in the CO 2 flow rate (223 t/d), in the molar concentration of CO 2 (10%) and in the permeability of the central conduit (stepwise decrease from 3 × 10 -12 m 2 to 1.5 × 10 -12 m 2 ) is simulated to reproduce the lowering of the temperature observed from 14 January to 30 April 2005 (Tables 3, 4).
In Fig. 8, we report the temperature and pressure variations reached in January 2005 and in April 2005. Temperature variations affect the central conduit throughout the period. In January 2005, when the temperature reaches the maximum value, the system is depressurized (due to the sharp increase in rock permeability) and is characterized by positive surface temperature variations, which contribute to small thermo-elastic ground deformation (Fig. 8g, i). In April 2005, the central conduit cools down and a buildup of pressure (about 5 MPa) affects the inlet zone, which contribute to small ground deformations (Fig. 8h, j). The temperature contributes to the deformations close to the source. Gravity changes values around 0.3 µGal and 0.2 µGal are modeled in January 2005 and April 2005, respectively ( Fig. 8 k-l).

Discussions and conclusions
The complex interplay among fluid transport, geochemical, and mechanical processes controls the dynamics of volcano-hydrothermal systems and regulates the spatial-temporal evolutions of geochemical and geophysical observables at the ground surface. Quantitative analysis of this interplay can be deduced through numerical modeling of fluid transport and rock deformations. Knowledge of model parameters is of paramount importance to achieve reliable simulation results. Among the parameters, permeability is the most sensitive but the least known. While over the decades, the large-scale imaging of elastic/anelastic rock properties has become feasible with tomographic studies, estimates of hydrological properties are still limited to single-point measurements in wells and rock samples. Since it is not feasible    Tables 3 and 4 to estimate permeability distributions over the entire volcanic edifice, numerical simulations can aid in exploring the dynamic response of the hydrothermal systems under several configurations and highlight those that are in agreement with the observational constraints.
With this in mind, we performed a parametric study at Vulcano. Thanks to the availability of long-term records of temperature, CO 2 emissions rates, and geophysical signals (i.e., seismicity, ground deformation, and gravity), Vulcano is an ideal site to acquire quantitative insights into the dynamics of volcano-hydrothermal systems. The comparisons between the different simulations allow discovering the dominant role of different permeability configurations that best explains the fluid flow dynamics and mechanical response, and also provides a better understanding of the hydrological properties. The simulation results highlight that the two investigated configurations at lower (LP) and higher (HP) permeability engender small ground deformation, which are in general less than 2 cm in the vertical and horizontal components. These outcomes are supported by the field measurements performed in the periods spanning the last crises. The buildup of pressure, since it is confined at depth and around the conduit, prevents large deformation at surface. Moreover, unlike other volcanohydrothermal systems, where a time delay is observed between the increase in seismic activity and the increases in gas emission at the surface (Peiffer et al. 2018), at Vulcano Island, a simultaneous increase in the occurrence of seismo-volcanic events and fluid emission has been observed during the recent crises period (Milluzzo 2012). In fact, this behavior indicates that the gas escape toward the surface is efficient, and hence, despite the presence of an impermeable layer (Madonia et al. 2019) emulating a sort of shallow cap rock, there is not an efficient sealing. This is due to the presence of permeable fractured zones in the crater area, which avoids strong pressure buildup at shallow depth. Also, the lack of significant deformation and volcano-tectonic events support this hypothesis. Although HP and LP configurations produce deformation of comparable amplitudes, differences of one order of magnitude or more arise when gravity changes are investigated (Fig. 2). Except for the HP configuration with a higher permeability in the main conduit (Fig. 3), small gravity changes of a few µGal, that are, however, beyond the accuracies of the modern state-of-the-art instruments (Schäfer et al. 2020), are obtained for almost all the simulations. This is mainly due to the high-temperature regime reached with the steady state that generates a wide gas saturation area around the main conduit. Further increases of temperature cannot engender significant changes in the fluid density. In the gravity model, we disregard density changes induced by internal deformation and internal density boundary displacements (Bonafede and Mazzanti 1998;Bonafede and Ferrari 2009;Currenti et al. 2008Currenti et al. , 2014, since the deformation is too small to produce significant related gravity changes. The lack of gravity variations in the periodic campaigns carried out at Vulcano Island (Di Maio and Berrino 2016) seem to support the LP configuration. Moreover, only LP is capable of generating fast temperature variations at the surface as observed by continuous temperature monitoring (Diliberto 2011). Therefore, the outcomes of the simulations suggest that permeability values of the geological formations above the Vulcano Island hydrothermal system range between 10 -15 and 10 -16 m 2 , falling within the range of low and intermediate values, as defined by Ingebritsen et al. (2010). These values are also in agreement with the permeability deduced from production tests performed on the VP1 well (AGIP).
Episodic gas emission and temperature transients observed at Vulcano Island at various time scales (Diliberto 2011; Granieri et al. 2006) suggest that transient variations in permeability may occur in the highly fractured crater area. Generally, the permeability, that can vary up to 17 orders of magnitude in common geologic media, tends to change over time Ingebritsen et al. 2010). Transient variations in permeability have been reported both in laboratory experiments and in field observations (1 order of magnitude over 1-10 years; Ingebritsen and Manning 2010). Different processes, both physical and chemical, can indeed cause sudden changes in permeability. Increases in temperature under constant effective stress, for example, can lead to the closure of the fractures (Polak et al. 2003), while hot fluids circulation in the host rock, causing dissolution and subsequent mineral precipitation, can alter this physical property Nicolas et al. 2020) resulting in a decrease of the permeability. By contrast, faulting and fracturing, produced by local and regional tectonics, can create new flow paths as well as an increase in pore pressure associated with the migration of fluids (Farquharson et al. 2017;Peiffer et al. 2018) and have the potential to progressively increase the rock permeability over time.
Even if we do not explicitly search for the processes responsible for the permeability transients, our attempt to simulate the cyclic episodic unrest at Vulcano island (Diliberto 2011) shows the important role of permeability transients within the highly fractured crater area in controlling fluctuations in gas emission and temperature. Therefore, the observed fluctuations in degassing rate and temperature changes may be alternatively explained by permeability transients and not only by increases in the flow injection rate at depth. In this study, we have shown that rapid temperature changes at