Gravimetric Constraints on the Hydrothermal System of the Campi Flegrei Caldera

The Campi Flegrei caldera (Italy) has been undergoing unrest over the past five decades including episodes of rapid ground deformation, seismicity, and variations in gas emissions. Hydrothermal fluids and gases are released most vigorously in the central sector of the caldera at the fumarolic fields of Solfatara volcano and Pisciarelli. We conducted a high‐precision gravity survey coupled with inverse modeling to image the shallow (<2‐km depth) structure of the hydrothermal feeder system. Results indicate the presence of three low density bodies beneath Pozzuoli, Astroni volcano and the Solfatara/Pisciarelli fumarolic fields. The first two are inferred to be sealed hydrothermal systems trapped beneath impermeable cap rock, while the latter depicts a plume‐like geothermal feeder system reaching the surface via a combination of Solfatara's maar‐diatreme structure and the intersection of NW‐SE and NE‐SW trending regional faults. The density contrasts of the reservoirs from background values are best explained by a multiphase mixture of caldera fill containing a secondary and interconnected void volume fraction of between 0.2 and 0.3 that hosts a vapor volume fraction ψv of between 0.38 and 1 and a liquid volume fraction ψl fraction of between 0 and 0.62. This work highlights the control of volcano‐tectonic structures on fluid movement in the shallow crust of hydrothermally active volcanic systems undergoing sustained or periodic unrest.


Introduction
Volcanic unrest is often characterized by anomalous seismicity, gas emissions, and surface deformation and is usually attributed to subsurface magma movement (Sparks, 2003). Volcanic calderas have complex subsurface structures resulting not at least from the vertical collapse of a preexisting volcanic edifice and often host both extensive hydrothermal and magmatic reservoirs (Gottsmann & Battaglia, 2008). Hydrothermal systems are a complex interface between magma reservoirs and the surface (Todesco, 2008) and not only produce measurable unrest signals but also modulate geophysical and geochemical signals from underlying magma reservoirs (Chiodini et al., 2002(Chiodini et al., , 2016Gottsmann & Battaglia, 2008;Ingebritsen et al., 2010;Todesco, 2008).
Campi Flegrei caldera (CFc) is a well-documented restless caldera where the separation of the signals from magmatic and hydrothermal sources has not been trivial (Troise et al., 2019). Solfatara volcano and neighboring Pisciarelli (Figure 1) host the main surface features of the hydrothermal system at CFc and are located ∼2.5 km to the NE of Pozzuoli, the center of ground deformation over the last 50 years  of unrest. Our current understanding of the structure and dynamics of the hydrothermal system at CFc is informed predominantly by geochemical constraints, geophysical data, and resulting models (Bruno et al., 2007;Caliro et al., 2007;Troiano et al., 2019): A multiphase plume of vapor and liquid fuelled by the interaction of magmatic and meteoric fluids at depth rises to feed fumaroles and mud pools at the surface (Chiodini et al., 2015). While significant subsurface density variations are expected from this model, gravity data have not been used to contribute to the understanding of the shape and size of the shallow-seated (<1.5 km) hydrothermal feeder system. Here, we present results from a new gravimetric survey of the central sector of the CFc including a high-resolution gravity survey of Solfatara volcano coupled with data inversion to image the density structure of the uppermost part of the hydrothermal system.

CFc Structure and Recent Unrest History
CFc is a ∼13-km-wide volcanic caldera in the Campanian Plain near Naples, Italy (Vitale & Isaia, 2014), formed by two major vertical collapses at ∼40 ka (Giaccio et al., 2017) and ∼15 ka (Deino et al., 2004). ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2019JB019231
Key Points: • Hydrothermal plumbing of the central sector of the Campi Flegrei is controlled by volcano-tectonic structures • Three shallow-seated (<1.5 km depth) hydrothermal feeder systems are imaged at Pozzuoli, Solfatara/Pisciarelli volcano, and Astroni volcano • Low densities of feeder systems are explained by porous caldera fill material with between 0.38 and 1 vapor volume fraction and between 0 and 0.62 liquid volume fraction in secondary void space Postcollapse eruptive activity over the last 15 ka generated at least 70 eruptions, mainly concentrated in the central eastern sector of the caldera (Bevilacqua et al., 2015;Smith et al., 2011). More than 20 eruptions occurred in the epoch of activity from 5.8 to 3.8 ka forming landmarks such as Astroni volcano and Agnano caldera (Isaia et al., 2009;Smith et al., 2011), with the latest magmatic eruption in 1538 CE creating Monte Nuovo volcano (Barberi et al., 1984). The vertical collapses and long-term ground deformation have divided the caldera floor into a block structure . The dominant fault trends within the caldera are NW-SE and NE-SW (Di Florio et al., 1999;Isaia et al., 2015;Vitale & Isaia, 2014) in addition to caldera ring faults (Berrino et al., 2008(Berrino et al., , 1991Gottsmann et al., 2006;Zollo et al., 2003).
The caldera fill is composed of intercalated lava flows, pyroclastic material, and marine and continental sediments (Piochi et al., 2014;Rosi & Sbrana, 1987). Gravity data depict the fill as a broad (∼6-km wavelength) low density anomaly (Barberi et al., 1991;Capuano et al., 2013). Drilling encountered a zone of thermometamorphic rocks below the fill at depths between 2.5 and 3.1 km, several small igneous intrusions, high thermal gradients of 100-170 K km −1 , and locally raised isotherms in the central eastern part (Piochi et al., 2014) of the caldera.
Episodes of rapid uplift and subsidence have been recorded at CFc since Roman times (Bellucci et al., 2006;Parascondola, 1947), though subsidence has been the dominant long-term trend (Barberi et al., 1984). Recent rapid ground uplift occurred in 1969Recent rapid ground uplift occurred in -1972Recent rapid ground uplift occurred in and 1982Recent rapid ground uplift occurred in -1984.5 m near Pozzuoli (Barberi et al., 1984). Bucking a trend of slow ground subsidence since 1989, numerous episodes of mini-uplift have occurred producing deformation on the order of centimeters accompanied by seismic activity (Chiodini et al., 2010(Chiodini et al., , 2015 followed by an episode of sustained uplift beginning in 2005. Fumarolic flow rate, discharge temperature, and seismicity increased at Solfatara volcano and Pisciarelli from 2006 onward, and an increased magmatic contribution was inferred from the composition of the fumarolic gases (Chiodini et al., 2016;Giudicepietro et al., 2019). The cause of the rapid meter-scale uplifts is still controversial with magmatic and hydrothermal sources or a mix of both identified as contributors (Troise et al., 2019). Many authors suggest that the more recent mini-uplift events have an exclusively hydrothermal origin (Amoruso et al., 2014;Chiodini et al., 2015;D'Auria et al., 2011;Gottsmann et al., 2006;Manconi et al., 2010). NW-SE faults are reactivated during uplift and subsidence and may be important pathways for the upward movement of gas and magma to the surface (Vilardo et al., 2010)

Solfatara Volcano and Its Fumarolic Fields
Solfatara volcano is the most thermally active part of the caldera . It releases 10 times more thermal energy than the conductive heat flux across the entire caldera floor (Chiodini et al., 2001). Hydrothermal gases and fluids are released most vigorously at Solfatara volcano's crater floor and its easter inner wall as well as at the Pisciarelli fumarolic field located on its NE flank ( Figure 1) (Caliro et al., 2007). The 100°C isotherm resides only a few hundred meters below the surface of Solfatara's crater (Piochi et al., 2014).
The main hydrothermal features within Solfatara's crater include a mud pool named La Fangaia and two main fumaroles named Bocca Nuova and Bocca Grande (see Figure 1). Detailed geological mapping by Isaia et al. (2015) identifies Solfatara's structure as a maar-diatreme. High-angle normal faults characterize the crater edges, and fault planes are predominantly NW-SE. Two lava domes reside in the NE and S crater walls. The crater itself is embedded in the older structure of the Agnano-Monte Spina Complex, and fumaroles and hydrothermal activity are concentrated in the fault zones and their intersections, where highly fractured rocks act as preferred pathways for fluid ascent.

The Hydrothermal Plume
The presence of a hydrothermal plume beneath Solfatara volcano was first suggested by Cioni et al. (1984), who, based on geochemical data, proposed that dry steam separates from a geothermal liquid at 236°C in a highly fractured zone to feed the fumarole fields. Further compelling evidence was presented by Chiodini et al. (2015) based on fumarole geochemistry, CO 2 flux, water table heights, seismic velocity, and interferometric synthetic aperture radar (InSAR) data, as well as by thermo-hydro-mechanical modeling (Coco et al., 2016;Todesco & Berrino, 2005). The numerical models require the multiphase flow of ascending hot fluids (H 2 O and CO 2 ) from depth through a porous medium to reproduce measured fumarole emissions, ground deformation, and gravity changes. The current conceptual model of the plume suggests that rising magmatic gases flash hydrothermal liquids in a deep "mixing zone" and form a gas plume, which ascends to the surface (Caliro et al., 2007;Chiodini et al., 2015). A summary of relevant geophysical and geochemical surveys of Solfatara volcano and the wider CFc is given in Table 1.

Data Acquisition
We performed a new static gravity survey from 8-12 July 2015 using a Scintrex CG-5 Autograv gravimeter (serial number: 572) in tandem with a TOPCON HiPer Pro Dual-Frequency GNSS base and rover system. The survey area encompassed the highly urbanized central sector of the CFc (Figure 2) and contained a total of 85 benchmarks laid out in two different spatial networks.
Benchmarks within Solfatara crater were ordered in a dense irregular grid with a minimum spacing of 17 m and included a local gravity control point. The remainder of the benchmarks were spaced more widely along the roads around Solfatara volcano with an average spacing of 1 km and a maximum spacing of 2.5 km. The different spacing permitted us to investigate the expression of the hydrothermal plume at Solfatara volcano at a similar scale to several local geophysical studies while also exploring the spatial distribution of the hydrothermal system across the central sector of the CFc ( Figure 2). Both the GPS reference receiver and the main gravity base station were located near Monte Nuovo, and all gravity measurements were tied to this reference. The entire network covered an approximate area of 36 km 2 , and the precision of repeat measurement was ±15 μGal (average of 12 cycles of 30 s long readings of 6 Hz raw data at each benchmark). Urban noise led to an average standard error of individual gravity measurements of ±8 μGal, which is a factor of between 3 and 5 higher than usually attainable during quiet conditions.
We recorded GNSS data for 5-20 min at 1 Hz at the survey benchmarks using a roving receiver/antenna unit. The base receiver/antenna unit recorded continuously at 1 Hz during the survey period. The derived precision of the benchmark locations was generally under 0.05 m in the vertical and better than 0.04 m in the horizontal after baseline processing of the benchmark locations against the base station, which in turn was processed against three permanent reference stations of the local INGV Osservatorio Vesuviano Permanent Global Navigation Satellite System (GNSS) network and three regional International GNSS Service (IGS) references (NOT1, MAT1, and MEDI).

Gravity Data Reduction and Correction
The objective of a static gravity survey is to obtain information about the subsurface density distribution. The magnitude of gravity at any point is influenced by latitude, elevation, topography of the surrounding terrain, Earth and ocean tides, subsurface density variations, and instrumental drift (Telford et al., 1991). Raw gravity data are therefore composed of several contributions and require careful corrections to obtain the component reflecting subsurface density variations only, known also as the Bouguer anomaly (BA). Earth tides and instrumental drift are removed first to obtain the observed gravity (g obs ) from which the BA can be obtained.
where g n is the normal gravity, FAC is the free-air correction, BS is the Bouguer slab correction, and TC is the terrain correction. A detailed description of the data reduction is given in the supporting information. Note. Question marks (?) indicate undisclosed or unclear survey dates.

Data Detrending
The regional BA is controlled by both shallow and deep-seated density distributions. Long-wavelength features (e.g., spatial variations in deep-seated bedrock thickness) must be removed to reveal the local BA caused by shallow-seated structures. We calculate a regional gradient of 0.86 mGal km −1 with a strike of N50°E from the regional BA data and derive the linearly detrended residual anomaly (LRA) data for further investigation. Our regional trend compares to a regional gradient of 0.5 mGal km −1 and a strike of N35°E presented by Cassano and La Torre (1987) who use a much larger and wider-spaced data set. The Topex gravity data (Sandwell et al., 2013) with a spatial coverage and average station spacing matching more closely with our survey give a regional gradient of 0.2 mGal km −1 with a strike of N37°E.
To test the robustness of our results, we detrended our data using the regional trend from the Topex data set. While the amplitudes of the resultant anomalies of course change, the presence and location of the main anomalies remain. Therefore, even using the lowest quoted regional gradient, we obtain model results that are reproducible and robust.
Removing a linear trend may not be appropriate in structurally complex areas such as collapse calderas to investigate anomalies associated with a shallow-seated hydrothermal system. Large-scale gravity surveys at CFc have consistently shown a negative gravity anomaly associated with low-density caldera fill (Berrino et al., 2008;Capuano et al., 2013;Nunziata & Rapolla, 1981). We explore the effect of the fill on our data by constructing a forward model based on the most recent gravity data presented by Capuano et al. (2013) and borehole density data (Barberi et al., 1991;Piochi et al., 2014). The caldera fill is simulated by stacked spheroids within a cylindrical volume of 2 km thickness and 3 km in radius with a density contrast of −300 kg m −3 , centered offshore of Pozzuoli (see supporting information Figure S5). Capuano et al. (2013) suggest that the uppermost part of the caldera fill contains remnant high-density feeder systems, as well as postcollapse lava flows and domes. We therefore set the top of the model at a depth of 1 km and subtract the simulated caldera fill anomaly from the regional BA data. We thus obtain a second local anomaly: the caldera fill detrended residual anomaly (CRA).
Values for the LRA and CRA anomalies are reported relative to the base station at Monte Nuovo. Secondary anomalies of Solfatara volcano have their values calculated from average background values of Solfatara's crater floor. Relative values are provided so that the anomalies are comparable across the two detrended data sets.

Total Horizontal Gravity Gradiometry
The first and second horizontal derivatives of Bouguer gravity data are useful to study structural controls on gravity anomalies (Cooper & Cowan, 2008). The first derivative highlights boundaries of buried bodies or faults. The second derivative yields inflection points of the first gradient and reveals absolute maxima/minima, which provide information on the shape of buried bodies or inclination of density interfaces. The total horizontal gradients are obtained from ( 2) where THD1 is the first total horizontal derivative, ∂g/∂x is the change in gravity in the x direction, and ∂g/ ∂y is the change in gravity in the y direction (Cooper & Cowan, 2008). Similarly, ( 3) where THD2 is the second total horizontal derivative (Fedi, 2002).

Inverting the Local BA Data
We invert the resultant Bouguer gravity data (LRA and CRA) to image causative density contrasts at depth using GROWTH2.0 (Camacho et al., 2002(Camacho et al., , 2011. GROWTH2.0 divides the model space into 3-D parallelepiped elements and obtains a 3-D anomalous density model using prescribed (a priori) density contrasts. Inherent nonuniqueness in the inversion is addressed by using a mixed minimization condition, which selects a solution based on least squares model fitness and model smoothness, or the total anomalous mass. Model inputs include the Bouguer gravity data, cell size, the density contrast with background density, and a balance factor. The balance factor determines the complexity of the model of positive and negative density contrasts with high balance factors producing simple models. Densities with too high a contrast produce isolated skeletal bodies and densities with too low a contrast produce inflated and interconnected bodies. Both constraints must be explored to find an appropriately complex model with a plausible density contrast and low autocorrelation (Camacho et al., 2011). As a result, the minimization of residuals is insufficient to establish the suitability of the model for given density contrasts. A mixed minimization procedure that balances the goodness-of-fit criterium (model fitness) with the total anomalous mass (model smoothness) minimization condition is applied to select the optimal model. We explored the model space for each data set iteratively, searching for suitably complex anomalous bodies and low autocorrelations for a low cell resolution and repeatedly increasing the resolution and retesting of the density contrast and balance factor at each iteration. A 50 kg m −3 km −1 increase in background density was implemented to prevent oversizing of anomalous bodies with increasing depth during the inversion.
We tested density contrasts in the range of ±300 to ±600 kg m −3 and balance factors from 10 to 40 (producing in total 119 model solutions) and selected the model with the lowest autocorrelation for given model smoothness. This methodology effectively uses a classic trade-off between model misfit and model simplicity (Gubbins, 2004). The best solution balances a compromise between adequately fitting the data and producing a suitably simple model. While the model with the best goodness of fit has an autocorrelation of 0.06 and a balance factor of 20, it yields an array of skeletal bodies of anomalous densities and does not satisfy the mixed minimization criteria for an optimal solution. Our optimal model of the CRA has an autocorrelation of 0.13 and a balance factor of 40 after 58 iterations, while the optimal model of the LRA has an autocorrelation of 0.14 and a balance factor of 40 after 61 iterations. Details on the inversion procedure and sensitivity tests are given elsewhere (Camacho et al., 2002(Camacho et al., , 2011.

The BA
The amplitudes of all anomalies are orders of magnitude above the uncertainties associated with individual measurements or the terrain correction and are therefore robust indicators of subsurface density variations.

Journal of Geophysical Research: Solid Earth
Figures 3a-3c show the distribution and amplitudes of the regional and local Bouguer anomalies (LRA and CRA). Linear detrending ( Figure 3b) reveals a broad (∼4 km wide) and negative (∼ −6 mGal) anomaly centered northwest of Solfatara volcano. It is composed of three distinct lows near Pozzuoli, Solfatara volcano, and Astroni volcano. The Pozzuoli anomaly is not present in the CRA data (Figure 3c), and the overall negative anomaly is significantly reduced in both wavelength and amplitude (∼0.4 km and ∼ −4 mGal, respectively) and its center shifted toward the north. The gravity lows between Solfatara volcano and Astroni volcano persist. Anomalies at the periphery of the survey are poorly constrained and hence ignored.
When looking in more detail at the Solfatara area, the patterns of the LRA and CRA anomalies are similar but with noticeable differences in the negative amplitudes (from average values of the crater floor) of the respective gravity lows. The main gravity low in the eastern part of the crater is −1.1 mGal in the LRA from an average value of −3.15 mGal and −0.76 mGal in the CRA from an average value of 1.42 mGal (Figures 4a and  4b). In both cases, the lows are located close to the fumaroles of Bocca Nuova and Bocca Grande and extend eastward toward Pisciarelli. There are gravity highs on the north-northeastern and southern edges of the crater in both data sets, but the north-northeastern high is strongest in the LRA data.
While there is a small gravity low (∼ −0.3 mGal amplitude from background levels in the crater floor) in the vicinity of La Fangaia in the LRA data (Figure 4a), this anomaly is only very weak in the CRA data ( Figure 4b).

Horizontal Derivatives
The first horizontal derivative of the LRA (Figure 5a) reveals strong gradients along the northeastern crater wall of Solfatara, around the edge of the low gravity region between La Fangaia and Bocca Grande and more subdued gradients around La Fangaia and elsewhere in the crater. Prevailing NNE-SSW and NW-SE trends are highlighted by the gradients (Figure 5b). The second horizontal derivative suggests similar fault trends (Figures 5c  and 5d). The structural trends obtained from total horizontal gravity gradiometry closely match field observations (Figures 5e and 5f).

Subsurface Distribution of Anomalous Mass
The optimal LRA and CRA models have a balance factor of 40, an autocorrelation of 0.14, and an a priori density contrasts of −450 to +450 kg m −3 ( Figure 6). The models image three main bodies of negative density contrast beneath Pozzuoli, Astroni volcano, and Solfatara volcano. Although it is difficult to directly relate mathematically derived density contrasts with rock density contrast, the optimal density range matches the 1σ range in rock densities about an average of 2,300 kg m −3 encountered in boreholes from Campi Flegrei (Piochi et al., 2014). Despite their inherent nonuniqueness, the models consistently provide robust results on the density variations at depth for different a priori density contrasts (supporting information Figures S6 and S7). The dominant anomalous negative density bodies persist in all inversions, although as expected they become larger and more interconnected with decreasing a priori density contrasts.
The optimal LRA inversion images the Solfatara/Pisciarelli anomaly as approximately 0.5 km wide and extending from close to the surface to 0.8 km below sea level. The anomalous body beneath the Pozzuoli Figure 3. Regional and local Bouguer anomaly maps. (a) shows the regional Bouguer anomaly, while (b) and (c) show the local Bouguer anomaly after detrending the data presented in (a) for two different components: (b) a regional long-wavelength trend to obtain the linearly detrended residual anomaly (LRA) and (c) an additional trend caused by the caldera fill sediments to obtain the caldera fill detrended residual anomaly (CRA). The color bars are scaled to the maximum absolute gravity of the three data sets. Contours are in mGal, and benchmarks are shown in white. The data are overlain on the 10 m DEM and Pozzuoli, Solfatara volcano, and Astroni volcano are marked by letters P, S, and A, respectively. area is 1 km in diameter at its widest and extends from ∼0.5 km to 1.2 km depth. It is slightly elongated in the NNE-SSW direction. The Astroni anomaly is elongated E-W, 1.75 km across its widest point and extends from 0.5 km to 1.4 km depth ( Figure 6).
The optimal anomalous bodies imaged by the CRA inversion are similar to those found for the LRA. However, the Pozzuoli anomaly vanishes and the anomalous bodies are imaged at a slightly shallower depth ( Figure 6). The long axis of the Astroni anomalous body is shifted slightly toward the north with respect to the LRA body. Figure 7 shows the surface traces of the −600, −450, and −300 kg m −3 density isosurfaces. The inversions of both the LRA (Figure 7a) and the CRA data (Figure 7b) consistently image the Astroni and the Solfatara/Pisciarelli anomalies in the same locations. The Solfatara/Pisciarelli anomaly covers the SE edge of Solfatara crater and extends to Pisciarelli in both cases. The Astroni anomaly is centered SW of Astroni crater and covers its SW wall. The imaging and colocation of the Astroni and Solfatara/Pisciarelli anomalies in both models is an indication of the robustness of the inversion, while the veracity of the Pozzuoli anomaly remains uncertain.

Imaging of Distinct Reservoirs: Subsurface Controls on Fluid Distribution
We present the first gravimetric image of the hydrothermal system at CFc. Inversions of two differently detrended data sets (LRA and CRA) provide robust and reproducible results and image two low-density reservoirs beneath Astroni and Solfatara volcanoes, which we interpret as shallow-seated, fluid-rich hydrothermal reservoirs. An anomaly beneath Pozzuoli is only imaged by one of the models, which may be attributed to the lack of offshore gravity data in this survey, potentially preventing us to properly account for the effect of the caldera fill on the data at Pozzuoli. Both model results for Pozzuoli are plausible, and alternative evidence is required to support the existence or absence of a low-density reservoir beneath Pozzuoli (see below).
The optimal modeled negative density anomalies indicate an ∼20% reduction in subsurface density from background values. This can be explained by a porous and fractured caldera fill containing hydrothermal fluids. Borehole data indicate drained host rock (dominantly volcanic tuff) densities ρ r between 1,600 and 2,200 kg m −3 in the top 1 km beneath the caldera containing between 5 and 40 vol% inherent void space (Piochi et al., 2014). To explain the modeled negative density contrasts, a reduction in background bulk host rock density is required. In the hydrothermally active areas imaged in this study this can, for example, be achieved by the generation of additional (secondary) void space by fracturing and/or hydrothermal dissolution (Scenario 1) or replacing the liquid phase in undrained porous host rock by a vapor phase (Scenario 2). In the former case, the background bulk densities will be those reported above, while, in the latter case undrained bulk densities of the caldera fill are in range of 1,650 to 2,500 kg m −3 for given porosities.
We first explore Scenario 1 of bulk density reduction from an average background host rock density ρ r of 1,900 kg m −3 . Assuming that the reduction in density is primarily driven by the creation of new void space φ that is fully connected and can host hydrothermal fluids in either vapor (density ρ v = 1.5 kg m −3 ) and/or liquid (density ρ l = 1,000 kg m −3 ) form, the optimal anomalous density contrast Δρ of the reservoirs of ∼ −400 ± 25 kg m −3 can be explained by a multiphase mixture of caldera fill containing an additional interconnected void volume fraction of between 0.2 and 0.3 that contains a vapor volume fraction ψ v of between 0.38 and 1 and a liquid volume fraction ψ l fraction 0 and 0.62. The parameter space of conceivable fractions of solids and voids (filled with vapor and/or liquid) that fit the optimal model for this scenario is shown in Figure 8 and can be reproduced by Given the borehole rock porosity and density ranges, the second end-member scenario (vapor replaces liquid in undrained caldera fill) is only feasible in rocks containing inherent connected void fractions of 0.4 or more in order to explain the optimal density contrast and is hence if at all only relevant for the top few hundred meters beneath the surface (Piochi et al., 2014). While the most plausible interpretation of the optimal models is a combination of processes associated with both explored scenarios, to explain the modeled density contrast at depths >250-m Scenario 1 must be dominant and from our gravity contrast model alone we favor the creation of additional void space. However, additional constraints are available to help explore the scenarios further.
The permeability structure of the central part of the caldera is key to understanding the distribution of fluids at the time of the survey. Hot low-density fluids will rise from their source until they attain neutral buoyancy, reach the surface, or encounter a barrier to flow, that is, a zone of reduced permeability. The results suggest the presence of fluid-rich bodies trapped beneath the surface of the CFc at Pozzuoli and Astroni volcano, while one body discharges freely at Solfatara volcano and Pisciarelli. This implies the presence of an impermeable seal preventing access of fluids to the surface at Pozzuoli and Astroni volcano. Geochemical and electric data indicate the presence of a two-phase hydrothermal plumbing system at the CFc with a gas-dominated regime residing at shallow (few tens to hundreds of meters) depth beneath the center of the caldera (Chiodini et al., 2011;Gresse et al., 2017).
Permeabilities measured in situ in boreholes at the CFc vary over 4 orders of magnitude (<10 −18 to >10 −14 ) (Piochi et al., 2014). Total horizontal gradiometry of Solfatara volcano ( Figure 5) shows a correlation of the geometry of low-density bodies with the main fault and fracture systems mapped in the field . The combination of in situ rock permeability and fracture/fault permeability may explain the distribution of surface expressions of hydrothermal activity in the caldera and their spatiotemporal evolution. Alunitic alteration at the Solfatara and Pisciarelli hydrothermal fields increases rock porosity and permeability and reduces density . Critically stressed faults can be hydraulically conductive (Jasim et al., 2015), while mineralization can seal previously connected pathways within or around a fault (Sibson, 1994). Faults can thus be both permeable pathways and impermeable inhibitors for fluid flow, depending on the stress regime and degree of alteration. Fluids themselves can modulate permeability via thermally induced hydraulic fracturing (Cusano et al., 2008;Knapp & Knight, 1977;Saccorotti et al., 2007). Permeability in hydrothermal systems is not static and changes constantly due to fracturing and cementation from hydrothermal precipitation, meteroic water invasion, or tectonic stresses (Rowland & Sibson, 2004). Faults at the CFc present fluid pathways on timescales of 1-10 years and 10 2 −10 3 years (Vilardo et al., 2010). It is hence conceivable that geophysical surveys conducted over the past few decades at the CFc provide different snapshots of a constantly evolving hydrothermal system. We therefore compare and contrast the Pozzuoli, Astroni, and Solfatara/Pisciarelli anomalies imaged by our study with published results from other investigations (see also Table 1).

The Pozzuoli Anomaly
Several studies found evidence for a shallow-seated hydrothermal reservoir beneath Pozzuoli, but its nature is contested: • Vanorio et al. (2005) use seismic velocity tomography to delineate a zone of high V p /V s ratios (2.3) beneath Pozzuoli centered at approximately 0.8-km depth and 0.8 km in radius. Below this, they image a low V p /V s (1.4) anomaly at 4-km depth. They interpret these features as a brine caused by steam condensation and a gas enriched formation, respectively. Journal of Geophysical Research: Solid Earth • Chiarabba and Moretti (2006) find a high V p /V s region below Pozzuoli at 0-to 2-km depth, overlying a low V p /V s anomaly at 3-km depth, which they interpret as steam condensation and gas accumulation, respectively. • Seismic attenuation tomography by De Siena et al. (2010) in tandem with V p /V s ratios image an anomaly 0-2 km below Pozzuoli, with a ∼1-km radius. The nature of this anomaly (liquid or gas dominated), however, remains ambiguous in the study. • Chiodini et al. (2015) and Caliro et al. (2007) use geochemical models to predict the vaporization of fluids at 2-km depth, which then rise to the surface. • A stacked (gas-rich pockets beneath liquid-dominated systems) arrangement of fluids is predicted by fault-controlled fluid flow modeling for CFc (Jasim et al., 2015).
In summary, the available evidence is inconclusive regarding the nature of a shallow-seated reservoir beneath Pozzuoli with indications for either vapor-or liquid-dominated regimes. One aspect that needs consideration is the potential for temporal change in the subsurface phase relationships. Chiarabba and Moretti (2006) and M. Battaglia et al. (2006) use combined data from 1984 and from 2001, while Vanorio et al. (2005) and De Siena et al. (2010) used data from 1984, only. Not only do the combined data sets risk masking of temporal signals, but the system may change over the course of 15 years. Temporal changes in elasticity of the upper crust have been suggested for CFc (Di Luccio et al., 2015). Cycles of sealing, fracturing, and resealing are implied on timescales of decades to centuries, for example, at Yellowstone caldera where drill cores plugged almost completely after 25 years (Dobson et al., 2003;Ingebritsen & Sorey, 1988). Our modeling results are consistent with either a liquid-or vapor-dominated system with the caveat that a liquid-dominated regime requires a significantly higher connected porosity compared to a vapor-dominated system to explain the gravity data (Figure 8).

10.1029/2019JB019231
Journal of Geophysical Research: Solid Earth

The Astroni Anomaly
The Astroni anomaly is unexpected as there are no records of fumaroles or other geothermal manifestations in the area. The anomaly is located at the convergence of two crater walls, and one might expect high-density material here compared to adjacent low-density crater fill. However, there are several geophysical anomalies associated with Agnano caldera (Isaia et al., 2009), in which Astroni volcano is nested: • Astroni volcano is seismically active Saccorotti et al., 2007). de Lorenzo et al. (2001) andDe Siena et al. (2010) image a seismic anomaly beneath Agnano at 0-to 3-km depth, which they relate to a high-temperature aquifer. • Capuano et al. (2013) identify an E-W elongated low-gravity anomaly at Astroni volcano at 1-to 2-km depth, which they interpret as a low-density gas-rich reservoir with secondary mineral precipitation as a mechanism for preventing surface expression. They suggest that hot gases condense near the surface and generate a water table which in the Agnano plain forms lake Agnano. • Troiano et al. (2019) identify a highly resistive (vapor-rich) anomaly corresponding to the southeastern edge of the Astroni crater, while the Agnano plain is characterized by a mainly conductive anomaly. • Water from the Agnano well has high P CO2 values, less negative carbon isotope signatures than meteoric water and high HCO −2 3 concentrations, indicating a contribution of heat and hydrothermal fluids from a magmatic system (Venturi et al., 2017).
In light of these geochemical and geophysical observations we suggest that the gravity anomaly of Astroni volcano is formed by a liquid-dominated highly fractured geothermal reservoir trapped beneath a less permeable cap rock. Localized seismicity may be due to the movement of fluids along faults, and it is plausible that there are pathways linking Solfatara/Pisciarelli and Astroni volcano via a network of mainly NW-SE and NE-SW to NNE-SSW trending faults.

The Solfatara/Pisciarelli Anomaly
The key features of the Bouguer anomalies at Solfatara/Pisciarelli are (i) the gravity highs on the SW and NE crater walls of Solfatara(>0.5 mGal), (ii) the gravity low on the southeastern side of Solfatara's crater floor (> −1 mGal) and adjacent Pisciarelli, and (iii) the moderate gravity low near La Fangaia (∼ −0.3 mGal) (Figure 4). These compare to the findings from an earlier survey of Solfatara (Oliveri del Castillo et al., 1968): (i) an elongate ∼1 mGal gravity high on the NE crater wall, (ii) two connected gravity lows reaching from the southern edge of La Fangaia (∼ −0.3 mGal) to Bocca Grande and Bocca Nuova (∼ −0.4 mGal), and (iii) a small scale ∼ −0.3 mGal gravity low in the western side of the crater (see Figure 9). Bruno et al. (2007) demonstrated the spatial correlation of the first three anomalies of the 1968 survey with areas of maximum seismic noise, areas of high CO 2 degassing, and elevated temperatures. The horizontal derivative of the 1968 gravity data highlights the role of faults in concentrating these density anomalies (Bruno et al., 2007).
We divide the anomalies of our survey into three classes (relative to the average gravity of the crater floor), (i) low, (ii) moderately low, and (iii) high to discuss their relationship with other recent geophysical observations at Solfatara/Pisciarelli (see Figure 9 and Table 1).

Low Gravity Class
The low gravity anomaly extends from the eastern side of the crater floor and crosses the crater wall toward Pisciarelli. It coincides with some of the main fumaroles of the geothermal field, which at Solfatara are clustered on its eastern side, where the rocks are intensely fractured . We explain the correlation of the gravity low and areas of intense gas emissions in two ways.
First, upwelling vapor replaces water in fractures and pore spaces in the rock causing a local decrease in density.
Second, low-density material is expected from intense hydrothermal alteration near the fumaroles where connected porosities can reach values of up to 61 vol%  with little variation in hydraulic parameters expected to depths of ∼500 m . Proximal areas around the fumaroles should therefore experience a substantial reduction in density due to the coupled effects of fluid flux and alteration.
At the time of the survey the most vigorous fumarolic activity was located at Bocca Nuova, Bocca Grande, and Pisciarelli, with no fumarolic activity noted on top of the crater wall. These fumaroles are also the 10.1029/2019JB019231 Journal of Geophysical Research: Solid Earth hottest of those in the Solfatara/Piscarelli area (Chiodini et al., 2001). Isaia et al. (2015) have suggested a link between Bocca Nuova, Bocca Grande, and Pisciarelli via faulting and fracturing through the crater wall. This region, therefore, may be the main pathway for fluids to ascend from depth. Other fumaroles may be fed less voluminously and/or by narrower, subordinate fracture networks, which are below the spatial resolution capability of the gravity survey. The strongest first and second horizontal gravity gradients (see Figures 5a and 5b) are across the crater floor and along the NE crater indicating that the anomalies are strongly influenced by faults and fractures.
Di  found a high resistivity zone close to Bocca Nuova and Bocca Grande but little correspondence between resistivity and the other fumaroles. The authors attribute the offset between the resistive body and the 1968 gravity low (see Figure 9) to fluid migration over the time between surveys, but our low gravity anomaly encompasses both the high resistivity body and the 1968 gravity low. Similar to the 1968 survey, our gravity low increases eastward from the main fumaroles ( Figure 9) toward Pisciarelli.
Solfatara undergoes spatiotemporal variations in ground deformation (D'Auria et al., 2012) and has high levels of seismic noise in an arcuate band from the south to the northeast of the crater floor (Bruno et al., 2007). Saccorotti et al. (2007) show maximum likelihood locations of long-period (LP) earthquakes clustered at 500-m depth beneath the SE rim near Bocca Grande. They interpret the LP signals as due to vibrations of fractures in a buried cavity filled by a water-steam mixture. An interpretation of ultrahigh-resolution seismic imaging of the center of Solfatara divides the first 30 m into a shallow zone of aerated tephra underlain by a liquid saturated layer that is deepening in the direction of La Fangaia and a deeper gas accumulation sloping upward toward the eastern side of the crater (De Landro et al., 2017). Seismic attenuation, a shear wave velocity anomaly, and low V p /V s ratios (Chiarabba & Moretti, 2006;Chiodini et al., 2015;De Siena et al., 2010) indicate a shallow gas reservoir around sea level beneath Solfatara. Byrdina et al. (2014) show high resistivity anomalies in both the crater walls and beneath the crater floor. Bocca Grande is directly above a narrow zone of moderate resistivity and surrounded by a zone of high temperature, high CO 2 , and low self-potential. To the east of Bocca Grande is a high resistivity anomaly, and to  (2015) and gravity (1968) anomalies of the Solfatara crater with the LRA data from this study. Dotted lines denote low-amplitude anomalies, dashed lines denote moderately low amplitude anomalies, and solid lines denote high-amplitude anomalies. The crater outline, Pisciarelli, La Fangaia, fumaroles, and lava domes are shown. The color scheme to identiy the LRA gravity anomalies is given in the legend. the west is a low resistivity zone extending toward La Fangaia. Magnetotelluric (MT) and electrical resistivity (ER) data (Troiano et al., 2014(Troiano et al., , 2019) depict a moderately resistive anomaly below the eastern crater wall and a high resistivity anomaly beneath the main fumaroles and Pisciarelli to 2.25 km depth with a radius of ∼0.15 km. High-resolution ER tomography images a gas-dominated reservoir at 60 m depth beneath Solfatara's crater floor that feeds Bocca Grande (Gresse et al., 2017).
In summary, we propose that the main gravity low is caused by a shallow-seated (<1,000 m depth below sea level) accumulation of a two-phase fluid within highly fractured and porous host rocks (Figure 8). It is plausible that the eastern crater wall is composed of highly altered rocks with elevated porosity compared to the rest of the rim, indicating relict and/or current fluid pathways. The imaged feeder system appears to encompass the most dominant pathway for ascending fluids in the central sector of the caldera through a combination of Solfatara's maar-diatreme structure (Troiano et al., 2019) and its intersection with the dominant fault systems (NW-SE and NE-SW to NNE-SSW) of the caldera.

Moderate Low Gravity Class
A moderate gravity low is located near La Fangaia (Figure 9) and matches the extent of the 1968 gravity low. The exact location of the shallow-most La Fangaia feeder system is hard to establish as it appears to change with time. Dried up pits were present during the survey, which must previously have been mud pools. We therefore use the fenced-off area to delineate La Fangaia (brown shape in Figure 9), although its most active portion was the southern part of the area at the time of the survey. Gentle first and second horizontal gravity gradients bound La Fangaia and show the same general direction as faults mapped in the field   (Figures 5a and 5b). The moderate gravity low of La Fangaia is also characterized by high CO 2 flux, low resistivity (with a long lobe extending westward beneath the surface), low self-potential, high temperature, elevated seismic noise, and earthquake clustering (Byrdina et al., 2014) (Table 1). These authors report a positive correlation between CO 2 flux and ground temperature, which are both anticorrelated with selfpotential. The lobate geometry of low resistivity is matched by the moderate gravity anomaly beyond the boundary of La Fangaia. Seismic noise is high near the gravity anomaly and has been correlated with anomalous CO 2 degassing (Bruno et al., 2007). Byrdina et al. (2014) interpret the La Fangaia geophysical anomalies by a liquid saturated plume with both a downwelling condensing liquid water and an upwelling vapor and CO 2 mixture. The water table is locally raised at Solfatara, outcropping at La Fangaia (97 m above sea level) and only 7 m below the surface at the OAK well nearby (see Bruno et al., 2017). A two-phase (gas and liquid) flow regime feeding La Fangaia is also proposed by numerical modeling (Rinaldi et al., 2011). This suggests that background densities for the crater floor are influenced by the presence of liquid water.
The moderate gravity low is thus likely formed by a CO 2 -bearing hot aquifer contained within altered and high-porosity crater-fill.

High Gravity Class
While the depicted gravity highs are constrained only by a low number of benchmarks, the anomalies coincide with the location of the Solfatara cryptodome (northeastern high), and the Mount Olibano lava dome (southern high) ) (see Figure 9). The northeastern high matches a gravity high detected by the 1968 gravity survey. Although poor accessibility prevented us from obtaining more measurements on the southern and northeastern rims of the crater, the transition from low-density crater fill to the high-density Solfatara cryptodome is well marked by a strong first horizontal gradient of the BA. We therefore interpret the gravity highs as remnant domes forming part of the crater rim.

Conclusions
The combined use of high-precision gravity and GPS measurements, high-quality digital elevation models (DEMs), and 3-D data inversion has shed light on the shallow-seated hydrothermal system at Campi Flegrei. The results complement a wealth of existing geophysical and geochemical data for the CFc and lend support to a model of a complex subsurface hydrothermal structure feeding the active fumarolic areas of the central sector of the caldera. We were able to delineate the shallow-seated two-phase hydrothermal plumbing system beneath Solfatara and Pisciarelli and identified a hydrothermal reservoir beneath the Astroni crater. The main gravity anomalies of Solfatara volcano detected in the new survey broadly match those identified by a previous gravity survey conducted in 1968. However, we show that some smaller anomalies may have evolved in size and location over time. This may indicate that within the resolution capabilities of the Bouguer gravity surveys, Solfatara's main hydrothermal feeder system remained broadly unchanged over the past 50 years with the exception of an enlargement toward Pisciarelli, which over the past 15 years has seen a strong increase of hydrothermal activity. Whether or not a separate hydrothermal system resides beneath Pozzuoli cannot be unambiguously answered by our findings, but there are indications for a shallow-seated low-density hydrothermal reservoir during the time of our survey. We encourage additional geophysical and geochemical studies particularly at Astroni volcano and Pozzuoli to test our model results.