Caldera’s Breathing: Poroelastic Ground Deformation at Campi Flegrei (Italy)

Ground deformation at Campi Flegrei has fuelled a long-term scientific debate about its driving mechanism and its significance in hazard assessment. In an active volcanic system hosting a wide hydrothermal circulation, both magmatic and hydrothermal fluids could be responsible, to variable degrees, for the observed ground displacement. Fast and large uplifts are commonly interpreted in terms of pressure or volume changes associated with magma intrusion, while minor, slower displacement can be related to shallower sources. This work focuses on the deformation history of the last 35 years and shows that ground deformation measured at Campi Flegrei since 1985 is consistent with a poroelastic response of a shallow hydrothermal system to changes in pore pressure and fluid content. The extensive literature available for Campi Flegrei allows constraining system geometry, properties, and conditions. Changes in pore pressure and fluid content necessary to cause the observed deformation can then be calculated based on the linear theory of poroelasticity. The predicted pore pressure evolution and fluid fluxes are plausible and consistent with available measurements and independent estimates.


INTRODUCTION
Caldera deformation easily stirs scientific debates about the actual source driving ground displacement and the crucial implications for hazard assessment. The Campi Flegrei caldera is a unique case history, with an impressive data set that describes ground deformation over a century, including three significant uplift phases starting in 1950, 1969, and 1982(Del Gaudio et al., 2010. Each uplift episode lasted a couple of years or less and together totaled more than 3.5 m of ground displacement ( Figure 1A). After 1985, a prolonged subsidence phase began, interrupted by occasional, minor, and short-lasting uplift events in 1989, 1994, and 2000. Uplift resumed after 2005 initially at a slow rate, which then increased in 2012 ( Figure 1B), causing a shift in alert level from green to yellow (Attention). This uplift phase is still going on, uninterrupted, at the time of writing (2021). Ground deformation data have been interpreted as the product of shallow magmatic sources (Berrino et al., 1984;Amoruso and Crescentini, 2011;Trasatti et al., 2011;D'Auria et al., 2015), with the possible contribution of hydrothermal or aqueous fluids (Casertano et al., 1976;Bonafede, 1991;Chiodini et al., 2003;Amoruso and Crescentini, 2011). Numerical modeling showed that hydrothermal circulation could cause some deformation, but without justifying the order of magnitude of the large uplift phases observed in the 50's, 70's, and 80's (Todesco et al., 2004;Rinaldi et al., 2010). On the other hand, subsidence is often explained in terms of poroelastic relaxation, considering the extensive degassing that characterizes the caldera. Many authors have interpreted the subsequent uplift phase in terms of heating and pressurization of the hydrothermal system (Chiodini et al., 2011;Aiuppa et al., 2013;Chiodini et al., 2015;Moretti et al., 2017). The ongoing uplift phase proceeds at a rate that is similar to the rate of subsidence, and, according to some authors, this implies that deformation mechanism and mechanical properties remained unchanged (Moretti et al., 2018). According to Amoruso et al. (2014), both a deeper (magmatic) source and a shallower (hydrothermal) one contributed to the deformation history of the caldera. Modeling of magma sill intrusion confirmed that subsidence might also arise from magma emplacement (Macedonio et al., 2014), but without ruling out the role of hydrothermal fluids, which may enhance sill propagation by hindering magma cooling (Amoruso and Crescentini, 2019) and contribute to both uplift and subsidence phases .
The potential involvement of the hydrothermal system in the recent unrest phases is also suggested by the significant gas discharge that characterizes the caldera. It mostly occurs through soil and fumaroles within the Solfatara crater, roughly at the caldera center, and its eastern flank (Pisciarelli). The presence of fumaroles at this location has been known for centuries, and measurements of gas temperature begun as early as the end of the 19th century ( Figure 2) (Sicardi, 1941). Temperature values recorded during this entire time span (shown in Figure 2) display an average of 154.17°C, with a standard deviation of 8.27°C. Extreme values were recorded at the beginning of the period when both minimum and maximum values were measured (122°C, in 1899, and 215°C, in 1935). These first measurements may be unreliable, as we cannot establish the measurement method or the exact location of the measurements. However, even considering only data collected after 1980, the average temperature and standard deviation do not change significantly (154.21 and 7.71°C, respectively). An histogram of the data (Figure 3) reveals two modes at 162°C and 148°C. These values correspond to the average temperature recorded at the two main fumaroles within the Solfatara crater (BG and BN, respectively). The temperature difference between the two fumarolic vents was recently explained in terms of a very shallow cooling effect associated with condensed steam flowing from BG into BN fumarolic conduit (Gresse et al., 2018). All these FIGURE 2 | Gas temperature measurements carried out within the Solfatara crater since the end of the 19th century as reported by Sicardi (1941); Sicardi (1970);Dall'Aglio et al. (1972); Cioni et al. (1984); Martini (1986);Caliro et al. (2014). Data and references for older measurements, described by Sicardi and published in Italian jounals, are provided as Supplementary Materials. considerations suggest that the large hydrothermal system feeding the fumaroles acts as a buffer, preventing large temperature fluctuations.
Discharged gases are mostly steam and carbon dioxide. Data on gas composition are available since the 80's and can be expressed in terms of CO 2 /H 2 O ratio. The evolution of gas composition was exhaustively described and discussed by Caliro et al. (2014) and more recently by Chiodini et al. (2021). After a first period characterized by fluctuating behavior, a steady increasing trend has begun in 2000 and continues today, leading to the maximum observed values. Estimates of gas discharge rate became available in 1998 (Cardellini et al., 2017). Measurements of diffuse degassing through the ground of the Solfatara crater revealed for the first time the magnitude of gas discharge at Campi Flegrei, where more than a thousand tonnes of carbon dioxide are released every day. CO 2 degassing also displays some fluctuation through time, associated with changes in the size of the degassing area. The measurements constrain only a fraction of the actual gas discharge in the area, as they only refer to carbon dioxide, neglecting other gas components, and do not account for the gas discharged at the fumaroles. The first available estimates of fumaroles contribution (Aiuppa et al., 2013) suggests that fumaroles could contribute an additional 40%.
Aim of this work is to test the hypothesis that the post-1985 deformation is a poroelastic response to change in pore pressure and fluid content within the hydrothermal reservoir. Literature data are used to constrain the position, size, geometry, and properties of the rock volume undergoing poroelastic deformation. As the Campi Flegrei caldera has been extensively studied for decades, a large amount of information is available to build a solid conceptual model, that will be described in the following sections, together with the choice of rock properties and system conditions applied here (Table 1). Once the conceptual model is established, it is possible to calculate the changes in pore pressure and fluid content required to justify the observed ground displacement as a poroelastic response. Results confirm that poroelastic deformation is a plausible mechanism to explain the ongoing unrest phase and that calculated pressure history and fluids fluxes are consistent with independent evidence and the system's current understanding.

CONCEPTUAL MODEL
The linear theory of poroelasticity relates the deformation of fluid-saturated rocks to pore pressure and fluid content through a set of parameters that include the properties of both the porous rocks and the fluid. Some simplifying assumptions allow for an analytical solution: the poroelastic source is assumed to have a simple, cylindrical geometry, and to correspond to a very permeable (∼10 −13 m 2 ) volume of rocks, where pressure changes propagate quickly and can be considered uniform over the source volume. The crude simplification is partially justified by the impressive gas discharge that has been measured in the caldera over the past 20 years. While this approach certainly prevents a detailed description of the system evolution through time, it allows us to explore the orders of magnitude of pressure changes and fluid discharge that would be involved in case of poro-elastic deformation. In the following, available literature data are used to define such volume and its properties. Once the poroelastic source is defined, it can be used to compute the changes in pore pressure and fluid content required to justify the observed vertical displacement.

The Hydrothermal Reservoir and Its Fluids
The reconstruction of the subsurface structures is based on geological considerations, constrained by data and samples from geothermal wells drilled in the area. The collapsed caldera is characterized by alternating layers of tuff and marine sediments, intersected by lava and overlaying thermometamorphic rocks (Piochi et al., 2014). Further information is retrieved from inversion of seismic data, showing that the seismogenic region is confined between a high-velocity basement at 3 km depth, and a shallower cap-rock found at depths between 1 and 2 km (Vanorio et al., 2005). Vanorio and Kanitpanyacharoen (2015) showed that this shallow layer is characterized by high values of the elastic moduli (as determined from seismic velocity data), and ascribed this effect to metamorphic decarbonation of marine sediments, causing a particular hardening of the rocks at these depths. The presence of this seismic horizon was later confirmed by Calò and Tramelli (2018) who identify a 500 m thick layer separating the deeper portion of the system, where magma may be present, from the shallower part. Similar conclusions were reached by Akande et al. (2019), who identified a shallow domain, where aquifers surround a hydrothermal reservoir, that is separated from a deeper domain by a cap-rock that these authors locate at a slightly shallower depth. Joint analyses of seismic velocity, quality factors and scattering provided information on the nature of fluid saturating the rocks, suggesting the presence of shallow aquifers (at depths <800 m) and of a two-phase reservoir extending to the cap-rock (Calò and Tramelli, 2018). According to these authors, the cap-rock is extensively fractured, allowing magmatic volatiles to reach the hydrothermal reservoir and mix with meteoric water. Similar attenuation anomalies were also identified by Akande et al. (2019). The presence of a widespread hydrothermal reservoir where magmatic volatiles mix with meteoric waters is consistent with geochemical evidences, such as the lack of acid species (Chiodini et al., 2021) or the rather constant emission temperature. The observed changes in gas composition have been interpreted in terms of a variable contribution of magmatic volatiles through time (Chiodini et al., 2003). As the systems evolves through periods of increased or reduced magmatic input, distribution of fluid phases may change, due to steam condensation or water evaporation (Todesco, 2009). Numerical modelling of hydrothermal circulation showed that the presence, size and position of gas-dominated, liquiddominated or two-phase regions depends on the combination of several features, each of whom may change through time. Rock anisotropy, and in particular the distribution of permeability were shown to affect significantly the phase distribution of the system . But any process affecting the pressure and temperature distribution within the hydrothermal system may shift the relative proportions of gas and liquid. Phase distribution has critical implications for geochemical modelling, a technique used to infer pressure and temperature conditions based on the composition of gases discharged at the surface (Chiodini et al., 2021). The assumption of liquid water coexisting with steam implies a trend of increasing of equilibrium pressure and temperature (Chiodini et al., 2015(Chiodini et al., , 2021, whereas the relaxation of such assumption leads opposite result (Moretti et al., 2017).
The relative proportion of gas and liquid within the system also bears important consequences in terms of poroelastic behavior of the fluid saturated rocks.
Further information to constrain phase distribution derives from magnetotelluric data, collected along a 5.6 km long profile by Siniscalchi et al. (2019). These authors developed a large scale resistivity model that highlights the presence of a shallow conductive, water-saturated layer above a high resistivity anomaly (100-300 Ω · m) underneath the Solfatara crater, extending to a depth of 3,000 m. The anomaly was interpreted as the well developed gas plume that feeds fumarole emissions. Figure 4 illustrates a conceptual model that summarises the evidences described above: the hydrothermal reservoir is simplified as a vertical cylinder 1,500 m thick and 1,000 m wide, whose top is located at a depth of 500 m. While consistent with the caldera setting, where horizontal layers of pyroclastic materials are limited by ring faults, this basic geometry allows a simple mathematical description of the poroelastic problem, described in the Method section. The shallow hydrothermal reservoir is fed by hot magmatic gases that Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 702665 separate from melt at depth and rise through fractures in the caprock layer. Meteoric water permeates the shallower layer, while a gas-dominated region likely occupies the inner portion of the reservoir, where magmatic volatiles and steam of meteoric origin may mix. In this conceptual model, the deformation observed at the surface (observation point O) is due to changes in pore pressure associated with the fluids entering or leaving the hydrothermal reservoir.

Rock Properties
The poroelastic response to fluid and pressure changes largely depends on the choice of rock properties. Various estimates exists in the literature, for different rock types at Campi Flegrei, based on measurements of both on outcropping rocks and on core samples from geothermal wells (Vanorio et al., 2002;Giberti et al., 2006;Heap et al., 2014;Vanorio and Kanitpanyacharoen, 2015). Elastic properties change with lithology, in situ conditions (depth, thus confining pressure, temperature, water content), and may be affected by compaction and mineral alteration (Heap et al., 2014;Vanorio and Kanitpanyacharoen, 2015). Even though a detailed description of the spatial and temporal evolution of rock properties goes beyond the aim of the present work, a choice of reasonable values can be guided by available data. Large porosity and shallow depth favour low values of the bulk modulus. Temperature may also lower the bulk modulus, depending on specific rock fabric and state of alteration. Vanorio et al. (2002) published data on samples of Neapolitan Yellow Tuff (NYT) and Campanian Ignimbrite tuff (CI). For samples with higher porosity, ϕ 0.46, measurements yielded bulk modulus, K 5 GPa and Poisson's ratio, ] ≈ 0.3. The sample with lower porosity (ϕ 0.33) showed a higher bulk modulus (K 10 GPa) and lower Poisson's ratio (] 0.25). Estimates given by Giberti et al. (2006) include: rock porosity, ranging from ϕ 0.48, for samples at or near the surface, to ϕ 0.1 or lower, for samples at depths greater than 1,000 m) (saturated) bulk modulus ranges from K 6.3 GPa, at the surface to K 47.3 GPa at depths. Dry values are lower, ranging from 1.60 to 29 GPa; Dry Poisson's ratio ranges from 0.28 (at surface) to 0.14 (at depths). Wet values are higher, spanning from 0.43 to 0.27, with increasing depths. According to Manconi et al. (2010), the shear modulus in the area may vary from 2 to 14 GPa, with a Poisson's ratio ranging from 0.2 to 0.4. The corresponding bulk modulus ranges from 2.7 to 65.3 GPa.
Heap and coworkers (Heap et al., 2014) highlighted the difference between dynamic estimates of elastic moduli, based on seismic data, and static ones, obtained from stress-strain relations. These authors suggest that static values (typically lower than dynamic ones) are more appropriate to describe the mechanical behaviour of rocks undergoing volcanic deformation. According to their measurements, dynamic values of Poisson's ratio range from 0.28 to 0.33 (lower value for dry conditions) and shear modulus ranging from 1.97 to 3.20 GPa (lower values for dry conditions). These values correspond to a bulk modulus K ranging from 4.27 to 6.93 GPa. Static values yield lower values of shear modulus (from 0.66 to 0.81 GPa, under an effective pressure of 5 MPa). Considering a Poisson's ratio of ] 0.3, the corresponding static bulk modulus reaches even smaller values, from 1.43 to 1.75 GPa. More recently, a detailed compilation of laboratory measurements carried out on volcanic rock samples was published by Heap et al. (2020), who highlighted the importance of porosity in determining the value of Young's modulus. These authors used experimental data to constrain an empirical relation that allows to upscale laboratory estimates of elastic moduli to larger length scales, where macroscopic fractures and heterogeneities are expected to lower their values. According to this approach, tuffs characterized by large porosity (ϕ 0.45) and with a Poisson's ratio of ] 0.3 are expected to have a bulk modulus as low as K 0.6 GPa. This value increases to K 4.5 GPa if porosity drops to 0.15. Porosity data from core samples collected during deep drilling at Campi Flegrei provide a guide to porosity changes with depth, with values of ϕ 0.33 and 0.22 at depths of 810 and 1,420 m, respectively (Giberti et al., 2006). The corresponding values of the bulk modulus range from K 1.4 to K 2.9 GPa (Heap et al., 2020).
Based on these data, rock properties are chosen as follows ( Table 1): porosity is set to the value observed in a core sample taken at 810 m (ϕ 0.33); the bulk modulus K 1.5 GPa, the Poisson's ratio ] 0.3. The Biot-Willis coefficient, α, is set equal to 1, implying that the deformation is entirely due to change in pore volume, while the solid grains do not deform. Figure 4, the post-1985 observed ground deformation can be considered as an unbounded problem domain, where the poroelastic displacement at the free surface is generated by changes in pore pressure or fluid content within a cylinder, placed at shallow depth. As a first order approximation, we neglect the role of temperature changes. The solution can be borrowed from the analogous gravitational problem (Wang, 2000, p. 109).

Given the conceptual model in
Considering a cylinder of thickness h and radius a, located at a depth d (Figure 4), the vertical displacement u z at the surface, expected at the origin point O due to a uniform pressure change p 0 can be expressed as (Wang, 2000): where ] is the Poisson's coefficient and c m is the uniaxial poroelastic coefficient. Under conditions of uniaxial strain, this coefficient is the ratio of axial strain to pore pressure changes and can be expressed as a function of the Biot-Willis coefficient α, the bulk modulus K and the drained Poisson's coefficient ]: .
Pressure change p 0 can therefore be obtained from the observed vertical displacement u z : Displacement can also be expressed in terms of the uniform change of fluid content ζ 0 : Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 702665 where c is the loading efficiency. c can be expressed as a function of the Skempton's parameter B, relating the pore pressure changes to the change in applied stress under undrained conditions, and the undrained Poisson's ratio ] u : The undrained Poisson's ratio ] u can be expressed as a function of drained Poisson's ratio ν, the Biot-Willis coefficient and the Skempton's parameter: The change in fluid content ζ 0 can be therefore computed as: (4)

PORE PRESSURE CHANGES
Changes in ground elevation (Figure 1, measured at the Pozzuoli harbour at different times since 1985 (Del Gaudio et al., 2010;, are used here to compute the corresponding changes in pore pressure, according to Eq. 2. Recent deformation is described by GPS data , which provide a more detailed temporal representation of ground displacement. Being interested in the general trend, this time series was sampled in order to keep the time interval between subsequent measurement comparable to the average time interval (≈154 days) that characterized the older data set, based on levelling measurements. Higher sampling rates do not alter the general outcome, but produce noisier results. The two data sets overlap over the period from 2000 to 2010, when they show comparable trends. In the following, values from the most recent GPS data set are considered after the year 2000.
Results are plotted in Figure 5 where pressure changes corresponding to each measured displacement are shown in green ( Figure 5A), while the overall pressure evolution is shown in blue ( Figure 5B). Open symbols refer to the data set published by Del Gaudio et al. (2010), while solid symbols are based on De . During the considered time period, pressure declines, with only minor increments corresponding to the mini-uplifts (Gaeta et al., 2003). The maximum pressure drop since 1985 is −4.26 MPa and is achieved in November 2004. The following uplift was associated with a continuous pressure build up, leading to a pore pressure increment of 2.38 MPa in September 2018.
The choice of system geometry may affect this result. Figure 6 compares the pressure evolution computed with different choices of reservoir size and position. Larger reservoirs require smaller pressure changes. The depth of the reservoir's top also affects the pressure change necessary to produce the observed displacement, with deeper reservoirs requiring greater pressure variations. The estimates for overall pressure drop associated with the post-1985 subsidence vary from −2 to 7 MPa. Only in the case of a very thin vertical cylinder with a radius of 500 m, the pressure drop required to achieve the observed displacement reaches −11 MPa.

CHANGES IN FLUID CONTENT
Under the same assumptions described above, displacement data can be used to compute the corresponding changes in fluid content (Eq. 4). In this case, we need to define the value of the Skempton's coefficient B, relating the pore pressure changes to the change in applied stress, under undrained conditions. Most applications of the linear theory of poroelasticity describe systems where the pore fluid is incompressible (typically, liquid water), for which B l 1 is a reasonable approximation. If we consider liquid water as the fluid filling the reservoir pores, the change in fluid content corresponding to the observed deformation range from −2.6 × 10 -4 to 1.9 × 10 -4 . Negative values refer to water leaving the reservoir, during caldera subsidence, while positive values refer to water entering the hydrothermal system, causing its inflation. The change in fluid content is a dimensionless quantity, normalized with respect to a reference volume. If we consider the reference volume to be the initial pore volume available in the hydrothermal reservoir, we can compute the corresponding changes in fluid volume that range from −4 × 10 5 to 3 × 10 5 m 3 , with an average value of −1.7 × 10 4 m 3 . Knowing the fluid density and the time intervals separating two consecutive measurements, we can then transform these changes in fluid volumes into flow rates ( Figure 7A). Values are computed considering a constant liquid density of 960 kg/m 3 (corresponding to near-boiling temperature at atmospheric pressure), the liquid flow rates required to cause the observed displacement range from −24.27 to 45.23 kg/s.
In the case of Campi Flegrei, however, evidence suggests that a significant portion of the hydrothermal reservoir is occupied by gas. To compute the change in fluid content for a gas-saturated system, we need to define the Skempton's coefficient and the fluid density, which are both depending on the reservoir evolving pressure. Considering the pressure evolution depicted in Figure 5, we can arbitrarily set the initial reservoir pressure at 4.5 MPa and apply the computed pressure changes to reconstruct the pressure history. Following Wang (2000), the value of the Skempton's coefficient for highly compressible fluids can be approximated as: where α is the Biot-Willis poroelastic expansion coefficient, K f and K are the fluid and rock bulk moduli, respectively, and ϕ is rock porosity. For highly compressible fluid K f ≪ K and, under isothermal conditions, K f corresponds to the reservoir evolving pressure. The B g values corresponding to the computed pressure evolution range from 4.7 × 10 -4 to 8.83 × 10 -3 , with smaller values corresponding to lower reservoir pressures. The corresponding changes in fluid content ranges from −0.1 to 0.075. Again, we can convert these dimensionless quantities into volumes by considering the initial pore volume available in the hydrothermal reservoir. The Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 702665 fluid volume change in case of a gas-saturated rock ranges from −1.56 × 10 8 to 1.17 × 10 8 m 3 . To compute the corresponding mass flow rates we need to introduce the time interval among different measurements and to specify the gas density. The density of the gas mixture depends on both evolving reservoir conditions and gas composition. As a first order approximation, we are considering a constant reservoir temperature that we can set at 250°C. Pressure history is based on the evolution shown in Figure 5, and considering an initial value of 4.5 MPa. We can consider the gas as a mixture of steam and carbon dioxide, whose densities at the considered system conditions can be computed based on data from the NIST Chemistry Webbook (Linstrom and Mallard, 2020). Gas composition (Caliro et al., 2014) is then used to set the relative proportion of steam and carbon dioxide at different times and to compute the corresponding density of the gas mixture. Obtained values range from 0.70 to 30.36 kg/m 3 , with lower values corresponding to the minimum pressure (0.24 MPa reached in November 2004). The corresponding gas flow rates range from −101.64 to 66.30 kg/s ( Figure 7B). The mini-uplifts observed after 1985 are clearly marked by positive increments in fluid before the year 2000. As expected, a gas-dominated system requires higher discharge rates than a liquid-dominated system.
As discussed in previous sections, it is reasonable to assume that at least in some portions of the hydrothermal reservoir underneath Solfatara both gas and liquid phases coexist.
The change in fluid content in a two-phase system depends on the relative proportion of gas and liquid that can be expressed in terms of gas volume fraction, or gas saturation (S g ). For a twophase mixture (where gas and liquid move together), we can compute the Skempton's parameter and the fluid density as a function of S g , based on their values for the single phases: B mix S g B gas + (1 − S g )B liq ρ mix S g ρ gas + (1 − S g )ρ liq Figure 8 shows the effect of the increasing the amount of gas in the mixture on the corresponding flow rates. When liquid water dominates (S g < 0.5), a smaller change in fluid content is required to achieve the observed deformation, with computed flow values that almost correspond to those obtained with a complete liquid saturation. However, as the gas fraction in the reservoir increases, the amount of fluid that needs to leave or enter the reservoir to justify the deformation changes significantly, with the maximum absolute value ranging from 44.56 kg/s with S g 0.95, to 69.73 kg/s when S g 0.99, to 121.17 kg/s, when the system is gas saturated.
This comparison shows that even a small fraction of water almost halves the flow rate required to cause a given deformation. It should be noted that the gas fraction values indicated in Figure 8 refers to the proportion of gas and liquid in the fluid mixture that leaves or enter the reservoir. Such proportion does not necessarily reflect the phase distribution within the reservoir. Even when both phases are present, gas and liquid are easily separated by their different density, viscosity, relative permeability, and by the effects of capillary forces. Generally speaking, gas is more mobile and more likely to move than liquid. When this is the case, the poroelastic response only depends on the properties of the phase that moves in or out from the reservoir (and on the volume available to that phase).
Computed changes in gas content can be used to estimate the overall amount of gas that transited through the reservoir since 1985. Open dots in Figure 9 show the total mass of fluid leaving  Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 702665 the reservoir at different times. Solid dots illustrate the gas increment that is being stored inside the same reservoir, during uplift.
Considering the fumarole gas composition as representative of the ratio between steam and gas everywhere in the system, we can compute the total amount of carbon dioxide that passed through the reservoir. The total mass of carbon dioxide is shown in Figure 10, where the total mass of gases is also shown. Over the time interval considered, the total amount of carbon dioxide transited through the system is 10.99 × 10 9 kg while the total mass of gas is equal to 27.52 × 10 9 kg.

DISCUSSION
The post-1985 ground displacement observed at Campi Flegrei include a 20-years-long subsidence phase followed by a 15-yearslong uplift. Unlike previous uplift episodes, this deformation history could be entirely due to changes in pore pressure and fluid content in a shallow hydrothermal reservoir. This work tested this hypothesis, by introducing some simplifying assumption about the poroelastic source of deformation. Data available from recent scientific literature was used to constrain the size, geometry, depth, and properties of the shallow hydrothermal reservoir. Observed ground displacement was then used to infer the driving changes in pore pressure and fluid content. Despite the simplifications (in particular, the uniform pressure changes throughout the source) the calculations presented here provide a first order estimate of quantities involved and allows some interesting considerations. In the following, results are compared with available data; then the role of temperature changes is examined; and finally the implications of a shallow, poroelastic source are discussed.

Comparison With Available Data
According to the calculation presented above, the observed subsidence phase corresponds to an overall pore pressure drop of −4.26 MPa, while the subsequent uplift phase requires a pressure build up of about 2.4 MPa ( Figure 5). An independent assessment of the conditions of the hydrothermal reservoir is based on geochemical equilibria of the CO 2 − CO − CH 4 − H 2 0 − H 2 system. Based on fumarole gas composition, inferred equilibrium pressure range from 3 MPa, during the 1982-84 unrest (Cioni et al., 1984;Chiodini et al., 2001), to 0.3-0.7 MPa, at the end of the subsidence (Chiodini et al., 2001). Further data for the subsequent period led to equilibrium values ranging from 1.7 to 3.7 MPa (Chiodini et al., 2011), with higher values corresponding to periods of uplift. More recently, Chiodini et al. (2021) extended the range from 2.7 to 6 MPa (or from 3.7 to 7.8 MPa, depending on the choice of redox buffer). The order of magnitude of these independent estimates is comparable with the pore pressure variations depicted in Figure 5. Equilibrium pressure values do not represent the conditions of the entire reservoir, but only those of a shallow equilibration zone (Chiodini et al., 2021). The comparison with these data however shows how the estimated pressure history ( Figure 5) is plausible and spans over a similar range. Different choices of reservoir geometry may account for slightly higher values of equilibration pressure (up to almost 8 MPa), with greater pressure changes being favoured by a thinner, smaller, or deeper reservoir.
Measured ground displacement were also used to compute the corresponding changes in fluid content. Such changes depend on the choice of fluid phase and properties. Maximum changes in fluid volume during the subsidence range from −4 × 10 5 m 3 (for liquid water) to −1.56 × 10 8 m 3 (for dry gas). Slightly lower values characterize the subsequent uplift phase, ranging from 2.9 × 10 5 m 3 , in case of liquid, to 1.17 × 10 8 m 3 , in case of gas. If the pore fluid is liquid water, the post-1985 subsidence requires a maximum discharge rate of 25 kg/s of water, whereas a greater maximum flow rate (45 kg/s) is needed to justify the ongoing uplift phase. These are small values compared to the shallow groundwater fluxes estimated for the nearby city of Naples, where discharge rate toward the sea is quantified in approximately 300 kg/s and aquifer recharge is of the order of 500 kg/s (Ducci and Sellerino, 2015).
When a gas phase permeates the reservoir, greater quantities are involved: when a mixture of steam and carbon dioxide is  Open symbols refer to all gas, solid symbol to CO 2 .
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 702665 considered, the maximum flow rate required to justify the post-1985 subsidence is as high as 102 kg/s, while the uplift corresponds to a fluid injection rate of 66 kg/s. Contrary to the case of a liquid saturated system, in this case greater flow rates are required during the subsidence than during the uplift, reflecting the strong dependence of gas density on evolving pore pressure. Calculated gas flow rates can be compared with available measurements of CO 2 diffuse degassing (Cardellini et al., 2017). CO 2 flow rates can be obtained from the calculated total gas flux ( Figure 7B), based on gas composition at different times (Caliro et al., 2014). In Figure 11 Figure 11 shows the ongoing uplift phase can be generated by flow rate of the same order of magnitude of measured fluxes. These figures depend on the choice of fluid phase, and even minor fraction of liquid water may significantly reduce the fluid change required to cause the observed displacement ( Figure 8). The comparison between calculated and observed flow rate suggest that the fluid change required to justify the observed deformation is reasonable within the Campi Flegrei volcanic system, and is consistent with available information.

Temperature Effects
Presented results do not account for temperature changes. Chiodini et al. (2021) inferred a temperature increment of almost 50°C in 10 years since 2010, based on geochemical equilibria and under the assumption of a two-phase system. We can make a rough estimate of the effects of such a temperature change on the deformation of the porous matrix. Considering a coefficient of thermal expansion of the order of 10 -5 , such a temperature increment would result in a volumetric strain of 5 × 10 -4 . If we consider the hydrothermal reservoir in Figure 4 as the the initial deforming volume, the corresponding volume change would be of the order of 7.85 × 10 5 m 3 . Assuming that the entire deformation takes place only along the vertical direction, the resulting vertical displacement would be of the order of 0.25 m, i.e., approximately half of the deformation observed at Campi Flegrei since 2010. Temperature changes also affect fluid properties, with higher temperatures lowering the fluid density and the corresponding flow rates. In a dry gas system at the condition of interest, a temperature increment of 50°C would reduce the flow rate by 15%. Greater effects are expected if phase changes are involved. Temperature changes are also going to affect rock hydraulic and mechanical properties, modifying the way in which both the porous medium and the fluid would respond to changes in pore pressure. A full description of the entire phenomenon requires the use of a fully coupled thermo-fluid-mechanical model Akande et al., 2021;Nespoli et al., 2021).
Regardless to the exact magnitude of the displacement caused by thermal expansion, heating during uplift is certainly responsible for a fraction of the ongoing deformation, reducing the pressure build up and fluid change necessary to obtain a given displacement. The pressure history and flow rates computed neglecting temperature effects are already comparable with independent estimates and can be considered plausible for the hydrothermal reservoir at Campi Flegrei. Accounting for thermal effects would require even smaller changes to explain the observed deformation, leaving the poroelastic deformation a very reasonable candidate to explain the post-1985 deformation history at Campi Flegrei.

Implications of a Poroelastic Source of Deformation
The role of hydrothermal fluids as potential drivers of poroelastic deformation is not new and was proposed in the past to explain ground deformation observed at Campi Flegrei (Casertano et al., 1976;Bonafede, 1991;Todesco et al., 2004;Battaglia et al., 2006;Rinaldi et al., 2010;Akande et al., 2021;Nespoli et al., 2021) and in different contextes (Hurwitz et al., 2007;Currenti et al., 2017;Mittal and Richards, 2019). The simple calculation presented above shows that poroelastic deformation is a feasible mechanism to explain both subsidence and uplift observed at Campi Flegrei during the last 35 years. Unlike previous unrest episodes, characterized by fast and large ground displacement, the post-1985 deformation is consistent with the inflation and deflation of a shallow hydrothermal system, changing pore pressure and fluid content. The idea was already proposed by several authors (Chiodini et al., 2011;Aiuppa et al., 2013;Chiodini et al., 2015;Moretti et al., 2017Moretti et al., , 2018, but the hypothesis was not tested before. This mechanism is simple, plausible, and consistent with available information on the caldera volcanic system. In the simplest configuration, the hydrothermal reservoir is fed by volatiles of magmatic origin and releases fluids at the surface through fumaroles and diffuse degassing. Even in such a simple system, the pressure history of the reservoir does not simply depend on the actual amount of gases available (i.e., on the size and conditions of magma degassing at depth). The overall evolution is controlled by the balance between gas entering the hydrothermal reservoir and gas discharged at the surface. Even in the presence of a constant gas supply, deflation or inflation may occur depending on how fast gas is discharged at the surface. After 1985, deflation may have occurred because magmatic gas recharge declined or stopped altogether, while discharge continued. However, subsidence could also have occurred with a steady, or even increasing, magmatic input, provided it was accompanied by a faster fluid removal from the reservoir.
The implications of ongoing unrest at Campi Flegrei depend on the source that drives the deformation, causing the switch between subsidence and uplift. Many processes operate within the magma chamber, influencing the exsolution process and affecting the rate, amount, and composition of released volatiles. However, the widespread hydrothermal system that lies between the degassing magma body and the surface is an active interface that controls and modulates the transit of fluids toward the surface (Mittal and Richards, 2019). Here, the interaction with liquid water may influence the upward propagation of magmatic volatiles, because of phase interference (related to relative permeability and capillary pressure effects) and the possible dissolution of gas species. Fluid-rock interaction may also intervene to modify the porosity and permeability of the porous medium. These shallow features affect the ascent of magmatic volatiles, thereby controlling the pressure history and the amounts of fluids within the hydrothermal reservoir.
The pattern defined by gas leaving and entering the reservoir well describes a deflation-inflation cycle (Figure 9) that seems to be related to the total amount of fluids that can be stored within the reservoir: the rate at which the gas enters the reservoir (solid dots) only increases after fluid discharge (open symbols) stops growing. At this time, the reservoir has reached its maximum deflation and lowest pore pressure, and more fluid can enter the reservoir. As the figure shows, the total gas input has been approaching the total mass discharged in the last few years since 1985. The gradual shift toward conditions existing in 1984 was already highlighted with concern, as long-term stress and damage accumulation may favor crust failure and the onset of eruptions (Kilburn et al., 2017). We do not know the actual pressure within the magmatic system and the amounts of fluids available to feed the reservoir at this time, as we do not know the actual state of the rock and the maximum pore pressure it can withstand, so there is no way to predict how the system will evolve. However, if fracturing does not occur, the hydrothermal reservoir could act as a valve, allowing fluid entrance and inflation only as long as its pore pressure is lower than the pressure of external fluids (magmatic volatiles or shallower meteoric fluids). If the displacement reached in 1985 represents an upper limit for reservoir inflation, the corresponding pore pressure might have been high enough to prevent further input into the reservoir. Under these conditions, a subsidence phase began, thanks to a fast fluid discharge at the surface and high pressure in the reservoir. As fluids are progressively removed, reservoir pressure drops, and a minor but increasing fluid input modulates its decline. Both the rate of fluid injection and fluid discharge change through time, governed by the evolving pressure, until they reach similar values, and the deformation pauses, as occurred between 2000 and 2005. From this time on, fluid input is enhanced by low reservoir pressure. Eventually, the rate at which fluids enter the system overcomes the rate at which they are discharged at the surface and uplift resumes. Once pore pressure inside the reservoir overcomes that of surrounding rocks, no more fluid cans enter, and subsidence could begin once again. The results presented here do not exclude the possible involvement of magma in driving the recent deformation at Campi Flegrei. Magmatic and hydrothemal fludis could both participate, to variable extents, in determining the evolution of past and ongoing unrest. The present study however shows that hydrothermal fluids alone can be responsible for significant ground deformation. Repeated cycles of uplift and subsidence without eruption already happened in the caldera history (Morhange et al., 2006;Isaia et al., 2009;Todesco et al., 2014), yet without conclusive explanations. A mechanism of this sort could work on different spatial and temporal scales, providing a plausible justification for these observations. Performed calculations and available data also provide some clues on the amount of carbon dioxide that entered the hydrothermal reservoir during the considered period. When calculated changes in fluid content are close to zero, the amount of gas that enters equals the amount of gas that leaves the reservoir and corresponds to the measured discharge ( Figure 11). The CO 2 input from 2000 to 2005 must have been of the order of 12 kg/s. In the following years, uplift resumed and, toward the end of 2006, the calculated fluid input in the reservoir reaches 5.6 kg/s, accompanied by an observed fluid discharge of 14.4 kg/s. To sustain both surface degassing and the increment in fluid content within the reservoir, at least 20 kg/s of carbon dioxide must have entered the reservoir. Later on, the measured flux exceeds 32 kg/s, while computed fluid increment reaches a maximum of 30 kg/s: more than 60 kg/s of carbon dioxide (corresponding to 5,184 tonnes per day) must have entered the reservoir. The increasing amount of carbon dioxide points toward an increased magmatic degassing or increased efficiency of decarbonation reactions in the cap-rock.
In total, more than 27 million tonnes of gas passed through the hydrothermal reservoir, modulating the ground displacement ( Figure 10). Approximately 11 of these were carbon dioxide. Assuming that all gas has a magmatic origin, these amounts can be used to constrain the size of the magma body degassing at depth. Considering a magma with a density of 2,700 kg/m 3 and initial content of total dissolved volatiles of 2.5% and a CO 2 content of 1%, the volume of magma necessary to exsolve the calculated amounts of gases is slightly less than half cubic km (4.07 × 10 8 m 3 ). Considering a higher initial gas content (5% total gas, 2% of CO 2 ) the corresponding magma volume halves (0.20 × 10 8 m 3 ). These values are reasonable for a volcanic system like Campi Flegrei, where similar volumes were estimated for magma intrusions before the Monte Nuovo eruption (Di Vito et al., 2016). Smaller magma volumes would be obtained considering a two-phase scenario (that would reduce the amounts of fluids involved) or accounting for the presence of a non-magmatic source of fluids.
The simple calculations presented above showed that reasonable changes in pore pressure and fluid content might cause ground deformation of the same order of magnitude of the observed displacement. Pressure history and fluid content may be affected by phase transition, chemical reactions, or changes in rock porosity and permeability that control the available pore volume and the rate at which fluid can propagate and exchange heat with the porous rock. A detailed description of all these effects requires a more sophisticated approach, accounting for the coupled thermo-hydro-mechanical processes in the natural system, from the magma chamber level to the surface (Mittal and Richards, 2019). Indeed, all these changes may generate measurable signals, including ground displacement and changes in gas flow rate and composition. These signals bear very different implications than those associated with an increased magmatic degassing due to a shallow intrusion. The analysis of recent seismicity at Campi Flegrei confirmed the connection between deformation, shallow seismicity, and hydrothermal activity (Giudicepietro et al., 2020). On the contrary, geophysical signals do not appear to reflect magma movement. The possibility that the ongoing unrest is related to a shallow poroelastic deformation should be taken into account and carefully examined.

CONCLUSION
The results and considerations described above allow us to draw some conclusions: • The post-1985 ground deformation at Campi Flegrei can be described as a poroelastic response to changes in fluid content and pore pressure in a shallow hydrothermal reservoir • The subsidence phase is consistent with an overall pressure drop within the reservoir of −4.26 MPa, while the following uplift required a pore pressure build-up of 2.38 MPa; • Corresponding changes in fluid content were transformed into flow rates that reached maximum values at the beginning of the subsidence (−102 kg/s) and the end of the studied period (66.30 kg/s). • Computed changes in pore pressure and fluid content are plausible, comparable with independent estimates and observations, and describe a system evolution that is consistent with the well-established conceptual model of the caldera. • During the last two decades, surface degassing increased during uplift periods implying a growing input sustaining both reservoir inflation and gas discharge.
• A shallow hydrothermal origin for the ongoing deformation may explain non-eruptive cycles of subsidence and uplift, driven by the balance between magmatic input and fluid discharge at the surface.
This work shows that, unlike previous unrest episodes, poroelastic deformation alone may provide a simple and plausible explanation for the volcanic unrest after 1985 that is not accompanied by significant evidence of magma movement. Further studies are necessary to address temperature effects, constrain the role of phase changes, and account for the rock's heterogeneous and evolving properties. Results presented above highlight the role of the hydrothermal system as a potential source of significant unrest episodes. A comprehensive volcanic hazard assessment should account for this possibility. This study confirms once again the importance of careful and continuous monitoring of the hydrothermal fluid circulation.

DATA AVAILABILITY STATEMENT
All datasets analyzed in this study can be found in the references cited in the text. Older temperature data that appear on Italian documents are provided as Supplementary Material.

AUTHOR CONTRIBUTIONS
MT confirms being the sole contributor of this work and approved it for publication.

FUNDING
This research was supported by INGV, Sezione di Bologna.