Nonlinear Processes in Geophysics A 2-D FEM thermal model to simulate water flow in a porous media : Campi Flegrei caldera case study

Volcanic and geothermal aspects both exist in many geologically young areas. In these areas the heat transfer process is of fundamental importance, so that the thermal and fluid-dynamic processes characterizing a viscous fluid in a porous medium are very important to understand the complex dynamics of the these areas. The Campi Flegrei caldera, located west of the city of Naples, within the central-southern sector of the large graben of Campanian plain, is a region where both volcanic and geothermal phenomena are present. The upper part of the geothermal system can be considered roughly as a succession of volcanic porous material (tuff) saturated by a mixture formed mainly by water and carbon dioxide. We have implemented a finite elements approach in transient conditions to simulate water flow in a 2-D porous medium to model the changes of temperature in the geothermal system due to magmatic fluid inflow, accounting for a transient phase, not considered in the analytical solutions and fluid compressibility. The thermal model is described by means of conductive/convective equations, in which we propose a thermal source represented by a parabolic shape function to better simulate an increase of temperature in the central part (magma chamber) of a box, simulating the Campi Flegrei caldera and using more recent evaluations, from literature, for the medium’s parameters (specific heat capacity, density, thermal conductivity, permeability). A best-fit velocity for the permeant is evaluated by comparing the simulated temperatures with those measured in wells drilled by Agip (Italian Oil Agency) in the 1980s in the framework of geothermal exploration. A few tens of days are enough to reach the thermal steady state, showing the quick response of the system to heat injection. The increase in the pressure due to the heat transport is then used to compute ground deformation, in particular the vertical displacements characteristics of the Campi Flegrei caldera behaviour. The vertical displacements range from 1 cm to 10 cm in accordance with the mini uplift, characterizing the recent behaviour of the caldera. The time needed to move fluid particles from the bottom to the upper layer (years) is compatible with the timing of the mini uplift.


Introduction
In an area where both volcanic activity and geothermal resources are present, a complete understanding of the thermofluid dynamic processes characterizing a viscous fluid in a porous medium is required.Thus, the convective heat transfer process becomes strongly relevant in the framework of volcanic risk definition and geothermal exploration.
A region where both the phenomena are present is the Campi Flegrei area (Fig. 1).Campi Flegrei is a caldera complex located in the Campanian plain region of south-central Italy, 15 km west of the city of Naples (Barberi et al., 1991).The plain has been volcanically active for at least the last 50 000 yr with the activity concentrated in Campi Flegrei, at Vesuvius and on the islands of Procida and Ischia.The main caldera at Campi Flegrei has a diameter of 12-15 km (Rosi and Sbrana, 1987;Lirer et al., 1987;Orsi et al., 1996), and its rim is thought to have been formed during a catastrophic eruption, about 39 kyr ago ca. which produced a deposit referred to as the Campanian Ignimbrite (Rosi and Sbrana, 1987;De Vivo et al., 2001).The following eruption, called Neapolitan Yellow Tuff, that occurred 12-15 kyr ago (Civetta et al., 1997;Deino et al., 2004) produced a further collapse of the inner part of the pre-existing caldera.The last eruption, testifying the persistent activity of the magmatic system, occurred in 1538, giving rise to an about 130 m spatter cone called Mt. Nuovo (Di Vito et al., 1987, 1999).The magma chemistry suggests a current magma chamber, which is currently very small, with a volume of perhaps only 1 km 3 .The area is also an example of a system dominated by phreatomagmatic eruptions, in which magma interacts with surface water, primarily sea water or deep aquifers beneath the caldera.A strong increase in knowledge on the Campi Flegrei Caldera was reached after the SER-APIS project (Zollo et al., 2003;Judhenerc and Zollo, 2004), through an extensive marine active seismic survey carried out in the Gulfs of Naples and Pozzuoli.Vanorio et al. (2005) detected a region of low Poisson ratio, which implies the presence of fractured, over-pressured gas-bearing formations.
During the early 1980s, Agip realized geophysical, geological and geochemical investigations in the Campi Flegrei area, in the framework of a geothermal exploitation project (Agip report, 1987).The deep drillings performed in the Mofete, San Vito and Licola areas (Fig. 1) showed the existence of strong heat flows with maximum temperature up to 400 • C at 3 km depth, as testified by the intense fumarolic/hydrothermal activity (Rosi and Sbrana, 1987).The thermal waters mainly consist of water steam and carbon dioxide, together with minor amounts of H 2 S, N 2 , CH 4 and rare gases (Chiodini et al., 2003).The relevance of the heat and mass transport at Campi Flegrei is confirmed by gas diffusivity measurements (Chiodini et al., 2001) that show the high intensity of this phenomenon.
Mechanisms accounting for the understanding of both volcanic and geothermal phenomena including ground deformation have generally involved two hypotheses: one purely mechanic, the other of thermo-fluid dynamic nature.The first is related to pressure generation exerted from the intrusion of magma into a magma chamber (Bianchi et al., 1987); the second is related to pressure generation in a geothermal reservoir, due to pore pressure variations and temperature, depending of the presence of an aquifer (Oliveri del Castillo and Quagliariello, 1969;Casertano et al., 1976).During the last years, the idea of the fundamental role of the hydrothermal fluids in triggering the activity at Campi Flegrei is rising.Nevertheless, signatures indicating cap rock formations, which are required to build up pore fluid pressure within reservoirs, have not yet been detected.Thus, many questions still remain unresolved and the mechanisms responsible for the Campi Flegrei activity not well constrained (Vanorio et al., 2006).
The high structural complexity of the Campi Flegrei area joint with always more evidence of a strong interaction between the magmatic chamber and shallow water table, also in quiescent periods, calls for a deep knowledge of the aquifers and magmatic source locations together with a deep characterization of the media's geophysical parameters and of physical properties of the rocks, in particular at high pressure and temperature conditions, as existing in the Campi Flegrei area.
To investigate the influence of the fluids at Campi Flegrei, several authors (Gaeta et al., 1998(Gaeta et al., , 2003;;Castagnolo et al., 2001;Troise et al., 2001;Todesco et al., 2003Todesco et al., , 2004) ) proposed thermo-fluid-dynamic 1-D and 2-D modelling, taking into account the anomalous thermal field that can induce variations in the rock rheology and in the elastic parameters of the propagation medium.The complex equilibria among the porous matrix and its fluid permeants are influenced either by heat transfer or by the coupled fluid flow (Ascolese et al., 1993a, b).Gaeta et al. (1998) (studying the possibility that hydrothermal fluids play a fundamental role in triggering activity at Campi Flegrei) proposed a model for a simplified analytical description of the transport equation within a Campi Flegrei-like geothermal system, modelled as a water permeated porous solid matrix, with temperature and pressure gradients, solving the equation at the stationary state, neglecting the transients following changes in the boundary conditions.They showed that changes of water flow can be very important in the genesis and evolution of unrest crises.Wohletz et al. (1999) developed a series of 2-D conductive/convective numerical models to evaluate the possible magma chamber configurations that predict the present thermal regime at Campi Flegrei, reproduced by the measured thermal gradients in boreholes at Licola, San Vito, and Mofete, reported by Agip (1987).They used a 2-D finite difference code including convective regimes showing that, for all reasonable models, a convective zone, developed above the magma chambers after caldera collapse, is necessary to achieve the high gradients observed by Agip (1987).Castagnolo et al. (2001) showed the coupled effect of mechanic and thermal perturbations induced in geothermal fluids, using a 1-D analytical approach.Troise et al. (2001) presented a 2-D approach to model fields of temperature and pressure in an aquifer, applying a finite difference scheme to solve the basic equation for the steady state.
As mentioned before, a peculiar characteristic recently discovered, with the help of the continuous monitoring network (Troise et al., 2007), is that the general downlift, beginning after the last uplift crisis, is interrupted by small uplift phenomena that seem to produce a reduction on the following downlift rate.With a thermal-fluid-dynamical model, Gaeta et al. (2003) investigated the occurrence of mini uplift episodes, characterized by relatively small positive vertical displacements.Gaeta et al. (2003) analyzed the joint mass and heat transport under pressure and temperature gradients, in a porous medium, generated by deep source changes and/or its coupling with shallower layers, finding a good agreement between the computed uplift and observed data by assuming realistic changes of the Péclet parameter (−5 ≥ Pe ≥ −1), which corresponds to the actual situation in the well, whereas Pe = −12 corresponds to an additional flux from the bottom to justify the observed ground uplift.In our modelling study we found Pe ranging in the interval from Due to the nonlinearity of the balance equations and the complexity of the solid-fluid interactions, all previous models consider simplifying assumptions depending on the major factors controlling the processes under investigation.The most common simplifications are the steady state, the domain dimensions reduced from 3-D to 2-D or 1-D, the fluid incompressibility, the low-Reynolds flow regime and the Boussinesq approximation to neglect the compressibility effects.To enhance the understanding of the phenomenology, we propose a thermo-fluid dynamic modelling, based on a finite element approach.Our model includes a transient phase, not evaluated in the analytical solutions, steady-state condition, a 2-D domain easily extensible to a 3-D domain, and fluid compressibility.In addition, we select a more recent evaluation of the medium's parameters (e.g.Giberti et al., 2006;Peluso and Arienzo, 2007), which represents a more reliable and accurate model for the thermo-dynamical property distribution in the propagation medium.
The proposed model allows us to pursue the objective of this study through the modelling of conductive/convective heat transfer processes with the following aims: (1) to reconstruct the thermal field and obtain information about the average advective fluid velocity by comparing the vertical temperature profiles with the temperature measured in the wells drilled by Agip; (2) to evaluate the elastic effect through vertical ground deformation caused by the pressure variations produced by hot fluid rising, as a possible cause of the bradisismo; (3) the timing of reaching the thermal steady-state condition.

Heat transfer modelling
Previous approaches to model the thermo-fluid dynamics of Campi Flegrei were generally based on simple analytic 1-D or finite differences 2-D models, solved at stationary state.In our case, trying to obtain a more robust basis for the development of a thermo-fluid dynamical model of the caldera unrest and the related geothermal field, the temperature distribution has been evaluated in unsteady state conditions, with reference to an idealised, two-dimensional system, constituted by a rectangular domain (Fig. 2) of porous rocks, permeated by a liquid solution, assumed in mono-phase condition and incompressible, heated from below and with heat exchange V. Romano et al.: A 2-D FEM thermal model to simulate water flow in a porous media through the side walls.Regarding the porosity and the permeability, the porous medium is assumed to be isotropic.
We propose modelling using a thermal source represented by a parabolic shape function and using more realistic and up-to-date (than those previously used) values for the medium parameters, inferred from the recent literature.In particular, thermal conductivity and specific heat have been evaluated by laboratory measurements performed by using core samples from the geothermal wells (Zamora et al., 1994;Yven, 2001;Vanorio et al., 2003) and for the absolute rock density (Russo et al., 2012).We use a parabolic shape to better represent a temperature increase applied at the central part of the bottom, simulating the contact of a localised magma chamber (generally approximated by a sub-spheric volume), which transfers heat to the shallower aquifers.This temperature profile allows to avoid the boundary effects produced by boxcar sources generally used in the previous papers.Recently, the knowledge of the rock's physical properties in the CF caldera has been considerably improved, at least in some particular areas (Dello Iacono et al., 2008) as well as the location of the possible thermal sources (Zollo et al., 2008).We study also some elastic parameters with the objective of an evaluation of ground displacement associated with convective heat transfer.These elastic parameters were computed starting from a 3-D P-and S-wave velocity model obtained from the SERAPIS project.In fact, a 3-D velocity model of the Campi Flegrei caldera has been reconstructed by seismic tomography using passive data collected during the last uplift crises (Vanorio et al., 2005) and active data collected during SERAPIS project.The two obtained models were merged (Satriano et al., 2006) to get a more detailed velocity model for a volume of 16 × 20 × 8 km, centred at the Bay of Pozzuoli on a 250 m regular grid A 2-D finite elements method (FEM) has been used to solve the equations for transient and steady-state temperature.In the Cartesian coordinate system (x = horizontal, y = vertical), the unsteady thermal energy equation (Bird et al., 2002) for a fluid at constant pressure (neglecting any contribution of radioactive decay, chemical reaction and latent heat of crystallization and fusion) can be written as (for symbols notation see Table A1): where T is the temperature, ε is the porosity, K is the equivalent thermal conductivity (Castagnolo et al., 2001), ρ * c * = 4.25 × 10 6 J m −3 K is the volume-averaged thermal capacity of the perfused rock (Yven, 2001), and ρ l c l = 3.15 × 10 6 J m −3 K is the thermal capacity of the permeant (Castagnolo et al., 2001).In Eq. ( 1) ρ * c * ∂T ∂t represents the rate of increase of thermal energy per unit volume, ∇ •(K∇T ) represents the rate of conduction heat transport per unit volume, ερ l c l (v • ∇T ) represents the rate of convective heat transport per unit volume, where v y = V 0 •(1 − y)•2• x − x 2 is the fluid velocity along the vertical axis and v x = V 0 • − 2 3 x 3 + x 2 − 1 6 is the fluid velocity along the horizontal axis.V 0 is the half of the maximum velocity at the caldera centre for each level, except for the top surface where v y = 0. We hypothetised that for the vertical component of the fluid velocity there is a parabolic shape in the horizontal direction characterised by decreasing values at shallower layers, according to the bottom temperature profile.
The V 0 value has been optimised by the best fit, with respect to the temperature measured in the boreholes, performing each simulations with a different V 0 values, ranging in the interval from 10 −3 to 10 −7 m s −1 .
Porosity values have been provided by Ascolese et al. (1993a, b) as inferred from laboratory measurements.The basal temperature (647 K) has been deduced from borehole temperature observations to simulate a possible magma chamber heating the shallower aquifers, while the upper temperature (293 K) is assumed from room conditions.The model geometry reproduces the characteristics of a geothermal system fed by a sub-spherical magma chamber.The domain sizes are L x = L y = 3 km, compatible with the geophysical information that shows no evidence of the presence of magma reservoirs with a size >1 km 3 at shallow depth (<4 km).In fact, the simulated system consists of a porous matrix heated by a small magma chamber refilled by a deeper magmatic storage zone; a chamber as small as the one assumed in this work would not be revealed by present-day tomographies (Zollo et al., 2003).The magma chamber volume is similar to those that likely fed several Campi Flegrei eruptions in the past, as well as to the small solidified magma bodies identified off-shore of the Bay of Pozzuoli.In addition, the aquifer we describe is effectively localised within the first 3 km of depth in Campi Flegrei where the critical point of water is reached (Troise et al., 2001;Gaeta et al., 2003).Thus, the region at which the permeant reaches the critical temperature is considered external to our system.
For this reason, as an initial condition, a temperature profile described by a parabolic function shape with respect to the horizontal coordinate (x) and linear with respect to the vertical one (y) has been selected: As boundary conditions a heat flux at the lateral surface has been assumed.
where h l is the heat transfer coefficient through the lateral boundary, and T o is the room temperature.At the bottom surface, a temperature profile with a parabolic shape with respect to the horizontal axis has been considered, while a heat flux has been assumed through the top surface: where h t is the heat transfer coefficient through the top boundary.
The thermal conductivity has been evaluated, considering saturated rocks, by the product of the thermal conductivities of water and solid matrix (Yven, 2001): Thermal conductivity varies with temperature.Using the relation by Brigaud and Vasseur (1989) for the solid matrix and by Ozbek and Phillips (1979) for the in situ fluid and neglecting the dependence from the salinity, we can take into account this dependence: (Yven, 2001).

Numerical simulations
The finite element method has been applied to simulate unsteady and steady-state temperature distributions of an aquifer subject to a temperature imbalance, due to a possible magma intrusion at a depth able to transfer heat to the aquifer from its top.The FEM approach tries to model the temperature in a geothermal system due to magmatic fluids inflow.Several tests have been performed by varying the parameters inside realistic ranges.In particular, the vertical fluid velocity has been varied in the interval [10 −3 , 10 −7 m s −1 ], trying to find the best value by inversion of the steady-state vertical gradient of temperature with respect to the temperature data measured in the borehole by Agip (1987).
One of the most important modelling issues is the assessment of the appropriate mesh density.Using a coarse mesh can generate inaccurate results.On the other hand, using a fine mesh leads to an increase of run times.To avoid such problems, in our work a sensitivity analysis of the appropriate mesh density was performed.After this analysis the domain was discretized into 32 620 elements of uniform sizes with a density of about 3624 elements per km 2 .
As an example of the results, the temperature field at different times (t = 0, 2, 4, 6, 10 days) and in steady-state conditions is shown in Fig. 3, for a velocity V 0 = 1 × 10 −5 m s −1 .It is evident that few tens of days are enough to reach the thermal stationary conditions.This relatively short time needed to reach steady-state conditions is clearly due to the relevant contribution to the heat transfer carried by the convective term.Figures 4a and b show the conductive and convective heat flow at steady state, where the convective term is the main contribution in the lower half of the model box, while the conductive term becomes predominant at shallower depth.
We compare the simulated vertical temperature with those measured in the wells.The best-fit model parameter V o was obtained by exhaustively searching for the minimum value of the variance between simulated and observed temperature.
In Fig. 5 the actual temperatures measured in wells in areas drilled by Agip at Mofete and San Vito are superimposed to the simulated steady state temperature distributions, having as curve parameter the distance from the centre of symmetry (x/L x = 0.5) of the model box in steps of 0.1.The simulation was performed using V 0 = 8×10 −6 m s −1 (Fig. 5a), while Figs.5b, c, d, e and f show, as examples, the simulations obtained using respectively, with the increase of the convective term of the heat transfer.The comparison between observed and simulated temperature is able to determine that a fluid velocity V 0 ∼ = 1 × 10 −5 m s −1 produces temperatures with the best fit to the real data on average.The Péclet number ranges between −1 and −10, in a good agreement with the values evaluated by Gaeta et al. (2003).
The simulated temperatures provide evidence of the respective rates of advective transport of the different sites, showing the influence of the relative abundance of the permeant fed into the system: abundant at Mofete 1 and 2, poor at San Vito 1, as reported in the log data of the Agip reports.The convective heat transport is able to justify the characteristic temperature gradients existing in the upper parts of the aquifer more steeply than within the caldera.In more detail, the data at San Vito are best fitted by a simulation with a lower advective fluid velocity in comparison with Mofete.This characteristic could be related to the abundance of fluid locally fed into the system and/or to the different degree of the local rock fracturing, both phenomena improving the mass transport.
Due to the relevance of ground deformation in the Campi Flegrei area, it is interesting to evaluate the effect that fluid circulation can have on the strain conditions of the porous medium, involved in the mass transport produced by a temperature gradient.In porous media, thermal and hydrodynamical characteristics have strong influence on the pore pressure.In fact, a pressure variation at the base of the geothermal system progressively migrates by water circulation to shallower depths, with a characteristic time depending on the velocity of the water flux and on the hydraulic parameters of the rocks.Pore pressure contrasts the lithostatic pressure due to the overburden rocks at different depths and causes elastic deformation of the rocks.The fluid circulation enhances the pore pressure ( p), determining a volume increase ( V ), where µ is the compressibility coefficient.
The linear deformation can then be approximated by The compressibility coefficient has been evaluated using the relationship: , where γ is the P-wave velocity, β is the S-wave velocity and ρ is the density.The bulk modulus has been evaluated using the seismic velocity as inferred from the seismic tomography realized during the SERAPIS project (Zollo et al., 2003), and the density has been evaluated from gravity inversion inferred by Russo et al. (2012).
The increase in pore pressure, produced by the fluid circulation and enhanced by the increased temperature at depth, can be evaluated by (Castagnolo et al., 2001) where H = L y is the total thickness of the domain, k is the hydraulic permeability, while < η > is the average dynamic viscosity of the liquid, and V y is the vertical fluid velocity.Using appropriate values for the parameters, we determine a vertical displacement in the range between 1 cm and 10 cm.This range for the vertical displacement agrees very well with the mini uplift characterising the recent behaviour of the Campi Flegrei caldera (Troise et al., 2007(Troise et al., , 2008)).

Discussion and conclusion
The objective of this study was the improvement of the description of the temperature profile within the Campi Flegrei caldera due to conductive/convective heat transfer processes, from both volcanological and geothermal points of view.The caldera has been modelled as a rectangular box formed by a water-permeated, porous, solid matrix.The geometry of the aquifer has been chosen to reproduce the main features of a geothermal system, overlying a localised magma source.
To get a description of the geothermal system, including the transient following the changes in the boundary conditions, the unsteady thermal energy equation has been resolved by a 2-D FEM approach.Because of the high permeability of the rocks, advective heat transport exceeds conduction through the porous matrix, justifying the observed temperature gradients.The vertical temperature profiles have been compared with the temperature measured in wells (Mofete and San Vito) drilled by Agip in the framework of a geothermal exploration.From this comparison a best-fit value for the average velocity of the advective fluid generated by the temperature gradient has been evaluated as V 0 ∼ = 1 × 10 −5 m s −1 .This velocity was found lower in the San Vito area in comparison with the Mofete area, in accordance with the abundance of fluid locally fed into the system, as inferred during the drilling.In fact, analyzing the well logs, Chelini and Sbrana (1987) demonstrated the relatively high abundance of permeant fed into the system at Mofete 1 and 2 and the poor presence of fluid at San Vito 1.This behaviour can be related to the presence of a more developed aquifer in the Mofete area.Alternatively, the different velocities found between the Mofete and San Vito wells can be related to the spatial distribution of the rocks' permeability, due to the different degree of rock fracturing.This suggests that the use of a 3-D numerical model, instead of a 2-D axisymmetric structure, should be preferred at Campi Flegrei (D'Auria et al., 2010).Our model and parameterisation are easily extensible to a 3-D domain.
The pressure variations generated by the hot fluid rising determine an elastic effect through ground deformation (vertical displacement) that was computed ranging between 1 cm and 10 cm.This range for the vertical displacement is well in agreement with the mini uplift characterizing the recent behaviour of the Campi Flegrei caldera.
Because we assume a depth of about 3 km for the bottom of the geothermal system, the inferred velocity of advection implies about 9.5 yr for deep fluids to reach the surface.In terms of its order of magnitude, this time scale is comparable with the duration of the last crisis (about 2-3 yr), confirming that thermo-fluid dynamical effects can contribute strongly to the mechanical behaviour of the Campi Flegrei caldera.
The relatively short time needed to reach thermal steadystate conditions (few tens of days) could generate a strong onset of the elastic deformation able to produce sharp vertical displacements.Therefore, the transient of a strong gradient in the ground deformation could be difficult to observe with a daily rate of measurements of ground deformation.This consideration highlights the problem of calibrating the time scales between phenomena (ground deformation)

Fig. 1 .
Fig. 1.Sketch of Campi Flegrei caldera, located in Southern Italy, inside Campania region.The red symbols indicate the areas where the geothermal wells of Mofete, San Vito and Licola have been drilled by Agip.

Fig. 2 .
Fig. 2. Schematic representation of the computational domain representing an idealized two-dimensional geothermal system.Domain sizes are L x = L y = 3 km.Boundary conditions are specified.

Fig. 4 .
Fig. 4. Conductive heat flow at steady state (A), convective heat flow at steady state (B).Vectors represent intensity and direction of heat flow.