Four Years of Continuous Seafloor Displacement Measurements in the Campi Flegrei Caldera

We present 4 years of continuous seafloor deformation measurements carried out in the Campi Flegrei caldera (Southern Italy), one of the most hazardous and populated volcanic areas in the world. The seafloor sector of the caldera has been monitored since early 2016 by the MEDUSA marine research infrastructure, consisting of four instrumented buoys installed where sea depth is less than 100 m. Each MEDUSA buoy is equipped with a cabled, seafloor module with geophysical and oceanographic sensors and a subaerial GPS station providing seafloor deformation and other environmental measures. Since April 2016, the GPS vertical displacements at the four buoys show a continuous uplift of the seafloor with cumulative measured uplift ranging between 8 and 20 cm. Despite the data being affected by environmental noise associated with sea and meteorological conditions, the horizontal GPS displacements on the buoys show a trend coherent with a radial deformation pattern. We use jointly the GPS horizontal and vertical velocities of seafloor and on-land deformations for modeling the volcanic source, finding that a spherical source fits best the GPS data. The geodetic data produced by MEDUSA has now been integrated with the data flow of other monitoring networks deployed on land at Campi Flegrei.


INTRODUCTION
Seafloor geodesy is now an established branch of geodesy for precise underwater measurements of position and displacements and is used for monitoring submerged volcanic areas, tectonic deformation and subsidence related to offshore oil extraction. The main analysis methodologies used are based on acoustic waves propagation and pressure measurements.
Acoustic methods use travel time measurements of acoustic wave packets between two reference locations together with estimation of the seawater sound speed. Two main acoustic techniques have been developed: 1) a GPS-Acoustic combination technique (GPS/A) where a kinematic GPS on a platform at the sea surface is combined with acoustic wave ranging to determine the position of a network of transponder on the seafloor (Bürgman and Chadwell, 2014) This technique has been used extensively to investigate crustal deformation at oceanic spreading centers and subduction thrust faults; usually the measurements are performed in campaign-style with research vessels or light inexpensive floating GPS stations using a glider. 2) horizontal distance measurements between two or more transponders; this is applied to relative deformation measurements across geological features (Petersen et al., 2019).
The main limitation of the acoustic techniques arises from the variation over time of properties of the marine environment, with the consequent difficulty estimating the path length starting from the measurement of the propagation time.
Another common technique in marine geodesy that is suitable for monitoring vertical ground displacement is based on the variation of hydrostatic pressure at the sea bottom measured by a bottom pressure recorder (BPR). However, BPR measurement suffers from greater uncertainties than GPS. To achieve the centimeter accuracy in the Up component need for seafloor deformation monitoring, BPR time series must be integrated with other complementary measurements such as sea level observations at local tide gauges, atmospheric pressure measurements and measures related to seawater column density (i.e. salinity and temperature). Moreover, an erratic instrumental drift of BPRs must be taken into account (Chierici et al., 2016).
Geodetic measurements on the seafloor are more frequently performed at large depths, typically greater than 1,000 m, than in shallow water. This is because the large variability over time of the sound speed strongly degrades the quality of acoustic measures. The sound speed in sea water depends on the temperature and salinity, which, in turn, are strongly affected by solar radiation, by seasonal cycles, and by water mixing due to sea currents. This variability is particularly pronounced in the shallow layers of coastal areas, where additionally, the presence of rivers or waste waters introduce further variability in water properties. Similarly, the variability of temperature and salinity over the time produces sea density changes detected as pressure variation by the BPR sensors. Consequently, long term pressure measurements are contaminated by a high level of noise that can be removed only with precise knowledge of the sea water density profile over the time (Chierici et al., 2016).
A new methodology for seafloor geodesy measurements suitable for shallow water environments, less than 100 m depth, was recently introduced in Southern Italy. A buoy, named CUMAS (Cabled Underwater Multidisciplinary Acquisition System, Iannaccone et al., 2009;Iannaccone et al., 2010) was deployed in November 2006, at a sea depth of 97 m in the Gulf of Pozzuoli, the marine sector of the Campi Flegrei volcanic area, for an experiment of acoustic data transmission. In 2012, the CUMAS buoy was equipped with a GPS station.
The CUMAS buoy consists of an underwater floating body attached to a concrete anchor ballast on the seafloor. The ballast holds the floating body fixed a few meters below the sea level so it is not affected by sea level variations (Iannaccone et al., 2009). A pole supporting a subaerial, top turret is inserted into the body, allowing a GPS station installed on top of the pole to be rigidly coupled to the seafloor and allowing direct recording of seafloor

MEDUSA Infrastructure Outline
The MEDUSA infrastructure consists of four spar buoys installed in the Gulf of Pozzuoli at 1.1 to 2.4 km from the coast and water depths less than 100 m (colored dots in Figure 1).
Two buoys (CFBA and CFBC), deployed where seafloor depth is about 40 m, have a long pole connecting the float and subaerial part to the ballast on the sea bottom. For these two buoys the distance between the GPS antenna and the ballast is about 47.5 m. The other two buoys (CUMAS prototype and CFBB), operating where seafloor depth is about 96 and 76 m respectively, use a steel cable to fasten the pole to the ballast ( Figure 2A). CFBB has a pole the same length as for CFBA and CFBC joined to a steel cable of 33 m length, while the CUMAS buoy has a pole of 26 m long and a steel cable of 83 m. The floats of CUMAS and CFBB provide a buoyancy force acting on the ballast and keeping the steel cable under tension. Functionally, each buoy behaves as a system rigidly anchored to the seabed so the top of the buoy where the GPS antenna is installed follows directly vertical movement of the seabed. Figure 2B, shows in a simplified view the length of the buoys compared to a reference monument and also the proportion of the pole length for each buoy type.
Each buoy is connected by a power supply and data transmission cable to a seafloor module equipped with geophysical and oceanographic sensors placed a few meters from the buoy ballast ( Figure 2A). Each module includes: a threecomponent broadband seismic sensor, two low-frequency omnidirectional hydrophones, a three-component MEMS accelerometer and an absolute pressure recorder. A titanium vessel contains acquisition systems, status and control sensors, and electronics. Power supply, data transmission and clock signals to the module equipment are provided through an armored, electro-mechanical cable, which also serves for module handling and installation/recovery. On the top of the buoy, a subaerial turret hosts the power supply system (solar panels and buffer batteries), and data transmission and control electronics. Real-time data transmission to a land acquisition center is provided by 5 GHz Wi-Fi. Each buoy is equipped with a Leica GR10 GPS receiver and AR20 antenna. GPS data are acquired with a sampling interval of 30 s and downloaded daily by the INGV-Monitoring Center in Naples for storage and processing. A digital compass measuring the attitude of the buoy is mounted on the turret. One buoy is also equipped with a meteorological station (temperature, humidity, atmospheric pressure, wind direction and speed).
MEDUSA has allowed extension into the Gulf of Pozzuoli of the Campi Flegrei geophysical monitoring system, which was exclusively land-based before 2016, improving the observation and assessment of seafloor deformation in the area and observation and location of earthquakes. MEDUSA is fully integrated in the Campi Flegrei surveillance system: all the data acquired are transmitted to a unique, land acquisition center where they are analyzed and stored together with the data provided by the other land monitoring networks.

cGPS Data Processing and Time Series Analysis
Continuous Global Positioning System (cGPS) data are available since April 2016 for the CFBA, CFBB and CFBC stations and since July 2016 for the CFSB station. The cGPS data has been processed in kinematic mode with the open source program package RTKLIB ver. 2.4.2 (http://www.rtklib. com) to obtain positions of the buoys every 30 s. The reference station for cGPS data processing is LICO (Figure 1), located at Note the different scale for the signal on the horizontal components (North and East) with respect to the vertical component (Up). This difference shows the influence of weather and sea conditions and marine water circulation on the horizontal components ( Figures 3B-6B) which produces large fluctuations of the GPS signal, about 10% of which are above about 30 cm. Consequently, we consider that under bad weather conditions the measurements of the position of the buoys is unreliable and we remove all data with fluctuations above ±30 cm in the North and East components. As shown in Figures  3C-5C, about 10-15% of horizontal positions at CFBA, CFBB and CFBC are outside ±30 cm of the corresponding origin (mean position in the considered period) and were discarded to reduce the noise in the North and East components for these stations. Moreover, Figures 3-6, shows that the amplitudes of the GPS horizontal displacements for the three newer buoys (CFBA, CFBB and CFBC stations) are about one third of those of the older CUMAS buoy (CFSB station). This difference is due to the improved design of the new buoys compared to CUMAS and also to the greater length of the CUMAS steel cable (see Table 1).
The CFSB station is characterized by a very high noise level ( Figure 6) and so is not been included in the horizontal component analysis. Its vertical component has been corrected  for vertical displacement induced by the horizontal movement of the buoy allowed by the long steel cable, as reported in De Martino et al. (2014a). This vertical correction is negligible for CFBA, CFBB and CFBC stations, due to the mechanical system adopted in the deployment of the three new buoys (Iannaccone et al., 2018).
In the Campi Flegrei caldera the deformation rates are relatively low (a few millimeters per month) and the trend of the displacement time-series can be best visualized with averaging to produce less frequent but more accurate positions (Larson et al., 2010). A weekly average is routinely applied (De Martino et al., 2014b) to data produced by the permanent Neapolitan Volcanic Continuous GPS (NeVoCGPS) network deployed on land (black circles in Figure 1. Here, for the MEDUSA data, we generate the weekly median time series of the cleaned kinematic solutions to which we applied an additional cleaning algorithm based on the interquartile range (IQR) statistic for outlier detection in the GPS position time series (Bock et al., 2000). An outlier is defined as having a residual with respect to the weekly median greater than 2 times the IQR. With this criterion, the algorithm removed about 5% of the kinematic data.
The weekly filtered GPS time series are shown in Figure 7. The uncertainties of the weekly positions (error bars) shown in Figure 7 represent the interquartile range (IQR) of each weekly median solution.
In the considered period, the weekly time series show that the four buoys are subject to uplift and to radial horizontal movements with respect to the center of the caldera, as previously found by Iannaccone et al., 2018. Table 2 reports the correlation between the cGPS horizontal components of the couples of the buoys CFBA, CFBB and CFBC, in the period between 1 March 2018 and 31 May 2020.
We interpret the high correlation between the horizontal components as due to the effect of the wind and marine currents acting coherently on the three buoys CFBA, CFBB and CFBC. The buoy CFSB is located at greater depth and is subject to larger horizontal fluctuations. For this reason we have not performed the correlation analysis for CFSB. Figure 8 shows the vertical displacement observed at the RITE cGPS station located in the town of Pozzuoli, the area of maximum vertical deformation of the caldera, and at the four MEDUSA cGPS stations on the buoys.
The consistency of the vertical displacement patterns between the stations is remarkable with a maximum uplift of the seafloor of about 20 centimeters at the CFBA station in the last four years and a decay of amplitude moving from the center to the border of the caldera.

RESULTS
Recently, Iannaccone et al. (2018) demonstrated that a single Mogi model point source can explain both the onshore and seafloor deformation from cGPS data. They compared the expected vertical displacement in the marine sector to that predicted by a Mogi model computed using only horizontal and vertical, on-land, cGPS measurements in the period between April 2016 and July 2017. This period, however, was characterized by small vertical movements.
From July 2017 to May 2020, the Campi Flegrei caldera has been characterized by larger vertical motion with a constant, quasi-linear trend over time (http://www.ov.ingv.it/ov/it/ bollettini/275.html). We use jointly the cGPS horizontal and vertical velocities of seafloor and onshore deformation for the modeling of the volcanic source.
The horizontal and vertical velocities, uncertainties, and noise properties were estimated from the weekly position time series of the 3 cGPS stations on the buoys (CFBA, CFBB and CFBC in Figure 7) and of the 21 cGPS stations deployed on land (De Martino et al., 2014b;Iannaccone et al., 2018), using the software package HECTOR (Bos et al., 2013). This software, based on the maximum likelihood estimation algorithm (MLE), estimates the linear trend in time series with temporally correlated noise. A flicker noise plus white noise model is assumed to take into account the existing noise signals in the time series and to compute more realistic uncertainties for the velocities (Mao et al., 1999;Williams 2008). An annual seasonal signal, partially caused by atmospheric loading, thermal expansion or multipath (Klos et al., 2018 and references therein), is estimated and removed.
The results of the velocity estimations are presented in Table 3 and Figure 9 (black arrows).
Next, we jointly inverted the horizontal and vertical deformation velocities measured at all the 24 cGPS stations (Figure 1) since July 2017 to investigate the caldera deformation. We used the dMODELS software package (Battaglia et al., 2013a;Battaglia et al., 2013b), which provides MATLAB functions and scripts to invert GPS data for the typical source shapes used in the volcanic environment (see eg: Gottsmann et al., 2006), such as spherical (McTigue, 1987), spheroidal (Yang et al., 1988) and sill (Fialko et al., 2001), embedded in a homogenous, isotropic, elastic half-space.
The results for different volcanic sources obtained by cGPS data inversion are listed in Table 4.
The statistics shows that a spherical source (McTigue, 1987) fits best the GPS deformation velocities. The sill-like source has larger values of Np (number of source parameters) and χ v 2 (chi-squared per degrees of freedom) and the spheroid has a smaller value of χ v 2 than the spherical source but also has 3 additional geometrical parameters (Battaglia et al., 2013b). Figure 9 shows the observed and calculated horizontal and vertical cGPS velocity field for the spherical source. The position of the best fit source to the cGPS data from July 2017 to May 2020 (red cross in Figure 9) is very similar to source position inferred by Iannaccone et al., 2018 in the time span April 2016 to July 2017.

DISCUSSION AND CONCLUSION
After a test period and four years of operation and improvement of the quality control and analysis  Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 615178 8 The deformation field observed in the Campi Flegrei area is generally axi-symmetric with the maximum of vertical displacement measured at the RITE GPS station (Figure 1), at the center of the caldera (De Martino et al., 2014b). The vertical deformation decreases radially from the center and vanishes on the caldera's border, about 6 km away. In this work we confirm a similar pattern is obtained with the joint use of marine and land GPS data. We demonstrated that a spherical source (McTigue, 1987) is in good agreement with both the vertical and the horizontal cGPS measurements.
This simple symmetric pattern would be perturbed if there were a lateral intrusion, such as a dyke, in the marine sector in the future. In this case, MEDUSA measurements would be essential for an early detection of an eventual magma rise in the marine sector of the caldera.
Despite the high noise level of the horizontal GPS components of MEDUSA, our results showing coherence with the corresponding onshore data indicate that useful signal can be obtained. However, the methodology we used for the recovery of the horizontal deformation components can be applied only after an observed monotonic trend of at least 1-year. This limitation implies that this method is not practical for real-time application, unless it is possible to achieve a higher accuracy and further reduce noise due to marine weather effects on the horizontal data.
Given the overall length of the buoys, from 47.5 to 108 m (see Table 1), it would be very difficult to decrease their movement in the sea. However, noise reduction through data position correction can be performed using heading, pitch and roll measurements produced by a digital compass on the top of each buoy. The MEDUSA infrastructure inspired the University of Southern Florida to build and deploy a similar buoy (SUBGEO-Shallow Underwater Buoy for Geodesy) in Tampa bay, Florida, at a water depth of about 23 m (Xie et al., 2019). Using a digital compass installed on the buoy, Xie et al. (2019) demonstrate the possibility to recover both horizontal and vertical components of the seafloor motion at the base of the buoy with an uncertainty of 1-2 centimeters.
Application of a similar procedure is under study for the MEDUSA buoys which also are equipped with digital compass to measure pitch, roll and heading data.
We have shown that MEDUSA provides continuous GPS measurements in shallow water (less than 100 meter depth) with a comparable precision to the land GPS vertical displacement and competitive with other systems such those using seafloor pressure gauges or GPS-acoustic positioning.  The MEDUSA monitoring infrastructure is now permanently integrated into the multidisciplinary local monitoring network of Campi Flegrei and its data are routinely published within the periodical bulletins produced for the Italian Civil Protection Department and in the warning protocols in one of the highest risky area of the world.  Table 3) and calculated displacements assuming a spherical source (Table 4), respectively. The location of the source is indicated by a red cross.