New Insights on Ground Deformation at Campi Flegrei Caldera Inferred From Kinematics and Dynamics Investigation of Borehole Tilt

The study of the ground deformation in active calderas provides valuable indications of the ongoing dynamical processes. In this framework, we analyze the borehole tiltmeter data recorded at Campi Flegrei (Southern Italy) from their first installation (1 April 2015), in order to retrieve the kinematics of the ground deformation and its main anomalies. Our approach involves the nonlinear dynamics analysis aimed at the identification of the degrees of freedom of the system and thus its complexity. Starting from the second part of the year 2020, the behavior of the dynamical system becomes collective, and a locally significant deviation of the deformative pattern from the background trend is observed. In particular, a series of 15 slow tilt variations (like jumps lasting a few days) appears in the ground deformation. They are associated with a very low‐dimensional system likely generated by a local second‐order source. The latter is related to fluid migration and it is superimposed on the primary one located in Pozzuoli town and driving the current uplift of Campi Flegrei. The effects of such a local, second‐order stress field are, indeed, evident in the area of the volcanic dome of Mt. Olibano, where they are measured as tilt steps. The superposition of primary and second‐order sources of deformation acting contemporaneously can justify the recent variations in geophysical and geochemical parameters. Our approach based on the joint dynamical and kinematical analyses of the tiltmeter data can be applied to other volcanic/tectonic areas or extended to other geophysical and geochemical variables.

Nonlinear behavior often characterized Earth's crustal deformation too. For example, high dimensional chaotic dynamics and synchronization phenomena have been observed in the GPS time series recorded in New Zealand, reflecting the high nonlinearity of the underlying plate tectonic processes (Hobbs & Ord, 2018). Nonlinear cyclical or quasi-cyclical dynamics characterized ground deformation in the Bucharest metropolitan area; such a pattern has been attributed to a tectonic stress field generated by NW-SE transpressional movements (Armaş et al., 2017). Recent analyses of GPS in Akutan and Okmok volcanoes and Piton de la Fournaise (La Reunion) also show a nonlinearity mechanism acting in the plumbing system (Walwer et al., 2022).
At Campi Flegrei caldera, the analysis of tide-gauge elevation data recorded between 1970 and 1990 shows that, in coincidence with the uplift crisis of 1982-1984, the system underwent a critical transition between phase-locking with tidal forces and low-dimensional chaos (Cortini et al., 1991).
In the present paper, we analyze the tiltmeter time series recorded at Campi Flegrei from 1 April 2015 to 30 April 2022, in order to retrieve the kinematics of the ground deformation and its main anomalies. In this area, the ground deformation, whose centroid is located in Pozzuoli , is a relevant phenomenon related to bradyseism and uplift (Del Gaudio et al., 2010). Significant ground deformation resulting in more than 150 cm of vertical uplift occurred during the crisis of 1982-1984. Since the last months of 2005, a new unrest started with an increasing rate of uplift (De Martino et al., 2021;Polcari et al., 2022) (about 90 cm up to now), seismicity (Bellucci Sessa et al., 2021), and degassing activity . Uplifts are related to both magmatic and hydrothermal sources (Cannatelli et al., 2020;Lima et al., 2021). Hot fluids exolve from a magmatic body in the deeper part of the plumbing system (∼8 km), accumulating at intermediate depths in a sill-like reservoir (3-4 km), and feeding the shallow (<2 km) hydrothermal system (Amoruso & Crescentini, 2022;Siniscalchi et al., 2019). A large amount of geothermal fluids is released at the surface through diffuse degassing at Solfatara-Pisciarelli fumarolic fields (Siniscalchi et al., 2019;Troiano et al., 2022). This hydrothermal activity is concentrated along faults and fractures, which are the main pathways for the fluid escape. Uplift episodes are accompanied by volcano-tectonic earthquakes (Bellucci Sessa et al., 2021), often occurring in seismic swarms. Most events are located at shallow depths (<4 km b.s.l) beneath the Solfatara-Pozzuoli area and are related to brittle shear failure mechanisms induced by the fluid FALANGA ET AL.

10.1029/2022EA002702
3 of 16 pressurization (Ricco et al., 2019). The close link between ground deformation, seismicity and geofluid circulation testifies the complexity of the Campi Flegrei volcanic system, which makes an accurate realistic modeling of the underlying processes difficult. As example, the modeling of the observed primary ground deformation (from SAR, GPS, and gravity data) has required the combination of simpler sources at deeper and shallower depths, alternating deflation and inflation phases (see, e.g., Amoruso & Crescentini, 2022;Camacho & Fernández, 2019;Castaldo et al., 2021;Samsonov et al., 2014;Tiampo et al., 2017, and references therein). Thus, such a complex system deserves effective methodologies of data analysis which can reveal the signature of different source processes.
In line with these thoughts, we characterize the deformation trend by adopting a nonlinear approach based on the reconstruction of the phase space. We find evidence of low dimensional deformation sources (and background trend) acting on different time scales. Looking at kinematic data from a nonlinear perspective allows getting more insight into the source process of ground deformation, providing valuable information on the actual dynamics state of Campi Flegrei caldera.

Materials and Methods
To take into account the eventual nonlinearity underlying experimental time series, appropriate tools are required, such as Average Mutual Information; False Nearest Neighbors and Grassberger and Procaccia's method.

Average Mutual Information
The Average Mutual Information (AMI) is the counterpart of the autocorrelation function in the nonlinear domain (Fraser & Swinney, 1986), which is introduced to evidence statistical dependencies in the data. It is defined as follows: (1) where P(x(i)) is the probability to measure x(i) and P(x(i),x(i + L)) is the joint probability to measure x(i) e x(i + L); L is related to the time lag τ. If x(i) and x(i + L) are independent, the joint probability is factorized into the product of the individual probability and I(L) goes to zero. In real cases, a certain amount of correlation among the variables always is present and the receipt is to take the first minimum of AMI as a function of L. If a first clear minimum is not found, L is chosen so that I(L) = I max /5 = I(0)/5 (Abarbanel, 1996).

False Nearest Neighbors Technique
Once the time lag has been estimated, it is possible to obtain information on the lowest embedding dimension by using the False Nearest Neighbors (FNN, Kennel et al., 1992). The method is based on the projection of the attractor onto subspaces of increasing dimension until it reaches the dimension for which the attractor is completely unfolded. In fact, if the space on which it is projected has a dimension lower than that of the attractor, the topological structure is no longer preserved. Some points are close (neighbors) to others due to the projection effect; this ambiguity is completely resolved in a higher dimension space. These points are called false neighbors. As the embedding dimension changes, the fraction of the false neighbors is calculated. The value of m for which this fraction goes to zero corresponds to the optimal. A brief description of the algorithm is the following. For each point of the time series we look for the next neighbors in a space with m dimensions. Then the distance − is calculated according to a certain metric. If the metric is Euclidean, then: Going to the m + 1 dimension: FALANGA ET AL.
10.1029/2022EA002702 4 of 16 If 2 +1 ≫ 2 then the value of the previous distance was false due to the projection onto subspace of insufficient embedding dimensions to unfold the attractor. Iteratively, we introduce the parameter: If is greater than a fixed threshold , the point is marked as having "false neighbors." Thus, the fraction of points is calculated for which > as a function of . The value of for which this function goes to zero is taken as the embedding dimension.
It is important to note that a deterministic dynamics is characterized by the fact that the curve goes to zero for a certain value of m and then remains zero. Low dimensional systems have low values of (a few time scales acting in the dynamics, in the range 2-4), whereas high dimensional systems experience high values of (many time scales acting in the dynamics, of the order of 5-8, see, Abarbanel, 1996). The embedding dimension provided by FNN is the lower limit for the dimension of the reconstructed phase space, as well as the upper limit for the attractor dimension. This will be a useful guide in estimating the dimension (i.e., the number of degrees of freedom) related to the experimental signals.

Grassberger and Procaccia's Method
The dimension is an important parameter that reflects the topological complexity of an attractor. It basically represents the amount of information necessary to specify a point on the attractor, that is, it provides approximately the number of degrees of freedom of the system.
Once the embedding space was reconstructed, Grassberger and Procaccia (1983) showed that the self-similarity properties of an attractor can be characterized by estimating their dimension (entire or fractal) through the correlation integral. Once obtained the m-dimensional signal ( ) = 1, . . . , with = − ( − 1) (M is the number of points in the embedding space obtained starting from the points of the scalar series original and is the appropriate time lag) the following quantity is defined as the integral of correlation: where is the embedding dimension, is the length of the hypercube (or the radius of the hypersphere) covering the m-dimensional space.
In the limit → 0 , → ∞ , and, in principle noiseless, ( ) scales versus : By inverting, the dimension 2 (an estimate of the attractor dimension a ) can be found: This dimension is called correlation dimension (Hilborn, 1994;Ott, 1993).
Note that the correlation dimension allows distinguishing between signals generated by deterministic dynamics and those generated by a stochastic process. In fact, given an experimental signal, if the correlation dimension is found to be smaller than the phase space dimension then the signal is supposed to be generated by nonlinear, deterministic dynamics. If, on the other hand, ( ) scales like , the dynamics are mainly stochastic. It is clear that, in linear cases, the fractal dimension is reduced to the topological dimension. There is, therefore, a first characterization of the dynamic systems according to the value assumed by 2 : it provides an indication if the experimental signals can be considered a linear superposition of normal modes, or it is essential to move in the direction of nonlinear, and/or possibly stochastic, systems.

Tiltmeter Data
In the present work, we analyze data from the tiltmeter network at Campi Flegrei. Specifically, we consider the three stations (CMP, ECO, and HDM) equipped with borehole tiltmeters operating since 2015. They record the uplift along NS (North-South) and EW (East-West) components, whose resultant is a ground tilt in the radial direction from the centroid of deformation. A first overview of the six signals recorded over a time span from April 2015 to April 2022 (more than 7 years) shows the broad similarity among the time series ( Figure 1). However, two abrupt tilt variations at ECO occur at the end of April 2016 and in mid-July 2018, and further minor changes affect HDM NS and ECO EW since mid-September 2020. Comparing the tilt with seismicity (Ricciolino & Lo Bascio, 2021), we notice that the greatest release of seismic energy occurs after the tilting reversal in July 2018, and since mid-September 2020 the earthquake intertime shortens while the tilt rate increases (area highlighted in yellow in Figure 1). In the same figure, grey lines represent earthquakes occurrence, also reported in Table 1.

Results
In each subsection, we briefly describe the main outcomes derived by the kinematics and dynamics analyses. Then, we provide our interpretations of the reconstructed deformation pattern in a conceptual scheme.

Kinematics of the Ground Deformation
The kinematics at the three sites can be visualized vectorially in a hodograph, in which the color curves represent the progressive tilt change recorded at each site ( Figure 2). The tilting directions of CMP and HDM (located respectively NNW and ESE of Pozzuoli) are generally congruent with those obtained from the deformation field induced by the Pozzuoli uplift (De Martino et al., 2021). In contrast, the ECO station, located NE of Pozzuoli, tilted consistently with the uplift only from March 2017 to July 2018. Before this time interval, its tilting direction was NW but, after the tilting reversal in April 2016, it became W. The second tilting reversal occurring in mid-July 2018 puts an end to the axisymmetric deformation of ECO which, from then on, tilts increasingly to W.
From January 2019 until September 2020, ECO changes direction again tilting toward N. Since then, a higher tilt rate is observed than in previous years, and ECO tilting direction also changes, following that of CMP, synchronizing with its trend in both azimuth and amplitude (highlighted with colored ellipses in Figure 3). By analyzing in more detail the increase in tilt rate since September 2020, highlighted in Figures 3a and 3b by the isochronous hatched curves of the extrapolated cumulative tiltmeter change, interesting details emerge. First, the increase in tilt can be quantified, at CMP, as 4.6 μ radians/month in the NNW direction from 11/09/2020 to 30/4/2022 compared with 3 μ radians/month in the same direction from 01/01/2020 to 10/09/2020. Second, a counterclockwise rotation of the tilting direction is observed first at ECO since 11/09/2020 and then at HDM since 22/04/2021. Third, the latter site exhibits a further interesting behavior: since July 2020 it shows a steplike pattern of the tilt evidenced by the occurrence of 15 slow tilt variations in the NNE direction ( Figure 3c and Table 2). These anomalies generate an asymmetrical areal deformation compared with that which would occur if they were not present (correct hodograph in Figure 3b). A zoom of the hodograph at HDM is plotted in Figure 3c. It shows 15 tilt changes: the first four increments concern only the NS component while from the fifth onwards the variation also involves the EW component. Their mean value is 3 μ radians in the N30.9°E direction with a period of 3.9 days. In comparison, the median values are respectively 2.6 μ radians in the N32.1°E direction with a period of 4.4 days. The intertimes between the anomalies have an average duration of 42.1 days (median value = 44.5 days) and tend to decrease over time. Furthermore, during the step-like pattern observed at HDM, the tilting direction slowly tends to align itself to the North (counterclockwise rotation) until November 2021 and then rotates in the clockwise direction. During the 15 polarized tilt increments, seismicity is rare, whereas it occurs systematically between one increment and the next (red dots with date and magnitude of the earthquakes).
A support to the evidenced unusual deformation detected by HDM is provided by another tilt station (OLB) close to it, which partially detected the same deformation pattern. OLB is 390 m away from HDM, which is a borehole sensor located on Mt. Olibano at a height of 112 m asl (−25 m from ground level), whereas OLB is a surface tiltmeter installed in a tunnel at the base of the volcanic dome (10 m asl and −1 m from ground level). The NS component of OLB recorded at least 9 of the 15 tilt anomalies (Figure 2 and Table 2).
The kinematics observed at HDM is compatible with the superposition of two processes. The first process is the response of the site to the bradyseism-related uplift that induces HDM to tilt to the ENE and is represented by the jump-corrected hodograph superimposed on the original one. The second process is probably due to the perturbation caused by a volcanic structure that, for more than a year, has been generating a jerky stress field in the proximity of the volcanic dome of Mt. Olibano. Such a field is producing a well-polarized deformation in the NNE direction (toward Solfatara).

Tidal Fortnightly Ground Tilt
Past studies of the ground tilt at Campi Flegrei show evidence that medium and long-period tidal constituents (diurnal, fortnightly, and monthly) induce ground oscillations that are superimposed to the normal deformation trend of the area (De Lauro et al., 2018). In particular, the analysis of tiltmeter time series filtered in the fortnightly (Mf) tidal band (Petrosino et al., 2020;Ricco et al., 2019) (2,6,7,9,10,11,12,13,14) were recorded simultaneously at OLB.
FALANGA ET AL.
10.1029/2022EA002702 9 of 16 For the first time, the LAOs of December 2019, April 2020, and March/ April 2022 (likely associated with an intense seismic activity that climaxed with earthquakes of maximum magnitude Md 3.1, 3.3, and 3.6, respectively) are also detected at HDM tiltmeter ( Figure 4a). Moreover, after the April 2020 LAO, HDM shows an amplitude modulation never observed before; since July/September 2020 the azimuths of ground oscillations in the Mf tidal band are pretty stable and roughly oriented in the NNE-SSW direction (Figure 4b). The comparison of the azimuths calculated between April 2015 and June 2020, mostly oriented along ESE-WNW, with those of July 2020-April 2022, shows a clear rotation of about 90° of the main tilting direction (Figure 4c). This breaks the exceptional stability of the azimuthal pattern observed throughout more than 5 years (Petrosino et al., 2020). In addition, starting from July 2020, the azimuth patterns at CMP and ECO appear less spread, with dominant direction along E-W at CMP. ECO also showed a slight rotation of the tilting plane from ENE-WSW to SSW-NNE, actually "synchronizing" it with the oscillation plane of HDM.

Nonlinear Dynamics of the Ground Deformation
The observations related to the kinematics combined with the simple-decaying waveforms of LAOs sustained for a long time in Mf band, and the occurrence of sudden jumps, especially at HDM NS since July 2020, evidence the presence of a complex dynamical system, whose level of complexity has to be evaluated. This can be achieved by estimating the degrees of freedom necessary to be taken into account to fully develop the dynamics. The recognition and the description of the dynamical system underlying the tilt along the time could shed light on the complexity of the system associated with the different observed events, giving information on the source mechanism and its evolution. On the other hand, we expect that the coupling between fluid flows and solid structure can produce nonlinear signals, also in the form of self-oscillations, meaning that the system can show a collective behavior, described by a few degrees of freedom (De Lauro et al., 2008;De Lauro, De Martino, Falanga, & Ixaru, 2009;De Lauro, De Martino, Falanga, & Palo, 2009).
Moving in the framework of the theory of dynamical systems, we adopt the standard procedures for reconstructing the asymptotic dynamics in suitable phase space from experimental series, by exploiting the time delay method (Abarbanel, 1996;Takens, 1981). The estimate of the suitable time lag used to construct the independent vectors in the reconstructed phase space is provided by AMI (Fraser & Swinney, 1986) able to catch the nonlinearity eventually present in series. A first indication of the dimensionality of the phase space is given by the application of the FNN (Kennel et al., 1992). The asymptotic dynamics embedded in this reconstructed phase space is then characterized by a proper attractor whose dimension could be also fractal. The latter is estimated by means of Grassberger and Procaccia method (1983). Notice that nonlinear systems can have attractors with a fractal dimension also when the dynamics is not chaotic.
The Non Linear Analysis (NLA-hereinafter) described in Section 2 was separately applied to all the series from the three stations starting from 01/04/2015 to 30/04/2022, along both directions NS, EW. The series were first filtered with a pass-band filter in the range [0.5-90] d and then partitioned in time periods following the hydrological cycle (about 14 semesters).
In Figure 5, the results related to all the stations are reported, where all the steps in the fractal dimension estimate are shown. The time evolution of the dimension reveals a transition mechanism that brings the system from higher dimensions (in the range 3-4) to self-organized oscillations characterized by a very low fractal dimension (about 1.5 or even less), highlighted at HDM NS since October 2020 and lasting also in the last semester. Indeed, it is noticeable the downing of the correlation dimensions versus embedding dimensions curves to the last time periods (last traces). It indicates that the system is undergoing a synchronised regime. In detail, the superposition of the correlation integral for all the traces is shown in Figures 5c1-5c3 and 5d1-5d3: the common linear region is used to estimate the fractal dimension. It is possible to observe at small scales steeper slopes, which indicates that the fluctuations are dominant at these scales. The curves of the FNN versus embedding dimensions are plotted in Figures 5e1-5e3 and 5f1-5f3: the percentage of FNN equal to zero indicates the suitable embedding dimension to be used to fully unfold the attractor. As it can be seen, relative to HDM in Figure 5e3, that is, regarding the NS component, the lower embedding dimensions are relative to the last traces, indicating a progressive synchronisation of the phenomenon toward a few degrees of freedom system. The inset shows the superposition of all AMI functions, whose first global minimum of around 300-400 hr provides an indication of the suitable time lags to be used for the Grassberger and Procaccia method application. Minima between 300 and 400 hr are, in some sense, associated with dynamics that evolve on time scales of the order of 45-60 days. The final result of the embedding estimates over time is displayed in Figures 5g1-5g3 and 5h1-5h3. For the NS component, as already indicated by FNN, it is clear that starting from October 2020, the system is more regular. The existence of the plateau suggests that the system is deterministic and nonlinear with different levels of complexity. The stability of the plateau over several embedding dimensions indicates a stable system; in particular, in Figure 5g3, the fractal dimension of about 1.5 evidences a nonlinear behavior of the low-complexity system in a dynamical regime producing self-sustained oscillations (cyan and dark red waveforms). HDM EW has a more regular behavior with dimensions in the range 1.5-2 independent of time ( Figure 5h). The other stations show higher fractal dimensions.
The analysis so far discussed is summarized in Figure 5i, in which we report the evolution of the estimated fractal dimensions at the three stations and for both the directions of motion. Further, such evolution is compared with the tilt as acquired by ECO EW, which well captures LAOs (Figure 5l). On average, the system is always low dimensional with dimensions lower than 4. CMP EW appears to have a dependence on the hydrological cycle and  Petrosino et al. (2020), who found clear evidence of an amplitude increase during the spring-summer time due to site thermoelastic effects. In addition, in the monthly and fortnightly tidal bands, the azimuths of the tilting planes switch their orientations according to the hydrological cycles. In correspondence with the greatest LAOs, the fractal dimension is very low (about 1.5-see the ellipses in Figure 5i in correspondence with LAOS in Figure 5l), indicating that the mechanism generating LAOs can be described by a very few degrees of freedom, underlying globally organized oscillations with respect to the background. In support of this result, the relative FNNs indicate embedding dimensions at most three, reminding us that this value corresponds to the dimension of the space in which the attractor has to be fully embedded.
In light of the previous results, that is, looking at the kinematics on raw signals, the evolution of the tilt in Mf frequency band, and the dynamical approach, it is interesting to notice that the tilt from HDM NS reveals a relevant variation with respect to the background, starting from the second semester of 2020 (see, Figures 3 and 4a) and still ongoing.

Discussions and Conclusions
The kinematics analysis evidenced an increase in the tiltmeter rate since September 2020. At the same time, the tilt direction of ECO rotates counterclockwise, synchronizing with CMP thus providing the same azimuth and amplitude. This mechanism may be the manifestation of a synchronization phenomenon obtained by means of phase locking (Pikovsky et al., 2001). In addition, at HDM, a peculiar deformation pattern consisting of a series of steps in NNE direction appears, eventually causing a smooth counterclockwise rotation since July 2020, which accelerated in April 2021.
In the Mf tidal band, the nearly regular seasonal LAOs, usually observed at ECO, after April 2020 disappear, and at HDM, we observe a 90-degree rotation of the tilting direction since July/September 2020.
NLA performed all over the time suggests that starting from October 2020 the system is making experience of more synchronized dynamics, especially evident at HDM NS.
Putting all the information together, we can hypothesize that the ground deformation, measured as steps by the borehole tiltmeters, can be the result of two sources: the main source (always active) and an additional second-order source, which locally affects the deformation.
The main source, which causes the uplift centered in Pozzuoli, gives origin to the background trend which at HDM results as a tilt mainly oriented toward ENE (thus mainly detected on the EW component). The main source corresponds to a medium-low dimensionality (<4). At CMP and ECO, the ground tilt related to the uplift trend is almost equally distributed on the NS and EW components due to the radial position of the source, providing similar attractor dimensions. NLA suggests that the main source is associated with a nonlinear dynamical system with few degrees of freedom. The background trend of the ground tilt corresponds to the primary ground deformation modeled in past studies (see, e.g., Samsonov et al., 2014;Tiampo et al., 2017;Castaldo et al., 2021, Amoruso & Crescentini, 2022.
At least up to April 2020, the LAOs superpose to the general trend as self-oscillations in Mf band, reducing the system's dimensionality down to 1.5, indicating a collective phenomenon, likely related to the global response to an external forcing of hydrological/tidal origin.
The second source, which was likely activated after April 2020, is associated with a nonlinear dynamics of very low dimensionality (∼1.2 at HDM NS). The low dimensionality indicates that the dynamics of the system has Olibano. This source generates the tilt steps, which are well polarized in NNE direction, and it is likely due to fluid accumulation east of a transfer structure identified by the polarization analysis of the seismic noise and oriented in the NE-SW direction (Petrosino & De Siena, 2021). Indeed, seismic noise polarization provides evidence of fluid redistribution in the Campi Flegrei area: while in October/ November 2019, the unpolarized anomaly due to the underground fluids was diffuse almost around the primary deformation source in Pozzuoli, in May/June 2020 it was at about 1 km east of its original location, thus evidencing fluid migration.
We hypothesize that the fluid accumulation toward the east sector continued after May/June 2020. In that case, it is likely that such fluids could have generated supplementary forces giving rise to a variation of the stress field and inducing the rotation of the ground tilt toward west, as observed at ECO and HDM. Indeed, the area between Mt. Olibano-Solfatara is much more fractured and thus prone to deformation (Isaia et al., 2015;Vitale & Isaia, 2014). In the Mf tidal band, the 90° rotation of the tilting direction at HDM and the absence of the LAOs at ECO could be ascribed to this second source of deformation. The delay between fluid accumulation and ground deformation is supported by the observation that threshold effects and mechanisms of mass and energy storage may induce lagged nonlinear responses in many natural systems (Phillips, 2003). The additional fluid emplacement could have generated a mass redistribution, resulting in a kind of dipole formed by the primary and secondary sources, just breaking the radial symmetry of the deformation field.
Our results indicate a change in the dynamics of the Campi Flegrei caldera starting from the second part of the year 2020, likely preceded by a preparatory stage (2018-2019) in agreement with other observations, such as variations in seismicity rate, ground temperatures of the fumarolic field, gravimetric and geochemical parameters. Indeed, an increase in both seismic swarms and background seismicity has been detecting since 2018 (Bellucci Sessa et al., 2021). Remarkable seismic swarms occurred in December 2019 and April 2020, likely promoted by hot fluid injections from depth coupled with cold water from heavy meteoric recharge (Petrosino & De Siena, 2021). Such "double-injection" mechanism produced local thermal effects which, combined with hydraulic pressurization, enhanced potential rock failure and fault instability (Akande et al., 2021). Starting from 2018, a clear endogenous forcing drives the ground temperature variations at Solfatara-Pisciarelli fumarolic field , with an enhancement in December 2019. Furthermore, after the seismic swarm of December 6, a significant rise of CO 2 concentration was observed, superimposing to an already growing trend. Such variation was related to a great increase in the fluid flux at Pisciarelli boiling pool. A significant fluid migration was also detected starting from 2018 by using polarization analysis of seismic noise (Petrosino & De Siena, 2021). In addition, further evidence of fluid mass movements comes from both magnetotelluric imaging (Troiano et al., 2022) and gravity measurements. Indeed, during the last gravimetric surveys of October 2019 (Berrino & Ricciardi, 2020), and September-October 2021 (available at the URL https://www.ov.ingv.it/ index.php/monitoraggio-e-infrastrutture/bollettini-tutti/bollett-mensili-cf/anno-2021-2/1013-bollettino-mensilecampi-flegrei-2021-10/file), a statistically significant gravity increase localized eastward of Pozzuoli of about 30 μgal was revealed. The authors interpret such variations as likely associated with the dynamics of the local hydrothermal system, such as fluid injections at shallow depth.
Although there is no evidence of spatial variation in the primary deformation source inferred from GPS data analysis (De Martino et al., 2021), changes in the deformation rate since 2012 have been recently observed (see, e.g., Tripaldi et al., 2022). This change of the Phlegrean deformation field is also confirmed by the InSAR measurements carried out in the area, compared with the measurements provided by precise leveling lines and cGNSS stations belonging to the Neapolitan Volcanoes (Polcari et al., 2022).
In addition, since 2015 an inflation in the deep deformation source has been resolved by a careful analysis of SAR data (Amoruso & Crescentini, 2022). Our results add a further piece of information from tiltmeter data, whose kinematic investigation combined with seismicity, as well as the dynamical analysis, provide indications of a local second-order deformation source. This derives from the sensitivity of the tiltmeter network in revealing small changes in the deformation pattern.
The integration of the results from NLA approach (applied for the first time to tiltmeter data) with those coming from other geophysical and geochemical observables is valuable and useful for a more realistic modeling of the Campi Flegrei dynamics, which has to take into account first-order mechanism and higher order effects. This FALANGA ET AL.

10.1029/2022EA002702
14 of 16 favors the interpretation of the ongoing dynamical phenomena and opens the route for a prompt detection of variation in the dynamic state, supporting the monitoring activity.
Finally, the new approach based on the application of nonlinear analysis to tiltmeter time series coupled with a kinematical analysis is a powerful tool able to extract the features of a second-order source. Even when the time interval spans a few years, the methodology is effective in evidencing small variations of the dynamical behavior of the system. Due to the generality of the proposed approach, it can be applied to deformation data recorded in other geological contexts (volcanic/tectonic areas) to monitor eventual changes from the background activity or extended to other geophysical and geochemical variables, without losing its efficacy.