Radial interpolation of GPS and leveling data of ground deformation in a resurgent caldera: application to Campi Flegrei (Italy)

This study presents a new method, called the Radial Interpolation Method, to interpolate data characterized by an approximately radial pattern around a relatively constrained central zone, such as the ground deformation patterns shown in many active volcanic areas. The method enables the fast production of short-term deformation maps on the base of spatially sparse ground deformation measurements and can provide uncertainty quantification on the interpolated values, fundamental for hazard assessment purposes and deformation source reconstruction. The presented approach is not dependent on a priori assumptions about the geometry, location and physical properties of the source, except for the requirement of a locally radial pattern, i.e., allowing multiple centers of symmetry. We test the new method on a synthetic point source example, and then, we apply the method to selected time intervals of real geodetic data collected at the Campi Flegrei caldera during the last 39 years, including examples of leveling, Geodetic Precise Traversing measurements and Global Positioning System. The maps of horizontal displacement, calculated inland, show maximum values lying along a semicircular annular region with a radius of about 2–3 km in size. This semi-annular area is marked by mesoscale structures such as faults, sand dikes and fractures. The maps of vertical displacement describe a linear relation between the maximum vertical uplift measured and the volume variation. The multiplicative factor in the linear relation is about 0.3 × 106 m3/cm if we estimate the proportion of the ΔV that is captured by the GPS network onland and we use this to estimate the full ΔV. In this case, the 95% confidence interval on K because of linear regression is ± 5%. Finally, we briefly discuss how the new method could be used for the production of short-term vent opening maps on the base of real-time geodetic measurements of the horizontal and vertical displacements.

Several techniques measuring ground deformation in calderas are available (Dzurisin 2006;Poland et al. 2006), including leveling (Dzurisin et al. 2002;Yokoyama 2013;Murase et al. 2014), Geodetic Precise Traversing (GPT; Barbarella et al. 1983;Punmia and Jain 2005), Global Positioning System (GPS, Dixon et al. 1997), tiltmeters, dilatometers, tide gauge (Mori et al. 1986), electronic distance measurements (EDM) and, more recently, interferometric synthetic aperture radar (InSAR, Massonnet and Feigl 1998;Burgmann et al. 2000;Parker et al. 2014;Tizzani et al. 2015;Stramondo et al. 2016). Ground-based techniques are in general sparse in space but almost continuous in time, while InSAR has good spatial resolution but time gaps of hours/ days. In this study, we focus on the most classical techniques such as leveling, GPT and GPS, which usually provide the longest temporal datasets of ground displacement, but that are often significantly sparse in space.
Ground leveling is the oldest technique that has been applied in volcanic areas for monitoring the ground deformation. It compares the elevation in fixed benchmarks with respect to a given datum at a sequence of discrete times. GPT networks involve placing survey stations along a line or path of travel and use the previously surveyed points as a baseline for observing the next point. This method produces more information than traditional leveling, including horizontal displacement measured by a triangulation between the stations, although measured on spatially localized and fixed paths, and still discrete in time. GPS stations using precise satellite positions and clocks have been more commonly utilized in recent years and can provide, unlike of the former methods, time-continuous data of vertical and horizontal ground displacement components, but often refer to a still relatively small number of measurement sites scattered on a large spatial region. A reconstruction of the full ground displacement field by interpolation of sparse GPS measurements was recently undertaken on the volcanic caldera of Santorini (Greece) by Lagios et al. (2013). They performed the comparison between GPS data and space-borne InSAR data without a specific focus on the interpolation method applied. A similar approach was applied by Amoruso et al. (2014a, b), D'Auria et al. (2015) and Trasatti et al. (2015) at the Campi Flegrei caldera. In these studies, a model including the intrusion at shallow depth of magma within a sill was tested using (i) leveling, GPT and InSAR data (Amoruso et al. 2014a, b); (ii) GPS and InSAR data (D'Auria et al. 2015); and (iii) GPS, leveling and InSAR data (Trasatti et al. 2015).
In this work, we present a simple method for interpolating sparse deformation measurements showing an approximately radial pattern, both in intensity and orientation, with respect to a localized central area. The method is called "Radial Interpolation Method" (RIM, in short). It allows the construction of both the vertical and horizontal components of the full displacement field, a task particularly challenging, especially for historical unrest episodes when InSAR data were not available or in places where InSAR is not possible (e.g., poor coherence) (e.g., Newhall and Dzurisin 1988a, b). The new method is applied to the spatially sparse ground deformation measurements recorded in the Campi Flegrei caldera ( Fig. 1) after 1980 by means of leveling and GPT  Gaudio et al. 2010). The plot shows also the analyzed periods and the type of dataset surveys and by GPS continuous stations. We include in the analysis a detailed uncertainty quantification of the new technique. Furthermore, we broadly estimate the volume variation of the ground both during uplift and subsidence episodes according to the RIM displacements maps, and we do a preliminary comparison with those calculated from the processing of InSAR acquisitions (see "Appendix A"). On these volume estimates, we test the effect of various assumptions on the radial extent of the interpolated region and on the offshore portion of the caldera. This volume variation may or may not be equivalent to the volume change at the source of deformation, depending on the source geometry and depth (e.g., Lisowski 2006).
Finally, we discuss how the interpolated data may be used as the input of analytical source models, possibly in conjunction with remote sensing observations, or become a key input for the Bayesian updating of vent opening probability maps, for example, at Campi Flegrei caldera.
This study has the main purpose to present the new interpolation method, and the secondary purpose to show its application to the Campi Flegrei special case. The interpretation from a volcanic processes point of view is beyond our scope.

Geological setting of the analyzed area
We test our approach on the Campi Flegrei caldera, one of the most studied and monitored volcanic areas in the world. We choose this case study because of data availability and of the significant implications in terms of hazard assessment. Campi Flegrei caldera (Fig. 1a) is located on the western coast of the Campania region (southern Italy; Vitale and Ciarcia 2018) in a very urbanized area, including the western part of the city of Napoli. Two caldera-forming eruptions, i.e., the Campanian Ignimbrite (~ 40 ka; Costa et al. 2012;Giaccio et al. 2017) and the Neapolitan Yellow Tuff (~ 15 ka¸ Orsi et al. 1992;Scarpati et al. 1993;Deino et al. 2004), produced two nested calderas, of which the largest is ~ 12 km wide (e.g., Acocella 2008). In the last 15,000 years, this area experienced at least 70 eruptions (Rosi et al. 1983;Di Vito et al. 1999;Orsi et al. 2004;Isaia et al. 2009;Smith et al. 2011;Isaia et al. 2015;Bevilacqua 2016;Bevilacqua et al. 2015;2016). In the last century, Campi Flegrei displayed several uplift and subsidence episodes called bradyseism. Three main phases of uplift occurred in the periods 1950-1952, 1968-1972 and 1980-1984, usually accom-panied by earthquakes at shallow depths and concentrated around the Solfatara crater and beneath the town of Pozzuoli and the Pozzuoli Bay (e.g., Kilburn et al. 2017;Castaldo et al. 2018). This produced a net vertical displacement of ca. 3.8 m from the minimum of 1950, measured in 1985 in the central part of the caldera (Del Gaudio et al. 2010;Amoruso et al. 2014a, b;De Martino et al. 2014a;Troise et al. 2019). Since 2005, after 20 years long phase of subsidence with a total vertical displacement of ca. − 1 m, the Campi Flegrei caldera has been raising again with a slower rate than in the previously occurred main phases of uplift, but slowly accelerating (Chiodini et al. 2017;Bevilacqua et al. 2018Bevilacqua et al. , 2019Patra et al. 2019); http://www.ov.ingv.it/ov/it/campi -flegr ei/ monit oragg io.html), with a total maximum vertical displacement in the central area of ca. 63 cm (November 2019).

Ground deformation data at Campi Flegrei caldera
In summary, the ground displacement measurements recorded at Campi Flegrei in the last 39 years generally show approximately radial patterns both in terms of intensity and orientation, i.e., axial symmetry (Corrado et al. 1977;Berrino et al. 1984;Bianchi et al. 1987;Orsi et al. 1999;De Martino et al. 2014a, b;Amoruso et al. 2014a, b;Iannaccone et al. 2018). The maximum vertical displacement (uplift or subsidence) is located in the town of Pozzuoli, in the central part of the Campi Flegrei caldera (benchmark 25, Fig. 1a).

Leveling and Geodetic Precise Traversing data
Campi Flegrei ground deformation data have been systematically collected from 1970, through the technique of the leveling lines (Corrado et al. 1977;Berrino et al. 1984;Bianchi et al. 1987). This technique can only provide measurements of the vertical displacement and consists in measuring height differences within a network of lines formed by permanent benchmarks, as a function of time. After June 1983, lines (from A to H, Fig. 2a, b) and benchmarks were increased in number (Amoruso et al. 2014a, b and reference therein). In June 1980 and June 1983, the leveling records were accompanied by GPT surveys-i.e., by distance and angular measurements (Barbarella et al. 1983) that allowed, for the first time, to describe also the horizontal component of the ground deformation. Horizontal and vertical component data were later used to infer and model the magmatic source by a number of studies (e.g., Berrino et al. 1984;Bianchi et al. 1987;Troise et al. 2007;Amoruso et al. 2014a, b). In our study, we analyze GPT data (also called EDM, i.e., Electronic Distance Measuring) from Amoruso et al. (2014a,b) and unpublished data collected by the leveling network of the Vesuvius Observatory (INGV Napoli). This includes leveling measurements recorded in the periods 1980-1983 (uplift), 1983-1984 (uplift), 1985-1988 (subsidence) and GPT measurements recorded in 1983GPT measurements recorded in -1984. We summarize these data in Supporting Information SI5 and SI6.

GPS data and InSAR
Campi Flegrei has been continuously monitored from 2000 by a GPS network named "Neapolitan Volcanoes Continuous GPS" (NeVoCGPS). A total of 25 stations are present (20 inland and 4 in the Pozzuoli Gulf) covering an area of about 130 km 2 . More details on data recording and processing are described in De Martino et al. (2014a, b) Martino et al. 2014a). These are called UP1-UP7. Specifically, we analyzed GPS measurements between April 2011 and April 2013. The civil protection authorities raised the alert level of Campi Flegrei to "Yellow" (warning) in 2012, and this level has been maintained for 7 years afterward (in 2019 at the time of writing). These series were recorded by GPS stations during four minor uplift periods named UP4 to UP7 (Fig. 2c, d). We report in Supporting Information SI1 The Campi Flegrei caldera has also been continuously monitored by InSAR technology starting from 1992 (ERS, ENVISAT and Sentinel-1 satellites). By exploiting the deformation measurements taken over the same point from two different acquisition geometries (namely the ascending and descending orbits), it is possible to retrieve the 3D displacement/velocity field by combining ascending and descending scenes (Dalla Via et al. 2012). In "Appendix A," we include two examples of vertical displacement maps based on the processing of InSAR data, focusing on the time window 2004-2007. The related estimates of volume variation of the ground and the vertical uplift observed at RITE station show the same linear relation of those based on our new interpolation method.

The Radial Interpolation Methods
Strain fields associated with magmatic inflation/deflation, and sill intrusion, generally show a radial pattern of the displacement vectors, both in intensity and orientation (e.g., Mogi 1958;Okada 1985;Davis 1986;Fialko et al. 2001;Miyagi et al. 2004;Ukawa et al. 2006;Palano et al. 2008;Chadwick et al. 2011;Hooper et al. 2011;Lagios et al. 2013;Ji et al. 2013;Grapenthin et al. 2013;Galgana et al. 2014). If the data are few and spatially scattered (e.g., GPS measurements-Figs. 2c, 3a), or sparse in localized areas (e.g., leveling data-Figs. 2a, 3b), the common interpolation methods used, such as linear, cubic, natural and spline, can provide inaccurate results because they do not preserve a radial symmetry in the deformation field. The new interpolation method proposed in this study instead exploits the locally radial geometric pattern to better reconstruct the strain field in the region, possibly covering areas far away from the data points.

The algorithm
Our method is based on a multipolar interpolation, i.e., we apply a linear interpolation between pairs of vectors, working in polar coordinates with respect to the intersection of the straight lines defined by their horizontal components (Fig. 3c). In particular, we first provide one-dimensional interpolations along the arcs connecting pairs of points, and then, we fill the rest of the space according to a twodimensional interpolation technique. In the considered examples, we utilize a linear interpolation, but this choice is not restrictive.
Ground displacement vectors are applied vectors, expressed by the coordinates = (x, y) of the point of application of the vector, in addition to the displacement coordinates = (d x , d y , d z ) . What we are interpolating is a field of 3D vectors defined on a 2D domain; thus, the interpolation technique is defined in ℝ 2 × ℝ 3 .
By considering a dataset of applied vectors (i) = (i) , (i) , i = 1, … , N , RIM can be described by the following steps: (A) Graph Construction We define an undirected graph G = (X, E). The set of nodes X ⊆ (i) ∶ i = 1, … , N includes all the acceptable spatial locations, i.e., excluding defective data. The set of edges E includes all the couples E (ij) = (i) , (j) of consistent nodes, i.e., excluding those associated with sub-parallel horizontal we define a pole P (ij) at the intersection between the st r a i g h t l i n e s (i) + k d (i) x , d (i) y |k ∈ ℝ a n d y |h ∈ ℝ . The spatial locations of x (i) and x (j) are expressed in polar coordinates with respect to P (ij) , i.e., x (i) = (r (i) , α (i) ), and x (j) = (r (j) , α (j) ), where r ∈ ℝ + and α ∈ [−π,π]. For simplicity, and without loss of generality, we assume α (i) = 0. A one-dimensional linear interpolation is then accomplished in ℝ + × [− , ] × ℝ 3 , between v (i) and v (j) . In particular, for all λ ∈ [0, 1], we define: In this way, we conjointly interpolate the location and the displacement vector. We set a number of discrete steps n ∈ ℕ , and hence, we define the family V ij including the vectors v (ij)t/n for each t in {1,…,n}. These vectors are placed along the elliptic arc that links x (i) and x (j) , and they are a convex combination of d (i) and d (j) (C) Cartesian Secondary Interpolation We finally apply a two-dimensional interpolation on the virtual displacements V, to fill the rest of the space ℝ 2 × ℝ 3 . Namely, we perform a linear interpolation without extrapolating the values outside the convex hull of V. (We implemented the function interp in the software R.) Alternative methods formulated in Cartesian coordinates have been tested, and their results are not significantly different. Figure 3c shows an example of the multipolar interpolation, considering the vectors V 1 = (x 1 , y 1 , d 1 ) and V 2 = (x 2 , y 2 , d 2 ). The spatial coordinates are (x 1 , y 1 ) and (x 2 , y 2 ), and the azimuth angles are α 1 and α 2 , respectively. The n-interpolated vectors (v (1−2)1 ,…, v (1−2)n ) are placed along the arc of circle centered in P (ij) = (C 1-2 ), at the intersection of the prolongations of V 1 and V 2 . The corresponding azimuth angles (α (1-2)1:n ) are equal to the nth part of the angle between the vectors V 1 and V 2 . Finally, the displacements d 1−2 t/n linearly interpolate those of V 1 and V 2 . The method provides m*ninterpolated vectors, where m is the number of edges in the graph.

Model validation and main sources of uncertainty
In order to quantitatively validate the RIM in an ideal case, we consider a synthetic ground deformation related to a spherical source in an elastic half-space, located at (0, 0, −D) (Mogi model;Mogi 1958;McTigue 1987) where R 2 : = x 2 + y 2 + D 2 , and with the following other properties: (1) depth of the center of the sphere from the surface (D = 3000 m); (2) radius of the sphere (α = 700 m); (3) pressure change in the sphere (ΔP = 20 MPa); (4) Young's modulus (μ = 10 10 ) and Poisson's ratio (ν = 0.25).
In Fig. 4a, we show the corresponding horizontal displacement map. Based on this, we randomly sample a number of uniformly distributed displacement values from this synthetic map, and we build a sequence of interpolation maps from these samples, based on 20, 15 and 10 of them. First, these data are linearly interpolated ( Fig. 4b-d). Then, we use the new RIM method. In Fig. 5, we report a graph of consistent nodes (Fig. 5a, d, g), the virtual displacement vectors (Fig. 5b, e, h) and the interpolated map of RIM (Fig. 5c, f, i). The point densities vary greatly, and the virtual displacement vectors can cluster along the arcs between stations close-by. However, the Cartesian secondary interpolation produces a map on a regular grid and removes these differences. Finally, to quantify the accuracy of the methods, for every interpolated map M we calculate the misfit with respect to the original displacement map A. That is, the error ratio: expressed as percentage. The values of (Err) are displayed in Figs. 4c, d and 5c, f, i, and RIM is associated with significantly lower Err values than in the linear interpolation. In summary, according to the chosen graphs of consistent nodes the new method performs better than a linear interpolation. This result is not depending on these particular graphs: Drawing all the possible arcs or linking the k-closest nodes also performs better than a linear interpolation.
It is important to emphasize that although the RIM can be applied to the vertical component of the ground deformation, it is always necessary to know the azimuths of horizontal components, because those define the elliptic arcs.
We remark that the results of the RIM are not equivalent to use an analytical source model. Indeed, our interpolation procedure is not aimed at accurately inferring the profile of the deformation along the radial direction better than a linear interpolation between the measured data. Instead, the key step of the RIM exploits the assumption of radial symmetry between each pair of measurements, operating along the direction orthogonal to the radii (see Fig. 3c, 5). Significantly, our results are independent of any specific geometry, location and physical properties of the source except for requiring a locally radial pattern, i.e., allowing multiple centers of symmetry. A global center of symmetry is not specified by the RIM, but a central region is approximated by the set of poles of symmetry P (ij) .
The main sources of uncertainty affecting the RIM include: 1. The propagation of the measurement error on the original dataset; 2. The choice of the set of consistent nodes defining the elliptic arcs; 3. The modeling error related to the interpolation itself.
The effect of the first source of uncertainty is calculated with a Monte Carlo simulation. That the measurement error affecting the horizontal displacement can have nonlinear effects on the interpolated maps, because it affects the position of the poles of symmetry P (ij) . We accomplished The second source of uncertainty is also explored with a sensitivity analysis made on the choice of the set of consistent nodes in the Campi Flegrei examples, and the results are attached in Supporting Information SI2. The definition of consistent nodes in the graph is essentially made to exclude sub-parallel horizontal displacements and to avoid elliptic arcs that overlap with clusters of other nodes. However, the comparison of Supporting Information SI2 to the figures in Sect. 5.1 showed that considering all the nodes and arcs is not producing a significant change in the results.
The third source of uncertainty, the modeling error, is quantified in the synthetic example shown in this section (see Figs. 4,5). In a general case, modeling error can only be loosely constrained according to the differences between measurements as a function of their spatial distance. If the radial symmetry breaks, then this error could be significantly large. The modeling error remains small if the assumption of radial symmetry is sufficiently robust, and the following analysis is conditional on that assumption.

Examples based on GPS data
The first class of real examples that we describe are related to four minor uplift episodes that occurred in 2011-2013 in Campi Flegrei caldera. Figure 6 shows the horizontal and vertical displacements measured at 14 GPS stations of the monitoring network over that period. The time series  Fig. 6, as well as the computed centers for each pair of GPS stations considered (poles of symmetry), are related to the average displacements with respect to the measurement error. So, the elliptic arcs related to points characterized by smallest horizontal displacements are very uncertain, i.e., the small-sized virtual displacement vectors have unreliable locations and directions. In Supporting Information SI4, we marked the stations in which the average measurement is lower than its standard error σ or than 2σ. The vertical deformations collected at the most distant GPS stations from the caldera center are not always significant. In particular, 9/15 stations in UP4, 8/15 stations in UP5, 2/15 stations in UP6 and 7/15 stations in UP7 did not collect measurements greater than 2σ. We remark that the RIM does not require to assume a unique center of symmetry since each pair of consistent nodes define their own pole of symmetry. Therefore, in Fig. 6, we show the mean location of these poles of symmetry, excluding from the computation few outliers located further than 5 km from RITE station that have a small horizontal displacement strongly affected by measurement error. Figures 7, 8, 9 and 10 show the interpolated maps of the GPS data for the uplift episodes according to our new method. All the GPS measurements in interval 2011-2013 are located inland, and hence, we are not considering the interpolated data in the Gulf of Pozzuoli, although the corresponding contour lines are also displayed. We also note that the RIM does not extrapolate the deformation outside the convex hull of the elliptic arcs that link the measurement sites. We performed a Monte Carlo simulation of 250 samples varying the observations, considering the measurement error as a source of uncertainty; increasing the number of samples does not significantly affect the result. All the figures include two pairs of plots, expressing displacement values at each location (x, y) on a rectangular grid of 250 m cell size. Plots (a-b) show the horizontal displacement, and plots (c-d) show the vertical displacement. The first plot of each pair shows the median values of displacement M(x, y) = 50th percentile (x, y); the second plot shows the uncertainty range U that we define as: U(x, y) = [95th percentile (x, y) − 5th percentile (x, y)]/2. The interpolated maps of upper and lower uncertainty percentiles are included in Supporting Information SI3.
In summary, the median values of horizontal displacement are in the range [0.1, 1.2] cm for UP4 and UP5, while they are in [0.1, 2.8] cm for UP6 and UP7. The pattern is consistent with the well-known dome-shaped geometry, with maximum of median horizontal displacement always located in a semi-annular region in the central area of the caldera, at 2-3 km from Pozzuoli. During UP4, the peak displacement is located in the center-east area of the caldera (close to ACAE and SOLO stations), for UP5 it is on the center-west area (ARFE and STRZ), and for UP6 and UP7 it progressively moves back to the center-east area (near STRZ and then SOLO stations). Uncertainty range is in [0.1, 0.25] cm for UP4 and UP5, while it is in [0.1, 0.4] cm for UP6 and UP7. These values are not large enough to significantly affect the pattern described above. Most uncertain regions for UP4 and UP5 are located close to the stations that reported least precise measurements. Instead, for UP6 and UP7 the uncertainty peaks are located in regions far from the measurement stations and related to the uncertainty affecting the elliptic arc placement. However, in the four uplifts the median values of the horizontal deformation are larger than their uncertainties in almost the entire region considered, and the deformation is significant with a 95% confidence.
Vertical displacement for UP4 and UP5 is characterized by elliptical isolines with WNW-ESE and NNE-SSW long axis directions, respectively. Instead, for UP6 and UP7 the vertical displacement shows more circular isolines.
The median values are in the range [− 0.3, + 2.3] cm for UP4, in [+ 0.1, + 2.5] cm for UP5, in [+ 0.5, + 6.5] for UP6 and in [− 0.2, + 5.0] for UP7. Maximum values are always close to RITE station, although in UP4 the peak displacement affects also the ACAE station. Uncertainty is in the range [0.3, 1.3] cm for all episodes, with peak values located either near some stations (mostly for U4 and U5) or far from stations (U6 and U7). In Figs. 7, 8, 9 and 10, we marked in red the regions where the vertical displacement was significant with a 95% confidence. In UP6, the interpolated vertical deformation field is significant in almost the whole region. Instead, in UP4, UP5 and UP7, in some of the most distant regions from the caldera center the uncertainty is larger than the median values obtained, and the vertical deformation is not significant with a 95% confidence.
In the future analysis, the same procedure could be adopted to compare either vertical or horizontal velocities over different time periods-indeed the displacements may vary pending on event duration, where the velocities could be more comparable. Moreover, vertical and horizontal displacements could be reconciled and used for volcano monitoring purposes, for example analyzing the 3D displacement field.

Examples based on leveling and Geodetic Precise Traversing data
The second class of real examples that we describe are related to the leveling and GPT data collected during the bradyseismic crisis of 1980-1985 and the subsequent phase of subsidence. We considered three time periods: 1980-1983, 1983-1984 and 1985-1988. The leveling data are based on a much larger number of measurements than the GPS data, but they refer to benchmarks that lie only on localized paths, leaving large areas uncovered. Moreover, leveling data do not provide horizontal deformation. The poles of symmetry are obtained from the processing of the horizontal displacement, and only the GPT data of 1983-1984 enable their calculation. Uncertainty analysis is not performed for these examples due to the lack of information available to us on the measurement error of GPT and leveling campaigns. The displacement values are expressed again on a rectangular grid of 250 m cell size. Figure 11a displays the horizontal displacements measured in 14 GPT stations. The displacement data and the graph of consistent nodes are reported in Supporting Information SI5. In Fig. 11a, we also show the mean location of the poles of symmetry, excluding few outliers located further than 5 km from RITE station. This mean center is located at 426481E, 4518885 N. The RIM was applied on all leveling data by implementing this location as a global center of symmetry. Figure 11b-d displays the vertical displacements measured along the leveling network of the considered periods. The network was significantly extended in 1983, with the setup of additional leveling paths (station labels 52-82). To improve clarity, some displacement vectors are not reported, and some leveling stations are not labeled in the figure. However, the displacement data collected in all the stations, and the graph of consistent nodes, are fully reported in Supporting Information SI6. This dataset also includes the periods of subsidence in 1989-1992 and 1995-2000, not detailed in this section, but considered in the next section for the depiction of the deformation boundary and for the calculation of volume variation of the ground. Figure 12 shows the interpolated maps of the GPT and leveling data according to our new method. In summary, the horizontal displacement measured in 1983-1984 is in the range [12, 34] cm. The GPT stations were all placed in the central part of the caldera, and the data outside of this region are not available. However, the pattern is still consistent with a dome-shaped geometry, with maximum horizontal displacement again located in a semi-annular region in the central area of the caldera, at 2-3 km from the town of Pozzuoli. The peak displacement is located in the center-west area of the caldera (close to the nowadays called ARFE station). This is similar to what observed in UP5, but the displacement scale is an order of magnitude larger.

Deformation boundaries
In this section, we further analyze the interpolated maps. We provide a quantification of the region affected by significant horizontal or vertical deformation, with the purpose to compare it to the past vent locations. In Fig. 13a, we mark the region of maximum horizontal displacement, as obtained from the GPS and GPT data. These lines only have a descriptive purpose and simply connect the points of maximum displacement observed over each radial line with respect to the approximate center of vertical displacements. Similarly, we also draw the external boundaries of vertical displacements, marking significant changes in the slope of the maps. We note that the eruptive vents of the last 5.2 kyr BP (i.e., Epoch III and Monte Nuovo eruption, with spatial uncertainty described in Bevilacqua et al. 2015;Bevilacqua 2016) are all close to these lines or inside the annular region between them. A further analysis of the relation of vent opening and deformation could have important implications for hazard assessment. Indeed, the Fig. 11 Plot a displays the GPT measurements related to 1983-1984. Black arrows show the horizontal displacements. Green dots mark the RIM computed centers for each pair of stations considered (poles of symmetry). A purple cross marks the mean of those centers. Plots bd display the leveling measurements related to 1980-1983 (b), 1983-1984 (c), 1985-1988 (d). Gray arrows show the vertical displacements. Red dots mark the leveling stations, with some label numbers displayed. Blue dots mark the RIM virtual displacement locations, along elliptic arcs. UTM WGS84, Zone 33 N coordinate system (m) 3D deformation field can help to constrain the expected radial distance of eruptive vents from the center of the caldera, and large values of horizontal deformation can be used as a spatial precursor to future vent opening (see Bevilacqua et al. 2015Bevilacqua et al. , 2017Rivalta et al. 2019;Patra et al. 2019).

Estimation of the volume variation of the ground
The interpolated maps of the vertical displacement allow us to evaluate the volume variation of the ground (ΔV) during the major bradyseismic episodes (either of uplift or subsidence). This quantity only estimates the ground deformation without considering any dependence on any specific source geometry, physical properties and location except for requiring a locally radial pattern. However, the volume variation of the ground may or, often, may not correspond to the volume variation of the source, depending on its unknown geometry  Giudicepietro et al. 2016). We remark that, in general, the proportionality between the volume changes at the source and at the surface is dependent upon many factors, such as the elastic parameters and the depth of the source (e.g., Lisowski 2006). In general, we calculate ΔV as the sum of the sub-volumes resulting from the product of the cell areas of the interpolated maps with the corresponding displacement value. This approximates the integral sum of the vertical displacement over the considered area. This integral sum is not significantly affected by cancelation of sub-volumes having opposite sign. Only in the case of UP4, the volume variation combines a possible small subsidence far from the center of the caldera and a larger uplift in the central part (see Fig. 7c, d).
We compare four ways of estimating ΔV for each of the studied episodes. The different approaches help to quantify the uncertainty related to the submerged sector of the caldera and to the large distal regions affected by a relatively small uplift. We obtain a first volume difference estimation (named ΔV1) by integrating the vertical displacement over the entire map, that is, the convex hull of all the RIM elliptic arcs. A second estimation (named ΔV2) is based on the integral sum over a circle of 6 km radius from the averaged center of symmetry (see Figs. 6,11a). This circle is consistent with the vertical displacement boundary in Fig. 13. Moreover, in UP4, UP5 and UP7, this excludes the regions in which vertical deformation is not significant with a 95% confidence. Then, we obtain a third estimation (named ΔV3) by calculating the volume in a circular sector approximately representing the northern side of Campi Flegrei caldera, where the density of measurements is higher (Fig. 13b), and assuming the same average volume/angle ratio over 360°. In particular, to obtain the total ΔV we multiply the resulting volume for 360°/140°. In substance, we are estimating the proportion of the ΔV that is captured by the GPS network onland and we use this to estimate the full ΔV. The center of this circular sector is the average center of symmetry calculated with the RIM, whereas the western and eastern boundary lines are defined by − 60° and 80° azimuth values (with a total angle of 140°, see Fig. 13b). The estimation ΔV3 does not rely on the interpolation of the vertical deformation occurred underwater, far from any GPS/leveling station before the recent placement of the GPS buoys Iannaccone et al. 2018). Finally, a fourth estimation (named ΔV4) considers the portion of this circular sector that is less than 6 km from the center.
Hence, for all the RIM maps described in the previous sections, we obtain four volume estimates reported in Table 1(A). Considering the leveling data, the volume variations ΔV3 and ΔV4 calculated on the northern circular sector are in general 12-15% larger than ΔV1 and ΔV2, respectively. Instead, the volume variations ΔV2 and ΔV4 in the circle of 6 km radius from the average center are slightly smaller than ΔV1 and ΔV3 (by about 2%). The volume changes based on GPS data include the effects of measurements errors. We assumed Gaussian errors and independence between the values randomly sampled at different stations. Under these assumptions, even if the vertical displacements are not always significant locally, the volume change is significant at 95% confidence, except for ΔV1 and ΔV3 in UP4. In particular, considering the mean values of the volume variation, UP4 shows a possible subsidence in the outer and northern part of the considered region, and this produces a 20% decrease from ΔV1 to ΔV3 and a 25% increase from ΔV3 to ΔV4. UP5 and UP6 are characterized by a wider region of uplift than other cases study, and this produces a 30% decrease from ΔV1 to ΔV2 and about 45% from ΔV3 to ΔV4. Moreover, ΔV3 and ΔV4 are about 50-60% larger than ΔV1 and ΔV2, respectively. Also UP7 shows ΔV3 and ΔV4 that are about 30-40% larger than ΔV1 and ΔV2. In general, the uncertainty reduces if the analysis is performed in the circle of 6 km radius from the average center, while the uncertainty increases if the analysis is performed on the northern circular sector.
We look at the changes in ΔV versus Δh max as a way to assess variations in source geometry and depth. If plotted against the maximum vertical displacement measured Δh max , approximated with the value measured at the RITE station, the volume variation data are significantly aligned, and Fig. 14 shows the linear regression least square fit: where the coefficients V0 and K are reported in figure and depend on the volume estimate considered ΔV1-4. In particular, K is about 0.26-0.27 × 10 6 m 3 /cm if based on ΔV1 and ΔV2, while it is 0.29-0.31 × 10 6 m 3 /cm if based on ΔV3 and ΔV4; the mean value of V0 is about 0.2 × 10 6 m 3 if (1) ΔV = V0 + KΔh max , Table 1 Volume variation of the ground obtained by the integral sum of vertical displacement ΔV1 considers the envelope of all RIM elliptic arcs, ΔV2 is the sum over a circle of 6 km radius from the averaged center of symmetry (see Fig. 6,11a), ΔV3 is based on the vertical displacement over a 140° circular sector on the northern side of Campi Flegrei caldera, multiplied by 360°/140°, and ΔV4 considers the portion of this circular sector than is less than 6 km from the center. Δh max is the maximum measured vertical displacement. The field of vertical displacement is obtained in (A) from the RIM; in (B) from InSAR acquisitions compared with least square fit (LSQ); and in (C) from least square fit only. Uncertainty range is based on 5th and 95th percentile values (~ 2σ) *InSAR processing method is detailed in "Appendix A" **Data from Del Gaudio et al. (2010)  based on ΔV1 and ΔV4, while it is 0.5 × 10 6 m 3 if based on ΔV3 and 0.1 × 10 6 m 3 if based on ΔV2. The substantial proportionality of ΔV and Δh max suggests that, apparently, the geometry and position of the volcanic source did not change with time (e.g., Amoruso et al. 2014a, b). Also, it suggests that the elastic mechanical behavior of rock did not change during the whole period considered. The value of Δh max in UP6 is not well explained by the linear regression if based on ΔV1 and ΔV3, that is, the calculated volume variation is larger than what the linear model would predict. According to a data processing of InSAR measurements detailed in "Appendix A," related to the ENVISAT SAR acquisitions on ascending (track 129, frame 809) and descending (track 36, frame 2781) orbits, we also calculated the volume variation of the ground during two minor uplift and subsidence events occurred in the period Dashed lines are 5th and 95th percentiles related to the measurement error on Δh max and the RIM uncertainty on ΔV when available, and to the confidence interval of the least square fit ΔV = V0 + K Δh max . A small box displays a zoom on the minor uplifts UP4-7 and shows uncertainty ranges on Δh max and ΔV. In plot a ΔV1 considers the envelope of all RIM virtual points, in plot b ΔV2 is obtained on a circle with 6 km radius from the averaged center of symmetry (see Fig. 6, 11a), in plot c ΔV3 is based on the vertical displacement over a 140° circular sector on the northern side of Campi Flegrei caldera multiplied by 360°/140°, and in plot d ΔV4 considers the portion of this circular sector than is less than 6 km from the center. All the plots show the mean value and the uncertainty range of the coefficients K and V0 11/2004-12/2007. These InSAR measurements are limited to the inland sector, so we produced volume variations with the methods ΔV3 and ΔV4. We report them in Table 1(B), and we compare them with the statistical inference in Eq. (1). We always approximated the Δh max by using the value measured at the RITE station. In this case, the value is closer to the maximum uplift in the InSAR data than the measurement error. Both the estimates of ΔV3 and ΔV4 are statistically consistent with the linear model, especially ΔV4, which does not include noisy data in the outer portion of the region.
Thanks to Eq. (1), we finally estimated the volume variation of the ground also for the uplifts in the periods 1950-1952, 1968-1972 (Δh max  We report them in Table 1(C). Under the assumption of this least square fit, the full time series of volume variation could be inferred from the time series of vertical uplift reported in Fig. 1b. The total volume variation ΔV4 corresponding to the historical maximum Δh max = 380 cm measured in 1985 from the minimum of 1950 is equal to [112 ± 3] × 10 6 m 3 .

Discussion
The presented interpolation method (RIM) allowed a detailed analysis of the spatially sparse data of ground deformation obtained from GPS, GPT and leveling measurements. In particular, RIM was applied to investigate the ground deformation patterns in the volcanic area of Campi Flegrei during the last 39 years. The resulting RIM maps of vertical and horizontal displacement components represent a useful tool to better understand the ground deformation history of this area, particularly before 1992 when space-borne remote sensing data were not available. The new method has also promising applications for retrospective studies in other caldera systems, especially in the context of counterfactual analysis of historical unrests, i.e., the stochastic modeling of past crises at the particular volcano of concern that had the potential for a dangerous event but did not ultimately result in a significant eruption (Selva et al. 2012(Selva et al. , 2015Hincks et al. 2014;Woo 2018, Aspinall andWoo 2019).
In fact, improving the spatial reconstruction of horizontal displacement is a matter of considerable interest in hazard assessment. It is a good proxy of the ground surface extension, which may eventually lead to an increase in the secondary permeability by fracturing and hence of degassing and, in extreme cases, to the opening of new volcanic vents/ fissures. Numerical models and experimental studies showed that dyke propagation is influenced by the local properties of the stress field and by preexisting fractures (Gaffney et al. 2007;Maccaferri et al. 2011;Le Corvec et al. 2013;Corbi et al. 2015Corbi et al. , 2016Rivalta et al. 2019). For the specific case of Campi Flegrei, the maps indicate that the maximum horizontal displacements are located over a semicircular annular region, with a radius of about 2-3 km from the town of Pozzuoli. This pattern is caused by the 3D radial nature of the deformation acting vertically close to the center. Figure 15 shows that faults with metric displacements and clastic dikes, not related to volcanic vents, are exposed in the annular area depicted by the maximum horizontal displacements, together with quite high values of fracture density (Vitale and Isaia 2014;Isaia et al. 2015;Bevilacqua et al. 2015. These findings suggest that the area was subject to extension at least since Epoch III of its activity (5.2-3.5 ka; Smith et al. 2011;Bevilacqua et al. 2016Bevilacqua et al. , 2017. In the RIM maps related to the bradyseismic crisis of 1983-1984 (Figs. 11a and12a) and episode UP5 (Fig. 8), the peak horizontal displacement is observed in the centerwest sector of the annular region. Instead, in the episodes UP4, UP6 and UP7 (Figs. 7,9 and 10) it is apparent that the spatial peak of the horizontal displacement, and hence the maximum extensional deformation, affected the center-east sector of the annular region, where most of the seismic and fumarolic activity of the Campi Flegrei caldera is currently observed (INGV periodic bulletins of Campi Flegrei; http:// www.ov.ingv.it/ov/it/campi -flegr ei/275.html). However, the change of the area of maximum extensional deformation may be related to the differences in the sites and methods of measurement and would require additional data to be confirmed.
On the base of quasi real-time geodetic measurements of the horizontal displacement, the new method and the associated maps could represent a useful input in the definition of short-term vent opening maps. Indeed, although long-term vent opening maps (Bevilacqua et al. 2015(Bevilacqua et al. , 2017 represent a probability field over the region, they do not consider monitoring information. The map of horizontal displacement could represent additional information in the construction of a vent opening forecast map, under the assumption that the source of the current deformation is going to affect the eruption occurrence and location. Similar considerations could be valid for the vertical displacement maps. The statistical combination of these two layers of information, i.e., a longterm map and a deformation map, is beyond the scope of this study, but it represents a promising task for future research (Bevilacqua et al. 2018;Patra et al. 2019).
A second output of our analysis is the reconstruction of the spatial distribution of the vertical component of ground deformation. In Campi Flegrei caldera, the resulting uplift patterns can be either quasi-circular (e.g., UP6 and UP7) or more elliptical in shape (e.g., UP4 and UP5). A nearly circular shape is also observed for the vertical displacement patterns related to the 1980-1983 and 1983-1984 periods. All the eruptive vents active in the last 5500 years are located within or near the edge of the vertical displacements, indicating that this quasi-circular area corresponds to the most active portion of the caldera, as suggested by different authors (e.g., Barberi et al. 1991;Capuano et al. 2013;Kilburn et al. 2017). The boundary of this area (see Fig. 13a) fits well with the data obtained from the gravity field horizontal derivative map of Campi Flegrei and the density anomaly maps at different depths (Capuano et al. 2013). This is also in agreement with the long-term distribution of the vent opening probability maps of the caldera as derived by Bevilacqua et al. (2015Bevilacqua et al. ( , 2017. The vertical uplift also enabled the calculation of volume variation of the ground for the analyzed periods and its approximated estimation even for not analyzed periods, on the base of linear regression. Indeed, a major interest of the method is to give the order of magnitudes of surface deformation volume, without having to perform source parameters inversion/optimization. On the other hand, this method should not be used to provide constraints for optimization purposes. These volumes may or may not be equivalent to the volume change at the source of deformation, depending on its geometry and depth, but represent significant information that could be used in further data inversion. The substantial proportionality of ΔV and Δh max suggests that the geometry and position of the volcanic source did not change with time, and that the elastic mechanical behavior of rock did not change during the whole period considered. Finally, we remark that, even though the RIM has been primarily devised for retrospective hazard analysis of past crises for which remote sensing data were not available, the RIM results could be applied to the latest satellite data available through the Sentinel-1 constellation from the European Space Agency (Guglielmino et al. 2011;Bonforte et al. 2013). For example, the method could be helpful where InSAR data are significantly sparse, that is, where most of the InSAR information is incoherent or scattered in disconnected patches or pixels (e.g., in vegetated areas). The method could provide enriched data to test analytical source models. However, this presents nontrivial technical complications which solution is outside the scope of this work. The additional uncertainty related to the interpolation and the local assumption of radial symmetry should be considered carefully in the inversion process. A possible strategy may be to employ the family V of virtual displacements, without performing the secondary interpolation.

Conclusions
Whenever there is reasonable evidence of radial symmetry between each pair of measurements, the new interpolation method RIM presented in this study enables the fast construction of ground displacement maps for both vertical and horizontal components. These maps use and rely only on ground-based measurements, which provide the longest temporal datasets in the past, and could be updated at a faster rate than satellite image acquisitions, even during a rapidly evolving crisis. The method includes uncertainty quantification related to the propagation of measurement error, as well as to the choice of the interpolation graph. The RIM is applied to diverse examples of uplift and subsidence episodes in the active Campi Flegrei caldera, either collected in the period preceding space-borne SAR missions, or more recently. This covers the bradyseismic crisis in 1982-1984, the subsequent subsidence in 1985-1989 and four minor uplift episodes observed in 2011-2013. Significantly, our results are independent on any assumption on the geometry, location and physical properties of the source except for requiring a locally radial pattern, i.e., allowing multiple centers of symmetry. We remark that our study is not equivalent to the ground deformation data obtained by using an analytical source model. The purpose of the RIM is therefore not to reconstruct the source, but only to exploit the additional information on ground deformation provided by the simple assumption of radial symmetry between each pair of nearby measurements. Nevertheless, if radial symmetry is not verified, the RIM results can be inaccurate. In fact, while a source model can physically reconstruct a full deformation profile, also along the radial direction, the key step of the RIM is an interpolation orthogonal to the radii. We envisaged that the new method has potential applications in: 1. The retrospective analysis of historical unrests, included those related to past crises that did not ultimately result in an eruption; 2. The efficient comparison and combination of remote sensing and ground-based measurements; 3. The production of enriched fitting data for the application of a variety of analytical and numerical source models; 4. The production of short-term hazard assessments based on the statistical combination of long-term vent opening maps and monitoring data. Fig. 15 a, d- 1972-1974, 1983-1984, and 2005-2015 (data from Kilburn et al. 2017), the boundary area of the presently active caldera (blue) and the maximum horizontal displacement area (yellow, ideally extended in the Gulf of Pozzuoli). UTM WGS84, Zone 33 N coordinate system (m) ◂ Finally, we summarize the main volcanological findings obtained from the application of the new method to the Campi Flegrei caldera as: • The interpolated maps of horizontal displacement component show maximum values in a semi-annular region located onland at about 2-3 km from the center of deformation (town of Pozzuoli), in all the uplift episodes considered. Faults with metric displacements, clastic dikes and high values of fracture density are observed within this region. This extensional area includes the Solfatara crater and the Pisciarelli locality, where most of the fumarolic activity is currently observed. The apparent location of the peak horizontal displacement values inside the semi-annular region depends on the considered uplift episode and can either be in the center-west or in the center-east area of the caldera. The semi-annular region could be ideally extended into the Gulf of Pozzuoli, based on the distribution of the epicenters of the earthquakes in the period 1972-2015, thus depicting a quasi-circular annular region. • The interpolated maps of vertical displacement provide a new estimation of the volume variation of the ground. Extrapolating the volume/angle ratio estimated in the northern circular sector over 360°, the result is generally larger than what obtained using the interpolated map over the Gulf of Pozzuoli. A linear regression is appropriate to compare volume variation ΔV and maximum vertical displacement. The multiplicative factor K in the linear relation is about 0.3 × 10 6 m 3 /cm if we estimate the proportion of the ΔV that is captured by the GPS network onland, and we use this to estimate the full ΔV. In this case, the 95% confidence interval on K because of linear regression is ± 5%. The linear model is consistent with InSAR data of the period 2004-2007. The model allows us to broadly quantify the volume variation during the uplift events of 1950-1952 and 1968-1972, supporting the hypothesis that the elastic mechanical behavior of the rock held during the whole period considered. Isolines of vertical deformation can either be circular or be elliptical with axes in NNE and NWW directions as observed in UP4 and UP5. Moreover, the slope of vertical uplift defines a line at the boundary of the well-known domelike shape, at about 6 km from the average center of deformation. Notably, the eruptive vents of last 5.2 kyr BP (i.e., Epoch III and Monte Nuovo eruption) are apparently contained between the semi-annular region depicted by maximum horizontal displacement and this external boundary based on vertical uplift. networks of the Future: INGV 2.0," funded by the Ministry of Education, University, and Research (Italy); the project "Collaborative research: Using precursor information to update probabilistic hazard maps," funded by National Science Foundation (USA), Grant #1821311; the Contract INGV-DPC 2019-2021, Annex B2, Task 2: "Development of a real-time monitoring system for the ground deformation in the Neapolitan volcanic region using HR-GNSS measurements, and development of statistical and numerical models for the short-term vent opening probability in the Campi Flegrei caldera." The manuscript does not necessarily represent official views and policies of the Dipartimento della Protezione Civile. The authors acknowledge the thorough reviews of Nico Fournier and of five anonymous reviewers.
Author Contribution SV conceived the conceptual idea of a multipolar interpolation. AB conceived the uncertainty quantification approach and set up the theoretical formalism. AB and SV developed the method, implemented and performed the simulations, interpreted the computational results and wrote the paper. SV and ANov performed the processing of InSAR data. All authors discussed the results, commented on the manuscript, provided critical feedback and gave final approval for publication.

Data Availability Statement
The datasets necessary to interpret and replicate the results in this article are all included as Supporting Information. In particular, the graphs and the analyzed displacement values of GPS data, the GPT data and the leveling data are in Supporting Information SI4, SI5 and SI6, respectively. For the sake of completeness, we include the GPS time series from De Martino et al. (2014a). We also include the results of a sensitivity analysis performed over saturated graphs (SI2) and the 5th and 95th percentiles of uncertainty of U4-U7 (SI3).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Appendix A. Volume variation of the ground based on InSAR acquisitions
In order to compare our estimates of ground displacement volume with independent estimations, we exploited also the InSAR results generated from the ENVISAT SAR acquisitions on ascending (track 129, frame 809) and descending (track 36, frame 2781) orbits and covering two minor uplift and subsidence events occurred in the period 2004-2007. Although a more extensive comparison would be feasible on other time windows, we limited our analysis to this preliminary example, because the required data processing is not trivial, and the uncertainty affecting the displacement is significant. Additional research aimed at the further reconstruction of the vertical and horizontal deformation from the processing of remote sensing data in combination with ground GPS is a challenging task that is outside the purpose of this paper (e.g., Guglielmino et al. 2011;Bonforte et al. 2013). The SAR images were processed by Persistent Scatterer (PS) Interferometry techniques (Costantini et al. 2017) and made available through the Extraordinary Plan for Environmental Remote Sensing funded by the Italian Ministry for the Environment, Land and Sea (Di Martire et al. 2017). For each PS, InSAR can only measure a projection of the real 3D deformation vector along the satellite Line Of Sight (LOS). Therefore, only by exploiting the deformation measures taken over the same point from two different acquisition geometries (namely the ascending and descending orbits), it is possible to retrieve the vertical and east-west directions of displacement, as shown in Dalla Via et al. (2012).
In particular, we generated a 100 × 100 grid and we assigned for each PS the mean value of the displacement data included in a circle with a radius of 300 m. Subsequently, we interpolated these data by a linear interpolation method. The obtained maps (Fig. 16) are similar to those constructed by RIM and show a radial pattern of the vertical displacement both for the uplift and subsidence phase. We calculated the volume variation of the ground (ΔV expressed as millions cubic meters [Mm3]) in the two bradyseismic episodes by the definitions of ΔV3 and ΔV4 previously explained, because these InSAR measurements are limited to the inland sector only. The obtained results are reported in Table 1(B).