Critical Transitions at Campi Flegrei resurgent caldera: A Novel Approach to Detect Early Warning Signals

Critical Transitions at Campi Flegrei resurgent caldera: A Novel Approach to Detect Early Warning Signals | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Critical Transitions at Campi Flegrei resurgent caldera: A Novel Approach to Detect Early Warning Signals Andrea Vitale, Andrea Barone, Enrica Marotta, Dino Franco Vitale, and 14 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7290139/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Understanding how volcanic systems evolve over time is a major challenge due to their complex behavior and constantly changing conditions. This study explores a new way to detect early warning signals of volcanic unrest by analyzing how different types of data, such as ground deformation, gas emissions, temperature, and earthquakes, interact with each other. Focusing on the Solfatara-Pisciarelli volcano system, that is a more active area in Campi Flegrei Caldera (Southern Italy), we used two advanced methods to identify critical transitions in the system: one to model the nonlinear relationships between variables, and the other to detect key moments when the system's behavior shifts. By including time delays between signals (LAG) we found that our model becomes much more accurate in identifying these changes. In contrast, models that ignore time lags show higher uncertainty. These results highlight the importance of considering how different warning signs may not appear all at once, and they offer a promising tool for improving the monitoring of active volcanoes and supporting early risk management. Earth and environmental sciences/Natural hazards Earth and environmental sciences/Solid earth sciences Volcanology Ground deformation Seismicity Gas emissions Critical transitions Early warning Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 1. Introduction Understanding how and when critical transitions occur in volcanic systems is a fundamental challenge in Earth Sciences. These systems are shaped by nonlinear and multivariable interactions between geological, geophysical, and geochemical processes that evolve dynamically over time. Anticipating changes in their behavior requires both continuous multidisciplinary monitoring and analytical frameworks capable of identifying subtle interdependencies and early warning signals. The Solfatara-Pisciarelli (SP) fumarolic area, located within the Campi Flegrei caldera (CFc) in Southern Italy, provides an exceptional natural laboratory to investigate these dynamics (Caliro et al., 2025 ) (Fig. 1 ). The CFc is a large resurgent caldera formed by two major explosive eruptions: the Campanian Ignimbrite (~ 39 ka) and the Neapolitan Yellow Tuff (~ 15 ka) (Steinmann et al., 2016 ; Orsi et al., 2022 ) and has experienced repeated unrest episodes marked by ground deformation, seismicity, and gas anomalies. The SP area, northeast of the town of Pozzuoli, is among the most active sectors of the caldera (Smith et al., 2011 ). It is characterized by intense hydrothermal activity, persistent degassing, and episodic seismic swarms, reflecting the complex coupling between deep magmatic sources and a shallow hydrothermal system (Barone et al., 2025 ). Geophysical and geological studies have revealed that the SP region overlies a diatreme-like structure extending ~ 3 km below the surface, composed of highly fractured volcanic rocks and fault systems that channel deep fluids toward the surface (Isaia et al., 2015 ; Troiano et al., 2022 ; Barone et al, 2025 ). Electrical resistivity tomography, magnetotelluric imaging, and 4D seismic tomography have confirmed the presence of fluid-saturated, high-permeability zones and documented dynamic changes in the caldera’s plumbing system during unrest phases (Siniscalchi et al., 2019 ; Giacomuzzi et al., 2024 ; Troiano et al., 2019 ). These structures play a critical role in controlling fluid migration and pressure redistribution, making the area particularly sensitive to magmatic or hydrothermal perturbations. The interaction between magmatic degassing and hydrothermal circulation also manifests in geochemical variations. Parameters such as CO₂ flux, CO/CO₂ and CO₂/H₂O ratios, and equilibrium temperatures and pressures of the hydrothermal system have shown temporal correlations with deformation and seismicity, offering key insights into the system’s evolution (Chiodini et al., 2003 ; 2012 ; 2015 ; 2022 ). Notably, an increase in fumarolic sulfur emissions since 2018 further suggests a reorganization of fluid pathways and magmatic inputs ( Caliro et al., 2011 ). Gas fluxes reaching tens of thousands of tons per day, and the coexistence of superheated steam and liquid water in the shallow reservoir, underscore the intense energy exchange occurring at depth (Chiodini et al., 2011 ; Bevilacqua et al., 2024 ; Giudicepietro et al., 2025 ). Between 2018 and 2024, a comprehensive monitoring effort has captured high-resolution time series of key parameters, including vertical ground deformation (Tizzani et al., 2024 ; Castaldo et al., 2021 ; Pepe et al., 2019 ; D’Auria et al., 2015 ), seismicity (Tramelli et al., 2024 ; Scotto di Uccio et al., 2024 ), surface temperature from UAV and satellite platforms (Marotta et al., 2024 ; Mercogliano et al., 2024 ), and gas compositions and fluxes. This rich multiparametric dataset enables an unique opportunity to assess how physical and chemical signals co-evolve and potentially anticipate transitions in the system. In this study, we address the need for an integrated methodology to identify and quantify critical transitions in volcanic systems. All analyses performed in the present study were carried out using several datasets acquired from this area of intense volcanic activity, as monitored by the multiparametric surveillance system of the INGV – Osservatorio Vesuviano. We employ a novel framework that integrates Multivariable Fractional Polynomial Analysis (MFPA) (Royston and Sauerbrei, 2008 ) and Global Critical Point Analysis (GCPA) (Rabinowitz, 1986 ). MFPA enables the modeling of nonlinear relationships among variables such as ground deformation, seismicity, heat flow, and CO₂ flux, while incorporating time lags that capture delayed cause-effect interactions. Building on the MFPA outputs, GCPA applies quadratic optimization techniques to identify temporal breakpoints across multiple variables simultaneously, without relying on thresholds or the prediction of specific events. Our aim is to identify systemic changes in the SP hydrothermal-magmatic complex by detecting associations between key parameters and their temporal structure. By uncovering lagged, nonlinear relationships and pinpointing critical transitions, this approach enhances our capacity to interpret volcanic unrest. Moreover, the method is scalable and adaptable to other high-risk volcanic areas, offering a promising tool for improving early warning capabilities and risk mitigation strategies. 2. Material and Methods The dataset analyzed in this study (Table 1 ) spans the period from 2018 to 2024 and includes high-resolution time series of key geophysical and geochemical parameters essential for monitoring the volcanic activity of the Solfatara-Pisciarelli (SP) system. These parameters comprise vertical ground deformation (InSAR), seismicity (from the SERENADE-GOSSIP catalog), fumarolic CO₂ flux, and thermal data acquired through UAV-based infrared imaging. Additional geochemical parameters—such as CO concentration, CO/CO₂ and CO₂/H₂O ratios, as well as equilibrium pressure and temperature—were obtained from a combination of automatic monitoring stations and monthly laboratory analyses conducted by INGV–Osservatorio Vesuviano. The monitoring focuses on the main active degassing vents: Bocca Grande (BG) and Bocca Nuova (BN) at Solfatara, and the Pisciarelli mud pool fumaroles, located less than 300 meters away on the external eastern flank of the Solfatara crater (Fig. 1 B). Due to the close spatial proximity and the consistent geochemical and physical behavior observed across these sites—such as similar gas compositions, temperature trends, and fumarolic dynamics—these areas are treated as a unified hydrothermal-magmatic system. This interpretation is supported by previous studies (Siniscalchi et al., 2019 ; Cardellini et al., 2017 ; Petrillo et al., 2013 ), which highlight the common origin and interconnected behavior of the two fumarolic fields. Accordingly, in this study, the Solfatara and Pisciarelli fields are collectively referred to as the SP system. To improve the spatial characterization of the seismicity and its temporal evolution, a clustering analysis was applied to the SERENADE-GOSSIP earthquake catalog using the DBSCAN algorithm (Ester et al., 1996 ; Cirillo et al., 2022 ). This density-based method identifies spatially coherent groups of seismic events and allows for the discrimination between background activity and structurally significant clusters. The analysis revealed a dominant vertical structure (Cluster 1), located beneath the Solfatara-Pisciarelli area, which is interpreted as a persistent fluid and/or fracture conduit related to the active hydrothermal system. Additional minor clusters were detected across the wider caldera, associated with transient or localized swarm-like activity. The identification of Cluster 1 was instrumental in isolating the subset of events most relevant to the SP system and in enhancing the temporal resolution of the multiparametric analysis. All time series were homogenized and normalized to facilitate cross-comparison and integration into the multivariate statistical models presented in the following sections. This section is organized in two parts: the first describes the sampling procedures and acquisition protocols for each dataset—including sensor specifications, thermal clustering methods, DBSCAN parameters, and data preprocessing steps—while the second presents the resulting time series for each parameter. Additional technical details, including spatial cross-sections of each seismic cluster, are provided in the Supplementary Materials (Figs. S1–S6). 2.1 Data collection Fumarolic Gas Composition (Bocca Grande - BG, Solfatara) The compositional data used in this study were collected monthly from Bocca Grande (BG), the main fumarolic vent at La Solfatara, using an automatic sampling station. Gas samples were analyzed for chemical composition at the Geochemistry Laboratory of INGV-Osservatorio Vesuviano. Samples were collected in under-vacuum flasks containing 4N NaOH solution (Giggenbach, 1975 ; Giggenbach and Gouguel, 1989 ) and analyzed for major gas species. Analytical procedures followed Cioni and Corazza ( 1981 ) , using gas chromatography with a single injection on dual molecular sieve columns (MS 5 Å capillary, 30 m × 0.53 mm × 50 µm) with TCD detectors and He/Ar as carrier gases. CO₂ absorbed in the alkaline solution was oxidized with H₂O₂ and quantified by acid–base titration. CO was measured in dry gas samples using a water-cooled condenser (20–30°C), chromatographic separation on a MS 5 Å 1/8 × 50 in column (He carrier), and detection with a high-sensitivity Reduced Gas Detector (HgO). CO₂ Flux and Thermal Anomalies (Pisciarelli Site) CO₂ emissions at the Pisciarelli fumarolic field were continuously monitored from 1 December 2018 to 30 May 2024, with daily resolution and point-based spatial data. Thermal surveys were carried out monthly from September 2019 to October 2023, with additional surveys up to May 2024. UAV-based radiometric thermal imaging was conducted during nighttime or evening hours to minimize solar interference. About 50 thermal mosaics were generated from flights at altitudes between 55 and 70 m a.g.l. (Marotta et al., 2024 ), with some mosaics combining data from different flight plans. Thermal mosaics were processed using the method of Marotta et al. ( 2019 ) to delineate surface thermal anomalies and quantify emitted thermal energy. Clustering algorithms (DBSCAN and k-means; Ester et al., 1996 ; Cirillo et al., 2022 ) were applied to segment images into thermally homogeneous areas. Each cluster was annotated with temporal metadata and thermal attributes (e.g., average temperature, thermal flux, size in pixels and m²), then reanalyzed with DBSCAN to track temporal evolution. Full details of data acquisition and processing are reported in Marotta et al. ( 2024 ) Seismicity (GOSSIP/SERENADE Catalog) Seismic data were sourced from the SERENADE database (Peluso, 2018a ; 2020 ) via the GOSSIP portal (Ricciolini et al., 2024), which provides hypocentral parameters, magnitude (Md), and station arrival times. Hypocenter locations were calculated using the Hypo71 code (Lee and Lahr, 1975 ) with local 1D velocity models for Vesuvius, Campi Flegrei, Ischia, and surrounding regions. The database integrates automatic Earthworm-based solutions (Johnson et al., 1995 ) and manual revisions by the INGV-OV Monitoring Room and Seismic Laboratory. Static event pages are generated with the Serewrap software (Di Filippo and Peluso, 2019 ), triggered with each new event insertion via the WESSEL portal (Peluso, 2018b ; 2020 ). Seismic events from 2018 to May 2024 with Md ≥ 0.2 were extracted and spatially filtered within a 1000 m buffer around the target area. Ground Deformation (Sentinel-1 InSAR Analysis) Ground deformation data were derived from Sentinel-1A (S1-A) C-band SAR images (λ = 5.54 cm) acquired in Interferometric Wide mode (IW) from January 2018 to April 2024, along ascending (Path 44) and descending (Path 22) orbits, each with 191 images. The Small BAseline Subset (SBAS) multi-temporal DInSAR method (Berardino et al., 2002 ; Jolivet et al., 2011 ; Rosen et al., 2000 ) was applied independently to both tracks. A total of 2101 interferograms were generated per track, with a maximum temporal baseline of 144 days. Interferograms were flattened using precise orbits and the NASA SRTM 1-arc-second DEM ( https://dwtkns.com/srtm30m ). Ascending and descending LOS displacements were combined (Manzo et al., 2006 ; Pepe et al., 2016 ) to obtain vertical and horizontal deformation maps, geocoded and referenced to a stable point. The final dataset achieved cm-to-mm accuracy (Tizzani et al., 2024 ), with a spatial resolution of ~ 90 m and a ~ 12-day temporal resolution, totaling 191 layers. Mean vertical deformation was computed from the time series of four pixels within 100 m of the target area. Table 1 A summary table of the considered dataset is provided. The observation time interval and frequency of each individual measurement are reported. Data Time-Interval Period Ground Deformation 2018–2024 12 Days Seismic Data 2018–2024 Continous Thermal Survey 2019–2024 Montly CO 2 Emission 2019–2024 Daily Fumarolic Gas 2016–2024 Montly 2.2 Data description The analysis of geochemical, thermal, seismic, and ground deformation data collected between 2018 and 2024 at the Solfatara-Pisciarelli system highlights a consistent multi-parameter response to evolving subsurface dynamics. The following sections briefly describe and interpret the observed trends. Temporal trends in fumarolic gas composition The SP hydrothermal system exhibits complex geochemical dynamics. The time series analysis of various geochemical parameters provided valuable insights into the system's temporal variations and underlying processes. Several key ratios were monitored, including T(CO-CO₂) [°C] which represents the equilibrium temperature for the hydrothermal system, CO [µmol/mol], CO/CO₂, CO₂/H₂O, CO₂/CH₄ ratio, and P(CO-CO₂-H₂-H₂O) [bar], which indicate the equilibrium pressure. These parameters provided insight into gas composition shifts, which are indicative of changes within the magmatic system (Fig. 2 ). CO measurements (Fig. 2 A): the compositional data used in this work were collected from Bocca Grande (BG), the main fumarolic vent of La Solfatara, using an automatic station. From 2018 to the beginning of May 2020, there was a notable increase in CO concentration. Post-May 2020, the concentration decreases rapidly until stabilizing around 2022. After 2022, the CO concentration exhibited limited fluctuations around a stable mean value. Post-2022 oscillations appear periodic and regular, while fluctuations between 2018 and 2021 were more irregular and pronounced. These observations align with studies indicating that high CO₂/CH₄ ratio in fumaroles can signal magmatic components influencing the hydrothermal system. Equilibrium Temperature of the CO-CO₂ System (Fig. 2 B): from 2018 to May 2020, there was a significant increase in equilibrium temperature. Post- May 2020, the temperature decreased rapidly, stabilizing around 2022. After 2022, temperatures exhibited limited fluctuations around a stable mean. Post-2022 oscillations appear periodic and regular, whereas fluctuations between 2018 and 2021 were less regular and more pronounced. The peak observed around 2021 likely represents a critical juncture for the system. Variations in equilibrium temperatures suggest underlying chemical or physical dynamics. CO/CO₂ Ratio (Fig. 2 C): from 2018 to May 2020, there was a significant increase in the CO/CO₂ ratio. Post-2021, the ratio decreased rapidly, stabilizing around 2022. After 2022, the CO/CO₂ ratio exhibited limited fluctuations around a stable mean, with oscillations appearing periodic and regular. Fluctuations between 2018 and 2021 were less regular and more pronounced. The peak around June 2020 likely represents a critical point for the system. Variations in the CO/CO₂ ratio suggest underlying geochemical and geophysical dynamics. Equilibrium Pressure of the Hydrothermal System P(CO-CO 2 -H 2 -H 2 O) (Fig. 2 D): from 2018 to 2021, there was a notable increase in equilibrium pressure. Post-2021, the pressure decreased rapidly, stabilizing around 2022. After 2022, equilibrium pressure showed limited fluctuations around a stable mean, with oscillations appearing periodic and regular. Fluctuations between 2018 and 2021 were less regular and more pronounced. The peak around 2021 likely signifies a critical point for the system. Variations in equilibrium pressure suggest underlying geochemical and geophysical dynamics. CO₂/H₂O Ratio (Fig. 2 E): from 2018 to 2021, there was a significant increase in the CO₂/H₂O ratio. Post-2021, the ratio decreased rapidly, stabilizing around 2022. After 2022, the CO₂/H₂O ratio exhibited limited fluctuations around a stable mean, with oscillations appearing periodic and regular. Fluctuations between 2018 and 2021 were less regular and more pronounced. The peak around 2021 likely represents a critical point for the system. Variations in the CO₂/H₂O ratio suggest underlying geochemical and geophysical dynamics. CO₂/CH₄ Ratio (Fig. 2 F): from 2018 to April 2020, there was a significant increase in the CO₂/CH₄ ratio. Post-2021, the ratio increased rapidly, stabilizing around 2022. After 2022, the CO₂/CH₄ ratio exhibited limited fluctuations around a stable mean, with oscillations appearing periodic and regular. Variations in CO2 fluxes at Pisciarelli and their implications CO₂ Flux measurements (Fig. 2 G): the Pisciarelli fumarolic site has been continuously monitored for CO₂ emissions over a data interval spanning from 1 December 2018 to 30 May 2024. Figure 2 G presents the time series analyzed as a monthly moving average, while the semiannual average (blue curve) highlights long-term trends and patterns. The time series reveals multiple peaks, most notably around mid-2020, late 2021, and early 2023, where the CO₂ flux exceeded 3.5 × 10⁴ g/m²/day. These peaks coincide with known periods of increased volcanic activity in the CFc. Significant reductions in CO₂ emissions were observed, with the lowest flow values occurring around late 2020 and late 2022, suggesting potential periods of reduced degassing or transient system stability. The semiannual average smooths out daily fluctuations revealing recurring cyclic patterns which may indicate a potential link between seasonal changes and degassing activity. Notably, periods of heightened activity are interspersed with phases of reduced emissions, underscoring the cyclic nature of the hydrothermal system's behavior. From 2020 onward, the amplitude of fluctuations increased, reflecting a potential intensification of subsurface processes, such as magmatic degassing or variations in hydrothermal reservoir pressurization. Over the monitoring period, the CO₂ flux shows a gradual increase in baseline values, which may signal ongoing pressurization within the magmatic-hydrothermal system or evolving permeability pathways. The detailed analysis of the CO₂ flux at Pisciarelli highlights both short-term variability and long-term trends. Peaks in emissions align with known episodes of unrest, while the semiannual average reveals broader patterns that provide crucial insights into the behavior of the volcanic system. In summary, the SP hydrothermal system has exhibited significant variations across multiple geochemical parameters between 2018 and 2024. The period from 2018 to 2021 is characterized by increasing trends, reaching peaks around 2021, followed by rapid declines and stabilization around 2022. Post-2022, the system shows regular, periodic fluctuations, suggesting the establishment of a new equilibrium state may be influenced by cyclic external factors. Spatial and temporal evolution of surface thermal anomalies The thermal data used in this study extend beyond those presented in Marotta et al. ( 2024 ) . Figure 3 integrates spatial and temporal analyses, illustrating the heating dynamics through land surface temperature (LST) and statistical trends. Specifically, Fig. 3 A shows mean and annual thermal maps from 2020 to 2024. In Fig. 3 B, the LST time-series plot illustrates measured temperatures (black), de-seasonalized temperatures (red), and seasonal trends (green) are reported. The data reveal a pronounced seasonal cycle (green) superimposed on broader temperature variations. De-seasonalized temperatures (red) isolate long-term thermal anomalies, providing insights into subsurface processes. Figure 3 C focuses on de-seasonalized temperatures (red) alongside their six-month moving average (black). The smoothed trendline offers a clearer visualization of gradual temperature changes, indicating potential long-term geothermal or hydrothermal dynamics in the study area. Figure 3 D shows the heat flux time series (blue) and the heat release rates (red). The consistent hotspots in mean map (red/yellow zones in Fig. 3 A) suggest persistent thermal anomaly activity, providing a comprehensive view of the thermal behavior over the analyzed period. Specifically, the mean annual map highlights two dominant patterns of thermal activity within the studied area. The first pattern, located on the right side of the maps, corresponds to the Pisciarelli mud pool area, exhibiting the highest thermal anomalies, with intense red zones persisting across all years, likely due to concentrated hydrothermal venting and localized geothermal heat flux. The second notable anomaly pattern is observed along the external flank of the Solfatara volcano, specifically near Via Scarfoglio area on the left side of the maps. This region consistently displays medium-to-low thermal anomalies (green-to-yellow zones), indicative of diffuse heat flow and potential subsurface thermal activity extending from the Solfatara hydrothermal system. Summarizing the observations: (i) Pisciarelli anomaly is a persistent and intense hotspot with slight interannual variability in intensity and spatial extent, representing the core of hydrothermal venting and geothermal heat emission. (ii) Via Scarfoglio Pattern consists of a diffuse and moderate heat flow along the external flank of the Solfatara volcano, showing gradual spatial expansion over time, potentially linked to subsurface dynamics. (iii) The temporal evolution of the LST field emphasizes that year 2022 is remarkable for an increased thermal activity, particularly in the Pisciarelli zone (Fig. 3 C) even if the overall spatial distribution of thermal anomalies remains consistent. Accordingly, a computed heat flow time series derived from the UAV-based thermal dataset from September 2019 to October 2023 is reported in Fig. 3 D, with a variable temporal resolution (minimum value 1 month). The estimate considers between 5 to 20 points per temporal instance, with a six-month moving average applied. Analysis of the heat flux and its variations: Fig. 3 D presents the time series of heat flux (blue) and heat flux variations (red), providing insights into the geothermal dynamics of the system. The heat flux series exhibits a relatively smooth trend interspersed with significant fluctuations. From 2020 to early 2022, a gradual increase is observed, indicating a progressive intensification of geothermal activity. This trend stabilizes around mid-2022, followed by a slight decline extending through late 2023. The six-month moving average highlights long-term patterns by filtering short-term noise and emphasizing overall heat flux trajectory in the region. A comparative analysis between heat flux and its variations (Fig. 3 D) reveals important correlations. The most prominent peak in mid-2022, is clearly evident in both heat flux (blue) and its variations (red), suggesting a system-wide thermal perturbation. However, while the heat flux series displays a smoother evolution, the variations series exhibits pronounced spikes, indicative of episodic fluctuations in the geothermal regime. Heat flux variations, represented on the right axis of Fig. 3 D, highlight relative intensity changes over time. The series of variations maintains a generally oscillatory pattern, with fluctuations contained within a predictable range. This consistency suggesting a dominant seasonal component likely influenced by atmospheric conditions and external factors affecting surface thermal measurements. In contrast, significant spikes in heat flux, particularly in mid-2020 and mid-2022, indicate abrupt events potentially linked to fluid injections, pressure changes, or localized thermal surges. These events temporarily disrupt the otherwise stable trajectory of the system. The synchronization of anomalies in both series, particularly around mid-2022, underscores a system-wide geothermal perturbation, likely driven by subsurface dynamics. By examining heat flux series as the central trend and analyzing deviations through the variations series, it becomes evident that the Pisciarelli geothermal system exhibits both stable long-term trends and episodic fluctuations. Seismic activity and clustering beneath the hydrothermal area In the present study, the CF seismic catalog was analyzed using the DBSCAN algorithm to identify and characterize spatial clusters of seismic activity (Supplementary Text). This method, particularly suited for datasets with varying density distributions, enables the detection of regions where seismic events are densely concentrated while distinguishing isolated events considered as noise. The results of this analysis provide a refined picture of the seismic dynamics within the caldera, offering critical insights into the spatial organization of earthquakes and their potential connection to underlying geological structures. The outcomes of the clustering analysis are illustrated in Fig. 4 (see also Figs. S1-S6), where the spatial distribution of seismic clusters is represented through contour lines depicting event density. The plan view (Fig. 4 A) offers an overview of how seismicity is distributed across the Campi Flegrei area, with different clusters identified by distinct colors. The contour lines effectively highlight the areas with the highest concentration of earthquakes, delineating regions of intensified seismic activity. To further investigate the depth distribution of these clusters, two vertical cross-sections were considered. The North-South cross-section (Fig. 4 B) provides a view of how seismicity extends in depth along this axis, revealing well-defined zones of high activity. Similarly, the East-West cross-section (Fig. 4 C) confirms the presence of deep-seated clusters and highlights their spatial correlation with shallow seismicity, suggesting a structured arrangement of events that may be influenced by fault networks and magmatic or hydrothermal processes. Among the various clusters identified, particular attention is given to Cluster 1, located directly beneath the SP area. The seismicity associated with this cluster presents a highly concentrated pattern, both laterally and in depth, and is likely influenced by subsurface fluid movements and stress accumulation within the volcanic system. Given its strategic location and relevance, this study will focus specifically on the seismic events belonging to Cluster 1, aiming to investigate their temporal evolution and potential implications for the dynamic behavior of the SP hydrothermal system. The spatial clustering results obtained through DBSCAN not only confirm the presence of localized zones of high seismicity but also provide essential clues regarding the interaction between tectonic and magmatic processes within Campi Flegrei. The density contours reveal that earthquake distribution is not random but follows specific structural trends, reinforcing the hypothesis that seismicity in the area is closely linked to the caldera’s ongoing geodynamic evolution. The identification and analysis of these clusters contribute to a deeper understanding of the mechanisms driving volcanic unrest, such as fluid migration and/or fault activities. Figure 4 D provides a representative time series of earthquakes per month for Cluster 1 (red in Fig. 4 A-C) over the period from 2018 to 2024, with normalized density values on the y -axis and time in years on the x -axis. From 2018 to the end of 2020, seismicity density remained relatively low and stable, indicating minimal seismic activity in the region. A gradual increase became evident starting in 2020, suggesting a progressive intensification of the underlying geological processes. This trend continues, culminating in a significant rise between 2022 and 2023, marking a period of pronounced seismic activity. In 2023, the seismicity density reached its peak, followed by fluctuations in early 2024. These variations may be attributed to changes in underground pressure, fluid migration, or small-scale eruptive dynamics. Toward the latest part of the dataset in 2024, a slight decline in seismicity density is observed, potentially indicating a temporary stabilization or a reduction in subsurface activity. Overall, the data highlight a phase of increased seismic activity in SP, which may be linked to ongoing volcanic or geothermal processes. A more detailed analysis integrating geophysical and geochemical parameters would be necessary to further elucidate the mechanisms driving these observed seismic trends. Ground deformation as an indicator of subsurface processes The vertical ground deformation at the Pisciarelli area, derived from InSAR data and based on the average of three coherent SAR pixels within the region of interest, exhibits notable spatial and temporal characteristics. Spatially, the mean ground velocity map (Fig. 5A) highlights a well-defined uplift area, with maximum velocities exceeding 12 cm/year. The deformation pattern decreases radially from the central uplift, forming a gradient consistent with subsurface processes such as hydrothermal activity. The vertical component was then used to compute the monthly deformation rate, which was smoothed using a six-month moving average. Temporally, the cumulative deformation time series (Fig. 5B) shows a continuous and steady increase from 2018 to 2024, reaching a total displacement of approximately 45 cm by the end of the observation period (for the pixels in Pisciarelli area). The near-linear trend in cumulative deformation indicates persistent driving forces in the system, likely related to pressurization or volumetric expansion within the subsurface structure. The time series of vertical deformation rates (Fig. 5C) provides a more detailed view of the temporal variability. The rates exhibit oscillatory behavior, with periods of acceleration and deceleration. Peaks in deformation rate occur around 2019, 2021, and late 2023, suggesting episodic intensifications in the processes driving deformation. The semiannual moving average captures overall trends while smoothing short-term fluctuations, showing a general increase in deformation activity punctuated by temporary slowdowns. The dataset, collected with high temporal resolution (~ 12 days) and spatial resolution (~ 90 m), effectively captures both the localized nature and dynamic evolution of the deformation. The observed patterns suggest a strong connection between subsurface processes and surface deformation, emphasizing the importance of continuous InSAR monitoring for understanding the ongoing dynamics of the Pisciarelli area. 2.3 Data Processing and methods The above-mentioned dataset were initially processed to allow comparisons across varying units and scales. In particular, normalization and temporal smoothing techniques (e.g., moving averages), have been applied to highlight significant trends. 2.3.1 Multivariable Fractional Polynomial Analysis (MFPA) The primary objective of this study was to investigate the association between the temporal behavior of deformation, seismic activity, and specific physical and chemical factors. This goal, along with the characteristics of the available dataset, guided the selection of specific statistical procedures and analytical approaches. A multivariable regression model was employed to assess the association between our main variable of interest (vertical ground deformation) and the set of measured physical parameters. The aim is to identify statistically significant correlation for each variable, independent of other factors. Since the data are structured as a single-point time series across multiple variables, an autoregressive model was adopted to mitigate statistical bias arising from the dependent variable’s non-independence, which can occur in regression models for this type of data. To address this, LAG value of the dependent factor was included as an additional independent variable in the model. The time lag order, i.e. the number of periods going back in time, was set to 1 unit (i.e., 1/10 of a year) based on the following considerations: (i) while increasing the LAG order may improve forecasting accuracy, it also increases model complexity, reduces interpretability, and raises the risk of multicollinearity. Since the goal of this work is not to develop a forecasting model, we keep this parameter low. (ii) A first-LAG model sufficiently captures the temporal dependence of the response on the dependent variable, as each value directly influences the subsequent one (Finkel, 2007 ), thereby addressing the temporal dependence bias. Additionally, by acting as a proxy for unknown factors that affect the dependent variable, the first LAG absorbs unmeasured effects, helping to stabilize estimates of the independent effects of other covariates and enabling more accurate evaluation. To assess the coherence of the relationships identified by each model, comparisons were made between models with and without the autoregressive (LAG) variable. The identification of the best final model followed the pragmatic approach suggested by Royston and Sauerbrei ( 2008 ) , summarized as follows: Selection of the final model by including factors significantly associated with the dependent variable, while assessing the functional form (linearity/non-linearity) of continuous factors using the multivariable fractional polynomial (MFP) algorithm. Evaluation of the overall and individual contribution of the factors selected in the final model to the dependent variable by calculating both the global variance explained (R²) and the individual variance partition for each significant factor (R²p). The Shapley-Owen decomposition algorithm (Shorrocks, 2013 ) was used to partition global R², quantifying each variable's impact independently of others, i.e., the impact with all other factors held constant. All models were checked for multicollinearity by calculating the Variance Inflation Factor (VIF) for each variable in the final model. A VIF value greater than 5 was considered indicative of multicollinearity (Kutner et al., 2005 ; Suleiman, 2015 ). To visualize the relationship between a dependent variable and a specific independent factor controlling for other significant variables in the final model, a scatterplot was created. In this plot, the dependent variable values were adjusted with respect to the mean of other significant variables considered in the model. Specifically, the relationship between y and one of the x variables was illustrated by adjusting observed y values for each item according to deviation of other observed x variables from their respective means, using estimated coefficients. Data were analyzed using Stata version 18.0 (StataCorp LP, College Station, TX, USA). Statistical significance was set at p ≤ 0.05 for univariate comparisons, variable selection, and testing between variable transformations within the MFP procedure. The factors (i.e., physical parameters) significantly associated in the performed model are used as input for a Global Critical Point Analysis to identify meaningful points in the time domain following described. The MFPA model converged in fewer than 30 iterations per run, with an average computation time below 1 minute on standard desktop hardware (Intel i7 processor, 16GB RAM). This performance confirms the suitability of the approach for potential near-real-time implementations in operational monitoring environments. 2.3.2 Global Critical Point Analysis (GCPA) GCPA method is designed to identify significant temporal transitions in multiparametric historical datasets. It leverages quadratic programming (QP) ( 57 ) to optimize the detection of a critical point, that represents the time when the combined dynamics of multiple variables exhibit significant changes. This critical point may correspond to a phase where the system transitions from stability to chaos or shifts between different equilibrium states. In the context of the SP system, GCPA serves as a valuable tool for analyzing variations in geophysical and geochemical parameters, helping to detect potential systemic changes without aiming to predict specific events, such as an eruption. Instead, the methodology focuses on understanding the system's evolution and identifying precursors to critical transitions, essential in volcanic systems where gradual or subtle changes may indicate an approaching tipping point (global critical point). The GCPA statistical technique identifies the temporal moment when significant changes occur across multiple variables simultaneously. These variables include physical or chemical parameters recorded over time, such as ground deformation, seismicity, temperature, heat flow, and CO₂ flux. By pinpointing a single temporal moment representing overall change across all analyzed time series, GCPA simplifies the complexity of multiparametric analysis and identifies a common breaking point explaining observed variations. Moreover, the application of quadratic optimization using the Karmarkar method (Karmarkar, 1984 ) has improved the efficiency of global critical point detection, allowing for a more precise assessment of systemic changes. This optimization framework enhances the reliability of transition detection by providing a mathematical basis for identifying system shifts beyond simple threshold-based analyses. The analysis was performed using R software (R Core Team, 2014 ) and involved five main steps: Data Preprocessing and Normalization: input time series are normalized to standardize datasets with different units or scales, (e.g., temperature, seismicity, gas flux), ensuring meaningful and robust comparison of variables while reducing the risk of bias in subsequent analysis. Quadratic Programming for Critical Point Identification: the problem is modeled using a quadratic objective function that maximizes deviations from expected or average behavior. The identity matrix Q represents the quadratic structure of the optimization problem, while the vector c captures normalized absolute deviations. Temporal constraints are applied to ensure that the identified critical point lies within the observed data range, preventing the selection of spurious points. To formally identify the time of critical transition, we structured the analysis as a quadratic programming (QP) optimization problem. Specifically, the critical point was defined as the time index minimizing the following objective function: $$\:\left(\raisebox{1ex}{$1$}\!\left/\:\!\raisebox{-1ex}{$2$}\right.\right)\bullet\:{\varvec{x}}^{\varvec{T}}\bullet\:\varvec{Q}\bullet\:\varvec{x}+{\varvec{c}}^{\varvec{T}}\bullet\:\varvec{x}$$ which represents a quadratic function in matrix form, commonly used in quadratic optimization (Quadratic Programming, QP); where \(\:\varvec{x}\in\:{\varvec{R}}^{\varvec{n}}\) is a binary indicator vector (with 1 at the candidate critical time and 0 elsewhere), \(\:\varvec{Q}\in\:{\varvec{R}}^{\varvec{n}\times\:\varvec{n}}\) is a symmetric matrix representing the covariance structure between the standardized input variables (often also defined as positive semi-definite, to ensure convexity), and \(\:\varvec{c}\in\:\varvec{R}\) is a column vector of linear coefficients collecting the normalized mean deviation of each parameter over time. More specifically, the function: $$\:\varvec{f}\left(\varvec{x}\right)=\left(\frac{1}{2}\right)\bullet\:{\varvec{x}}^{\varvec{T}}\bullet\:\varvec{Q}\bullet\:\varvec{x}+{\varvec{c}}^{\varvec{T}}\bullet\:\varvec{x}$$ is a typical objective function in quadratic optimization problems, where: \(\:\left(\frac{1}{2}\right)\bullet\:{\varvec{x}}^{\varvec{T}}\bullet\:\varvec{Q}\bullet\:\varvec{x}\) , is the quadratic term, and represents a paraboloid if QQ is positive definite; $$\:{\varvec{c}}^{\varvec{T}}\bullet\:\varvec{x},\:\text{i}\text{s}\:\text{t}\text{h}\text{e}\:\text{l}\text{i}\text{n}\text{e}\text{a}\text{r}\:\text{t}\text{e}\text{r}\text{m}.$$ Note that the factor 1/2 is conventional and useful because the derivative of the quadratic term is \(\:\varvec{Q}\bullet\:\varvec{x}\) which simplifies the gradient computation: $$\:\nabla\:\varvec{f}\left(\varvec{x}\right)=\varvec{Q}\bullet\:\varvec{x}+\varvec{c}$$ Finally, we remark that the matrix Q encodes the co-evolution and coupling between variables (e.g., deformation, seismicity, gas), while c provides a temporal profile of their joint anomalies. This formulation allows us to detect the time step where the system collectively exhibits the greatest departure from background behavior. Optimization Techniques: the R Optimization Infrastructure (ROI) package (Theußl et al., 2020 ), combined with the ROI.plugin.quadprog plugin , is used to solve the quadratic programming problem. This setup efficiently handles large datasets with multiple parameters and integrates advanced optimization techniques, including: Interior-Point Methods : ensure convergence by iteratively refining the solution within the feasible region. Conjugate Gradient Methods : particularly suitable for large and sparse datasets. Penalty and Barrier Methods : facilitate constraint handling and enhance computational efficiency. Critical Point Detection: the index corresponding to the maximum deviation within the normalized time series is identified as the global critical point, representing the moment when significant transitions occur across the analyzed parameters (Boyd and Vandenberghe, 2004 ). Visualization and Interpretation: the identified critical point is visualized within the context of normalized and smoothed time series using the ggplot2 package (Wickham, 2009 ) for trend enhancement. Additional visualizations compare raw, normalized, and smoothed data to provide a comprehensive view of the dataset's temporal evolution (Shumway and Stoffer, 2017 ). The GCPA method is particularly suited to complex systems with dynamic and interdependent parameters. By integrating significant geophysical and geochemical variables, such as ground deformation, seismicity, temperature, heat flow, and CO₂ flux, the algorithm identifies moments of significant systemic change. These insights are crucial for understanding volcanic behavior, assessing hazards, and refining monitoring strategies. The GCPA framework is adaptable to large datasets and diverse scientific applications, allowing the incorporation of various constraints and objective functions tailored to specific research needs. For our analysis, we use the time-series of physical parameters retrieved by the MFPA. To evaluate the method's response and its robustness to variations in multiparametric time series, we performed two analyses by considering different observation intervals for the parameters: the first (2018–2023) that represents a subset of the second (2018–2024). 3. Results 3.1 Best-fitting multivariable model for ground deformation using time-lagged analysis We applied MFPA to determine the best combination of geophysical and geochemical parameters that explain ground deformation at the SP site, considered as the dependent variable. This choice is motivated by the prominent role that vertical ground deformation plays in signaling pressurization processes at the Campi Flegrei caldera. Long-term geodetic observations have consistently linked episodes of uplift to periods of increased hydrothermal or magmatic activity (Del Gaudio et al., 2010 ; Tizzani et al., 2024 ; Vanorio et al., 2025 ). In the Solfatara-Pisciarelli area, deformation is often the first measurable indicator of subsurface changes, and it integrates complex interactions such as gas accumulation, heat transfer, and fluid migration. Therefore, it provides an effective dependent variable for exploring multiparametric cause-effect relationships in the SP system. In addition, deformation is measured with high spatial and temporal resolution via InSAR techniques, ensuring the reliability and continuity required for regression-based analyses. This makes it a valuable proxy for unrest modeling, particularly when combined with complementary parameters such as CO₂ flux, seismicity, and heat flow. The goal was to capture significant nonlinear associations and delayed effects that influence vertical deformation dynamics, offering a robust interpretation of the physical mechanisms behind observed surface displacements. Following the assessement of multicollinearity among all parameters (geochemical parameters: T(CO-CO₂) [°C], CO [µmol/mol], CO/CO₂, CO₂/H₂O, CO₂/CH₄, P(CO-CO₂-H₂-H₂O) [bar] and CO₂ flux; land surface temperature; heat flux and seismicity), only those showing no collinearity were included in the MFPA. Two models were compared: one that incorporates a 30-day time lag (LAG model), and one that does not (NO-LAG model). To enhance the interpretability of temporal dependencies among monitored parameters, we developed two regression models: one incorporating a time-lag structure (LAG model) and one without (NO-LAG model). In the LAG model, the dependent variable (vertical ground deformation) is regressed against its own past value (i.e., the value observed 30 days earlier), in addition to other explanatory variables. This autoregressive term captures delayed cause-effect dynamics that are common in volcanic systems, where processes such as pressurization, degassing, and heat transfer propagate over time through subsurface pathways. The rationale behind the 30-day lag was based on: the average response time of deformation to known geochemical variations at CFc as inferred from empirical monitoring (Chiodini et al., 2017 ); the need to introduce temporal structure into the model without overfitting or increasing multicollinearity; the goal of capturing the immediate memory of the system while ensuring robustness. This approach aligns with previous studies modeling dynamic interactions in volcanic environments, where LAGs in the order of weeks to months are commonly used to represent internal feedbacks and delayed surface responses. From a statistical perspective, the first-LAG structure offers a balance between capturing temporal autocorrelation and avoiding overfitting or collinearity issues associated with higher-order lags. Testing multiple lag values would increase model complexity and reduce interpretability, without providing substantial gain in explained variance for the current scope of identifying critical transitions (rather than forecasting). Nonetheless, future developments will explore multi-lag configurations and rolling-window approaches to further generalize the framework. Table 2 summarizes the coefficients, explained variance, and model selection metrics for both approaches. In the LAG model (Panel A), we found that four predictors significantly contribute to deformation: the first lag of deformation, seismic activity, heat flow, and CO₂ flux. This model explained 99.9% of the total variance (global R² = 0.999), demonstrating that it provides the optimal multivariable configuration for this dataset. The first-lag term alone explained 59% of the variance, followed by seismic activity (33.1%), heat flow (4.9%), and CO₂ flux (3.0%). The Akaike Information Criterion (AIC = -12.3) and Bayesian Information Criterion (BIC = -4.9) further supported the model’s. Panel A Panel B Global R 2 = 0.999 Global R 2 = 0.97 AIC= -12.3; BIC=-4.9 AIC = 114.9; BIC = 123.8 Parameter Funct. Form Coeff. p R 2 p (%) BIF (%) Funct. Form Coeff. p R 2 p (%) BIF (%) First Lag Lin. 0,98 ≤ 0.001 59 100 - - - - - Seismicity Lin. 9,3 ≤ 0.001 33,1 97,4 1/x 2 -0,02 ≤ 0.001 41,8 99,7 Temperature Lin. 0,005 = 0.58 - 19,2 Lin. -0,2 ≤ 0.001 1,0 84,7 Heat Flow Lin. 0,00023 = 0.003 4,9 24,4 x 2 2,7 ≤ 0.001 12,7 99,7 CO 2 Flow x 3 -0,01 ≤ 0.001 3,0 65,9 Lin. 0,0002 ≤ 0.001 9,8 97,0 P(CO-CO 2 -H 2 -H 2 O) Lin. -0,03 = 0.36 - 67,9 Lin. -0,7 ≤ 0.001 34,7 99,5 CO 2 -H 2 O Lin. -9,5 = 0.15 - 24,1 Lin. -27,2 = 0.44 - 18,1 Constant - 18 ≤ 0.001 - - - 18,8 ≤ 0.001 - - Table 2 . Deformation regression models with (PANEL A) and without (PANEL B) LAG autoregressive term. Bold face marks the significant factors of each final model. AIC = Akaike information criterion. BIC = Bayesian Information Criterion. R 2 p = Fraction of explained variance. Funct. Form = Functional form. Coeff. = Coefficient. Lin. = Linear. p = observed alfa error. BIF = Bootstrap inclusion frequency. Figure 6 A displays the observed deformation time series (blue dots) overlaid with the LAG model prediction (red line), showing excellent agreement. The residuals (green dots) are minimal and randomly distributed. In contrast, the NO-LAG model (Fig. 6 B) resulted in increased residuals, lower explanatory power (R² = 0.97), and significantly higher AIC (114.9) and BIC (123.8) values. Covariate-adjusted partial relationships (Fig. 7 ) revealed a strong positive association for the first-lag deformation (Fig. 7 A) and seismic activity (Fig. 7 B). Heat flow (Fig. 7 C) showed a negative linear relationship, while CO₂ flux (Fig. 7 D) displayed a nonlinear response. The bootstrap inclusion frequency (BIF) analysis confirmed the LAG model's robustness: the LAG term and seismicity exceeded a 95% inclusion rate, while CO₂ flux and P(CO–CO₂–H₂–H₂O) surpassed the 50% threshold. Conversely, parameters such as temperature and the CO₂/H₂O ratio displayed low inclusion frequencies and limited statistical significance. In the NO-LAG model, seismicity and equilibrium pressure of hydrothermal system emerged as dominant predictors, explaining 41.8% and 34.7% of the variance, respectively. Figure 8 illustrates nonlinear relationships and inverted trends among the five most relevant predictors, including negative associations for temperature and pressure, and positive trends for heat flow. Despite strong individual contributions, the NO-LAG model underperformed relative to the LAG model across all evaluation metrics. 3.2 Identification of critical transitions using GCPA The analysis was performed over two overlapping time windows: 2018–2023 and 2018–2024. The first interval was selected to identify transitional behaviors prior to the onset of the 2023 critical point (CP2), while the extended window allows for the detection of the additional system-wide reorganization that characterizes CP2. This approach supports a robust validation of the Global Critical Point Analysis (GCPA) framework across different temporal baselines. We used the outputs of the MFPA as input for the Global Critical Point Analysis (GCPA), focusing on six variables: ground deformation, seismicity, temperature, heat flow, CO₂ flux, and equilibrium pressure of hydrothermal system P(CO-CO₂-H₂-H₂O). We ran GCPA over two time intervals: 2018–2023 and 2018–2024. The analysis revealed two distinct critical transitions: CP1 – November 30th, 2020 ( Fig. 9 A ) : This transition was characterized by synchronized peaks in CO₂ flux (> 3.5 × 10⁴ g/m²/day), equilibrium pressure of hydrothermal system. Seismic activity showed a sustained increase, and the rate of ground deformation rose sharply, indicating intensified magmatic degassing and increased pressurization within the system. CP2 – April 1st, 2023 ( Fig. 9 B ) : This event coincided with a new phase of unrest. Seismicity peaked, ground deformation accelerated, and multiple geochemical indicators as CO₂ Flux and equilibrium pressure, reached critical values. Finally, we highlight that the integration of individual time series into composite metrics, as shown in Fig. 9 , enabled the identification of significant changes and evolutions within the volcanic system. Figure 9 A highlights the synchronization of variables leading to CP1, while Fig. 9 B illustrates the multivariable interaction underlying CP2. These panels effectively demonstrate how different parameters interact and highlight critical transitions. In this context, the interpretation of the acquired data in light of the performed analysis results highlights that: On the temporal evolution of monitored parameters: Fig. 9 A shows that, approaching CP1, CO₂ flux, equilibrium temperature, ground deformation, seismicity, land surface temperature, and heat flux all exhibited a coordinated increase. This synchronized behavior suggests a coherent systemic shift leading to the first critical point. On the behavior of variables near CP2: Fig. 9 B indicates that, prior to CP2, trends across the same set of variables became more differentiated. While CO₂ flux, seismicity, and ground deformation peaked distinctly, other parameters such as heat flux and land surface temperature showed more irregular evolution. This divergence suggests a more complex, possibly fragmented transition. On the contribution of the composite curve: The composite curves in both panels effectively integrate the normalized trajectories of all included parameters. These aggregated trends clearly delineate the timing and intensity of CP1 and CP2, serving as robust indicators of multiparametric transitions. On the identification of critical transitions: The temporal co-evolution of CO₂ flux, equilibrium temperature, ground deformation, seismicity, land surface temperature, and heat flux confirms the occurrence of two major system reorganizations CP1 and CP2 marking distinct phases in the recent dynamics of the Solfatara-Pisciarelli volcanic system. 4. Discussion 4.1 The role of time-lagged regression in ground deformation modeling The proposed methodology allows us to reveal a strong overall association between the variables of interest and the set of factors (i.e., physical parameters) tested in the models. Specifically, regarding ground deformation, seismic activity plays a significant role, accounting for one-third of the overall variance even when a lag term is included in the model. In other words, it explains a substantial portion of the variance in a model while also accounting for unmeasured associative factors. This evidence is further supported by the high stability of seismic activity in the final models, as indicated by its consistently high frequency of statistical significance (BIF above 97%, Table 2 ) in the bootstrapped datasets. In models without the LAG term, other factors gain prominence in the association, attempting to compensate for the 'unmeasured' influences. The high BIF values observed for almost all considered variables in these models reflect this compensatory effect. In particular, for a substantial number of factors (except seismic activity), we observed non-significant levels paired with high BIF values, and vice versa. This suggests that none of these factors, aside from seismic activity, establish a stable association. The LAG model effectively accounts for the delays between independent variables in the historical time series. It demonstrates a strong correspondence between observed data ("Raw data") and the model fit ("Fit"), with residuals ("Residual") remaining close to zero (Fig. 6 A). This indicates that the model successfully captures the temporal delays between independent variables (seismicity, heat flux, and CO₂ flux) and ground deformation. By including a LAG term, the model improves predictive accuracy by accounting for the dynamics of the cause-effect relationships between these independent parameters and the dependent variable. In contrast, in the NO-LAG model, the residuals exhibit greater variability and a systematic temporal oscillation, suggesting that the model fails to adequately explain the temporal delays between the independent variables and ground deformation. In complex systems like the CF volcanic field, processes such as fluid transport or gas release require time to propagate, leading to delays between the cause (variations in independent variables) and the effect (ground deformation). The inclusion of this time interval in the lag model allows for a more accurate representation of these dynamic processes, resulting in a closer fit between observed data and the model. The near-stationary residuals close to zero in Fig. 6 A further confirm the model's effectiveness in capturing temporal delays. 4.2 Identification of critical transitions through GCPA Furthermore, the application of the GCPA has been fundamental in identifying the timing and characteristics of critical transitions within the volcanic system. The GCPA method enabled us to detect two major critical points (CP1: November 30th 2020, CP2: April 1st 2023), which correspond to significant variations in multiple monitored parameters, including deformation rates, CO₂ flux, and seismicity. These transitions suggest that the SP system underwent two significant evolutionary phases, likely reflecting pressure accumulation and subsequent adjustments within the magmatic-hydrothermal environment. The construction of composite metrics by integrating multiple normalized time series allowed us to observe the synchronized evolution of key geophysical and geochemical parameters. This aggregation process revealed collective behaviors and system-wide responses, such as those leading to CP1 and CP2, that would have been difficult to identify through the analysis of individual variables alone. By emphasizing moments of convergence among signals, composite metrics provide a valuable tool for detecting emerging instability in complex volcanic systems. The first transition (CP1) marked the onset of an instability phase, characterized by an abrupt increase in gas fluxes and localized ground deformation acceleration. A subsequent stabilization phase followed, leading up to the second transition (CP2). It was associated with an additional shift in seismic and geochemical parameters, indicating renewed activity in the SP system. Identifying such transitions is crucial for volcanic hazard assessment, as it provides insights into the system’s evolution over time and highlights periods of potential increased unrest. 4.3 Methodological advantages over classical and non-parametric approaches One of the major strengths of this study is the integration of advanced statistical modeling (MFPA) with critical transition detection (GCPA), providing a novel framework for revealing instability phase in a volcanic system. This approach offers three main advantages: refined distinction between transient fluctuations and systemic transitions. The methodology may reduce the risk of false positives in hazard assessment by distinguishing between short-term variations and significant shifts in system dynamics. detection of major critical transitions. The identification of two major critical transitions, corresponding to significant changes in deformation rates, CO₂ flux, and seismicity, highlights the importance of data-driven approaches for monitoring evolving volcanic conditions. scalability to other volcanic systems. The validated methodology can be extended to other volcanic systems, enhancing global monitoring strategies by applying similar multiparametric analyses in high-risk calderas. 4.4 Link between critical transitions and overpressure source Our results allow to make a connection between the identified critical transition points (CP1 and CP2) and the evolution of the shallow overpressure source presented by (20) achieved using 4D geodetic tomographic inversion techniques. The GCPA results highlight a distinct critical point corresponding to a major shift in system dynamics, indicating a transition from stable conditions to an accelerated phase of pressure accumulation and deformation. The temporal coincidence between the GCPA-derived critical point and the inferred migration of the overpressure source suggests a strong linkage. This is supported by the observed increase in CO₂ flux and seismic activity recorded prior to the modeled overpressure source migration, as well as by the acceleration in ground deformation patterns detected through InSAR data. Specifically, the increase in CO₂ flux, deformation rates, and seismicity, key indicators in the GCPA model, precedes the modeled ascent of pressurized fluids. This supports the hypothesis that the critical point represents the onset of enhanced permeability and fluid mobilization, facilitating the upward migration of the over-pressured magmatic-hydrothermal mixture. 4.5 Applicability and perspectives for development The integrated MFPA–GCPA framework presented here is not restricted to the Solfatara-Pisciarelli area, but can be generalized to other volcanic systems with dense multiparametric monitoring. The core advantage lies in the model's ability to capture nonlinear and time-delayed interactions between diverse geophysical and geochemical signals, without relying on pre-set thresholds or specific event forecasts. For example, volcanoes such as Yellowstone, Long Valley, or Taal — where persistent deformation, degassing, and shallow seismicity coexist — could benefit from this approach, particularly in identifying early signs of systemic reorganization or slow-onset unrest. Similarly, densely monitored stratovolcanoes like Etna or Stromboli, where cyclic eruptive behaviors are observed, may provide ideal testbeds for implementing time-lagged multivariate models. The case of Campi Flegrei is particularly suited to demonstrate this methodology, due to the exceptional continuity of time series, the known recurrence of hydrothermal/magmatic transitions, and the availability of high-quality deformation, gas and seismic datasets. In this context, the proposed framework enhances our ability to recognize critical transitions that may remain hidden in traditional analyses. We anticipate that applying this approach across different volcanic settings will improve comparative analyses, foster cross-caldera learning, and contribute to harmonized early warning protocols at global scale. The accuracy of the identified critical point depends on the resolution and consistency of input datasets, and uncertainties in model assumptions regarding fluid migration dynamics may affect interpretations. Further validation through independent geophysical and geochemical observations is necessary to refine these findings for identifying and interpreting critical transitions in volcanic systems. This study demonstrates that the combination of MFPA and GCPA provides an effective tool for understanding the dynamics of complex volcanic systems. The identification of critical transitions in the SP area underscores the potential of this methodology in early warning systems and risk assessment frameworks. Furthermore, the introduction of a time-lagged model significantly improves the predictive power of ground deformation analysis, capturing key interactions between geophysical and geochemical parameters. These insights are particularly relevant for optimizing volcanic monitoring networks and anticipating hazardous events with greater accuracy. Compared to non-parametric approaches such as decision trees, random forests, or k-nearest neighbors (k-NN), the Multivariable Fractional Polynomial Analysis (MFPA) offers several advantages. First, it allows the identification of interpretable functional relationships between variables, including the assessment of non-linear trends within a parametric framework. This aspect is crucial when the physical meaning of coefficients and lag structures is necessary to relate model output to volcanic processes. Second, the MFPA enables the partitioning of explained variance and the decomposition of predictor contributions (via R²p and BIF), providing quantitative insights into which parameters are most strongly associated with observed dynamics. In contrast, non-parametric methods are often “black-box” in nature, lacking direct interpretability unless paired with post-hoc explainability tools Third, parametric models are computationally efficient and robust when applied to moderate-sized time series, as in our case, and they require fewer tuning parameters or training iterations. This makes them particularly suitable for real-time applications and for systems where interpretability and parsimony are as critical as predictive power. Nonetheless, future developments may incorporate ensemble strategies where MFPA is complemented by non-parametric classifiers or anomaly detectors to enhance forecasting capabilities and sensitivity to rare events. 5. Conclusions The joint application of GCPA and MFPA represents significant methodological advancement, as it enables both the identification of non-linear interactions between variables and the detection of critical transitions in volcanic activity. Modeling temporal delays and systemic tipping points is crucial for understanding the evolution of volcanic unrest and supporting decision-making processes in hazard mitigation. To further develop this approach, future research should focus on: (i) expanding the dataset to include real-time monitoring data and high-resolution satellite imagery for enhanced early warning signals detection; (ii) applying this methodology to other active volcanic systems to validate its robustness and scalability across different geological contexts; (iii) developing hybrid models that integrate AI-driven pattern recognition with our current framework to enhance forecasting capabilities. This study utilized a uniquely detailed dataset for the SP area, enabling a critical step in understanding volcanic system interactions through a specifically calibrated one-dimensional model. This approach paves the way for future multidimensional models and broader regional applications, potentially refining volcanic monitoring and hazard assessment by integrating geophysical and geochemical parameters. In conclusion, this research represents a step forward in volcanic hazard assessment, demonstrating the feasibility of a data-driven, multiparametric approach to detecting critical transitions in active volcanic systems. The integration of advanced statistical methods with real-world monitoring data can substantially improve early warning systems, ultimately enhancing the safety and resilience of communities living in volcanic regions. In this context, in recent years, the identification of early warning signals (EWS) for critical transitions in volcanic systems has become a central goal for improving monitoring capabilities and response preparedness. Several studies have explored this challenge through approaches ranging from tipping point analysis in complex systems to stochastic modeling and machine learning techniques. Early warning indicators based on statistical metrics such as increasing variance, lag-1 autocorrelation (critical slowing down), and kurtosis have shown promising results in specific case studies (Scheffer et al., 2009 ; Dakos et al., 2012 ; Lenton, 2013 ). However, their application in multiparametric systems is often hindered by high noise sensitivity and limited robustness under complex dynamics (Bernuzzi & Boers 2022 ; Gualandi and Liu, 2021 ; Segall, 2019 ; Boers, 2021 ). In comparison, the integrated MFPA–GCPA framework presents three key advantages: (i) It explicitly models nonlinear and time-delayed relationships between geophysical and geochemical parameters—an aspect rarely addressed in conventional EWS models, which often assume Markovian and instantaneous dynamics (Kéfi et al., 2014 ). (ii) It avoids arbitrary thresholds: critical points are identified through mathematical optimization over normalized signals, thereby reducing the risk of false positives due to random fluctuations ( Boettiger & Hastings, 2013 ). (iii) It enables the quantitative and temporal integration of multiple signals into a single composite metric, rather than treating them separately, thus facilitating the detection of systemic transitions rather than isolated local events—an essential feature in highly interdependent systems such as volcanic ones (Newhall & Pallister 2015 ; Behr et al., 2024 ). The MFPA–GCPA framework thus represents a hybrid approach: more flexible than classical descriptive models, yet more interpretable and statistically robust than many artificial intelligence-based methods (Bountos et al., 2022 ; Behr et al., 2024 ). This combination of interpretability and methodological rigor makes it especially well-suited for real-time volcanic surveillance, where both reliability and transparency of signals are critical. From an operational perspective, the application of this framework offers significant potential as a decision-support tool for civil protection authorities and monitoring agencies, such as INGV or global volcanic alert centers. Specifically: The inclusion of time lags enables earlier detection of transitional phases compared to conventional signals, thereby improving the responsiveness of early warning systems (Papale, 2017 ; Petrosino et al., 2023 ). The identification of systemic transitions (CP1 and CP2) provides a temporal map of structural changes within the system, which can inform dynamic updates to hazard levels (Petrosino et al., 2023 ). The construction of normalized composite metrics allows for seamless integration into monitoring dashboards, enhancing the interpretability of signals even for non-technical users (GVMID working group, 2024 ; Potter et al., 2014 ). The proposed approach could be implemented in near-real time by integrating existing data streams (e.g., InSAR, seismicity, gas emissions) with an automated MFPA–GCPA analysis module. This would enhance the predictive capabilities of surveillance networks, providing a robust mechanism for detecting critical phases, even in the absence of clear eruptive precursors (Ebmeier et al., 2018 ; Bevilacqua et al., 2020 ). Finally, the method's adaptability to diverse geological settings makes this framework a strong candidate for the standardization of next-generation multidisciplinary monitoring strategies, particularly in high-risk caldera systems such as Campi Flegrei, Yellowstone, Taupo, and Toba (Newhall & Pallister, 2015 ; Biggs et al., 2016 ). Declarations Acknowledgments We thank the DTA.AD004.065 project, "Development and applications of geophysical methods based on the use of environmental remote sensing data," and its subproject DTA.AD004.065.001, "Geophysical analysis and modeling of environmental processes," for partially supported this research activity. UAVs thermal data were supported both by the DPC-INGV B2 project WP2 Task 5 (2019-2021) and by the DPC-INGV agreement Annex A. Seismic data are supported by the DPC-INGV agreement Annex A. Funding The authors received no funding for this work. Author contributions: Vitale A (AV), Barone A (AB), Marotta E (EM), Vitale D F (DFV), Pepe S (SP), Peluso R (RP), Castaldo R (RC), Avino R (RA), Mercogliano F (FM), Pepe A (AP), Accomando F (FA), Avvisati G (GA), Belviso P (PB), Bellucci Sessa E (EBS), Carandente A (AC), Perrini M (MP), Sansivero F (FS), Tizzani P (PT) Conceptualization: PT, AV, DFV Methodology: PT, AV, DFV, EM Investigation: AB, AP, SP, RA, AC, PB, GA, EBS, FS, RC, FM, EM,FA, RP Visualization: PT, AV, AB, MP Supervision: PT, AV, EM Writing original draft: PT, AV, DFV, EM, AB, RC, FS, MP, SP, RP, EBS Writing review & editing: Data and materials availability: “The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.” References Avino, R., Peluso, R., Marotta, E., Belviso, P., Caliro, S., Carandente, A., Chiodini, G., Macedonio, G., Avvisati, G., Marfè, B., Thermal energy release measurement with thermal camera: The case of La Solfatara Volcano (Italy). Remote Sens. 11, 167 (2019). https://doi.org/10.3390/rs11020167 Barone, A., Gola, G., Caliro, S., Chiodini, G., Tizzani, P., Castaldo, R., Long-term thermo-fluid dynamic modeling of Solfatara hydrothermal system, Campi Flegrei caldera. J. Volcanol. Geotherm. Res. 459, 108277 (2025). https://doi.org/10.1016/j.jvolgeores.2025.108277 Behr Y, Cerminara M, Rinaldi AP, Azzaro R. Probabilistic, multi-sensor eruption forecasting: a Bayesian framework. Geophys Res Lett . 2024;51:e2024GL112029. https://doi.org/10.1029/2024GL112029 Bernuzzi S, Boers N. Critical slowing-down in dynamical systems driven by non-stationary correlated noise. Phys Rev Res . 2022;4:013230. https://doi.org/10.1103/PhysRevResearch.4.013230 Berardino, P., Fornaro, G., Lanari, R., Sansosti, E., A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 40, 2375–2383 (2002). https://doi.org/10.1109/TGRS.2002.803792 Bevilacqua A, Neri A, Genco R, Isaia R, Marzocchi W. Failure forecast on Campi Flegrei GPS data 2011–2020. arXiv preprint . 2020. arXiv:2007.02756. Available from: https://arxiv.org/abs/2007.02756 Bevilacqua, A., Neri, A., De Martino, P., Isaia, R., Novellino, A., Tramparulo, F., Vitale, S., Radial interpolation of GPS and leveling data of ground deformation in a resurgent caldera: Application to Campi Flegrei (Italy). J. Geod. 94, 13 (2020). https://doi.org/10.1007/s00190-020-01355-x Bevilacqua, A., Neri, A., De Martino, P., Giudicepietro, F., Ricciolino, P., Patra, A., Pitman, E. B., Bursik, M., Voight, B., Flandoli, F., Macedonio, G., Accelerating upper crustal deformation and seismicity of Campi Flegrei caldera (Italy) during the 2000–2023 unrest. Commun. Earth Environ. 5, 742 (2024). https://doi.org/10.1038/s43247-024-01865-y Biggs J, Pritchard ME. Global satellite observations of magmatic and volcanic deformation: implications for volcano monitoring. In: Pyle D, Mather TA, Biggs J, editors. Global Volcanic Hazards and Risk . Cambridge: Cambridge University Press; 2016. p. 1–34. Boers N. Observation-based early-warning signals for AMOC collapse. Nat Clim Chang . 2021;11:680–688. https://doi.org/10.1038/s41558-021-01181-1 Boettiger C, Ross N, Hastings A. Early warning signals: The charted and uncharted territories. Theor Ecol . 2013;6:255–264. https://doi.org/10.1007/s12080-013-0192-6 Bountos N, Andronopoulos S, Tzovaras D. Self-supervised contrastive learning for volcanic unrest detection. arXiv preprint . 2022. arXiv:2202.04030. Available from: https://arxiv.org/abs/2202.04030 Boyd, S., Vandenberghe, L., Convex Optimization (Cambridge Univ. Press, Cambridge, UK, 2004). https://doi.org/10.1017/CBO9780511804441 Caliro, S., Avino, R., Capecchiacci, F., Carandente, A., Chiodini, G., Cuoco, E., Minopoli, C., Rufino, F., Santi, A., Rizzo, A. L., Aiuppa, A., Allocca, V., De Vita, P., Di Vito, M. A., Chemical and isotopic characterization of groundwater and thermal waters from the Campi Flegrei caldera (southern Italy). J. Volcanol. Geotherm. Res. 460, 108280 (2025). https://doi.org/10.1016/j.jvolgeores.2025.108280 Caliro, S., Chiodini, G., Avino, R., Carandente, A., Cuoco, E., Di Vito, M. A., Minopoli, C., Rufino, F., Santi, A., Lages, J., Mangiacapra, A., Monteleone, B., Pappalardo, L., Taracsák, Z., Tramati, C., Vizzini, S., Aiuppa, A., Escalation of caldera unrest indicated by increasing emission of isotopically light sulfur. Nat. Geosci. 18, 123–132 (2025). https://doi.org/10.1038/s41561-024-01632-w Cardellini, C., Chiodini, G., Frondini, F., Avino, R., Bagnato, E., Caliro, S., Lelli, M., Rosiello, A., Monitoring diffuse volcanic degassing during volcanic unrests: The case of Campi Flegrei (Italy). Sci. Rep. 7, 6757 (2017). https://doi.org/10.1038/s41598-017-06169-5 Castaldo, R., Tizzani, P., Solaro, G., Inflating source imaging and stress/strain field analysis at Campi Flegrei Caldera: The 2009–2013 unrest episode. Remote Sens. 13, 2298 (2021). https://doi.org/10.3390/rs13122298 Chiodini, G., Vandemeulebrouck, J., Caliro, S., D’Auria, L., De Martino, P., Mangiacapra, A., Petrillo, Z., Evidence of thermal driven processes triggering the 2005–2014 unrest at Campi Flegrei caldera. Earth Planet. Sci. Lett. 414, 58–67 (2015). https://doi.org/10.1016/j.epsl.2015.01.012 Chiodini, G., Giudicepietro, F., Vandemeulebrouck, J., Aiuppa, A., Caliro, S., De Cesare, W., Tamburello, G., Avino, R., Orazi, M., D’Auria, L., Fumarolic tremor and geochemical signals during a volcanic unrest. Geology 45, 1131–1134 (2017). https://doi.org/10.1130/G39447.1 Chiodini, G., Todesco, M., Caliro, S., Del Gaudio, C., Macedonio, G., Russo, M., Magma degassing as a trigger of bradyseismic events: The case of Phlegrean Fields (Italy). Geophys. Res. Lett. 30, 1434 (2003). https://doi.org/10.1029/2002GL016790 Chiodini, G., Avino, R., Caliro, S., Minopoli, C., Temperature and pressure gas geoindicators at the Solfatara fumaroles (Campi Flegrei). Ann. Geophys. 54, 2 (2011). https://doi.org/10.4401/ag-5002 Chiodini, G., Caliro, S., De Martino, P., Avino, R., Gherardi, F., Early signals of new volcanic unrest at Campi Flegrei caldera? Insights from geochemical data and physical simulations. Geology 40, 943–946 (2012). https://doi.org/10.1130/G33251.1 Chiodini, G., Caliro, S., Avino, R., Bagnato, E., Capecchiacci, F., Carandente, A., Allocca, V., Aiuppa, A., The hydrothermal system of the Campi Flegrei Caldera, Italy. in Campi Flegrei: A Restless Caldera in a Densely Populated Area, G. Orsi, M. D’Antonio, L. Civetta, Eds. (Springer, Berlin, 2022), pp. 239–255. https://doi.org/10.1007/978-3-642-37060-1_9 Cioni, R., Corazza, E., Medium temperature fumarolic gas sampling. Bull. Volcanol. 44, 23–29 (1981). Cirillo, F., Avvisati, G., Belviso, P., Marotta, E., Peluso, R., Pescione, P. R., STARTED (StaTistical Analysis clusteRed ThErmal Data) (version 1.0.1). Zenodo (2022). https://doi.org/10.5281/zenodo.5886707 Dakos V, Carpenter SR, Brock WA, Ellison AM, Guttal V, Ives AR, et al. Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. PLoS One . 2012;7(7):e41010. https://doi.org/10.1371/journal.pone.0041010 D’Auria, L., Pepe, S., Castaldo, R., Giudicepietro, F., Macedonio, G., Ricciolino, P., Tizzani, P., Casu, F., Lanari, R., Manzo, M., Martini, M., Sansosti, E., Zinno, I., Magma injection beneath the urban area of Naples: A new mechanism for the 2012–2013 volcanic unrest at Campi Flegrei caldera. Sci. Rep. 5, 13100 (2015). https://doi.org/10.1038/srep13100 Del Gaudio, C., Aquino, I., Ricciardi, G.P., Ricco, C., Scandone, R., Unrest episodes at Campi Flegrei: A reconstruction of vertical ground movements during 1905–2009. J. Volcanol. Geotherm. Res. 195, 48–56 (2010). https://doi.org/10.1016/j.jvolgeores.2010.05.014 Di Filippo, A., Peluso, R., SEREWrap – SERENADE Wrapper (version 5.3.2). Zenodo (2019). https://doi.org/10.5281/zenodo.14844623 Dirac, P. A. M., The Principles of Quantum Mechanics, 4th ed., Oxford, UK: Oxford University Press, 1958. Ebmeier SK, Elliott JR, Naismith A, Biggs J. Synthesis of global satellite observations of magmatic and volcanic deformation: implications for volcano monitoring. Geol Soc Lond Spec Publ . 2018;446:305–322. https://doi.org/10.1144/SP446.5 Ester, M., Kriegel, H. P., Sander, J., Xu, X., A density-based algorithm for discovering clusters in large spatial databases with noise. in Proc. 2nd Int. Conf. Knowl. Discovery Data Mining (KDD-96), (1996), pp. 226–231. Finkel, S. S., Linear panel analysis. in Handbook of Longitudinal Research: Design, Measurement, and Analysis, S. Menard, Ed. (Elsevier/Academic Press, Burlington, MA, 2007), pp. 486–487. Fletcher, R., Practical Methods of Optimization (Wiley, Chichester, ed. 2, 2013). Gabriel, A. K., Goldstein, R. M., Zebker, H. A., Mapping small elevation changes over large areas: Differential interferometry. J. Geophys. Res. 94, 9183–9191 (1989). https://doi.org/10.1029/JB094iB07p09183 Giacomuzzi, G., Chiarabba, C., Bianco, F., De Gori, P., Piana Agostinetti, N., Tracking transient changes in the plumbing system at Campi Flegrei Caldera. Earth Planet. Sci. Lett. 637, 118744 (2024). https://doi.org/10.1016/j.epsl.2024.118744 Giggenbach, W. F., A simple method for the collection and analysis of volcanic gas samples. Bull. Volcanol. 39, 132–145 (1975). Giggenbach, W. F., Gouguel, R. L., Collection and analysis of geothermal and volcanic water and gas discharges. Chemistry Division Report, DSIR, New Zealand, No. CD 2401 (1989). Giudicepietro, F., Avino, R., Bellucci Sessa, E., Bevilacqua, A., Bonano, M., Caliro, S., Chiodini, G., Burst-like swarms in the Campi Flegrei caldera accelerating unrest from 2021 to 2024. Nat. Commun. 16, 1548 (2025). https://doi.org/10.1038/s41467-025-39587-1 Goldberg, D. E., Genetic Algorithms in Search, Optimization, and Machine Learning, Reading, MA, USA: Addison-Wesley, 1989. Gualandi A, Liu Z. Variational Bayesian Independent Component Analysis for InSAR displacement time-series: application to the 2009 Sierra Negra unrest. J Geophys Res Solid Earth . 2021;126:e2020JB020845. https://doi.org/10.1029/2020JB020845 GVMID Working Group. The Global Volcano Monitoring Infrastructure Database: toward integrated dashboards. Front Earth Sci . 2024;12:1284889. https://doi.org/10.3389/feart.2024.1284889 Isaia, R., Vitale, S., Di Giuseppe, M. G., Iannuzzi, E., Tramparulo, F. D., Troiano, A., Stratigraphy, structure, and volcanotectonic evolution of Solfatara maar-diatreme (Campi Flegrei, Italy). Geol. Soc. Am. Bull. 127, 1485–1504 (2015). https://doi.org/10.1130/B31183.1 Johnson, C. E., Bittenbinder, A., Bogaert, B., Dietz, L., Kohler, W., Earthworm: A flexible approach to seismic network processing. IRIS Newsletter 14, 1–4 (1995). Jolivet, R., Grandin, R., Lasserre, C., Doin, M. P., Peltzer, G., Systematic InSAR tropospheric phase delay corrections from global meteorological reanalysis data. Geophys. Res. Lett. 38, L17311 (2011). https://doi.org/10.1029/2011GL048757 Karmarkar, N., A new polynomial-time algorithm for linear programming. Combinatorica 4, 373–395 (1984). https://doi.org/10.1007/BF02579150 Kéfi S, Guttal V, Brock WA, Carpenter SR, Ellison AM, Livina V, et al. Early warning signals of ecological transitions: Methods for spatial patterns. PLoS One . 2014;9(3):e92097. https://doi.org/10.1371/journal.pone.0092097 Kutner, M. H., Nachtsheim, C. J., Neter, J., Li, W., Applied Linear Statistical Models, 5th ed. (McGraw-Hill/Irwin, New York, NY, 2005), p. 409. Lee, W. H. K., Lahr, J. C., HYPO71: A computer program for determining hypocenter, magnitude, and first motion pattern of local earthquakes. U.S. Geol. Surv. Open File Rep. 75–311 (1975). Lenton TM. Environmental tipping points. Annu Rev Environ Resour . 2013;38:1–29. https://doi.org/10.1146/annurev-environ-102511-084654 Manzo, M., Ricciardi, G. P., Casu, F., Ventura, G., Zeni, G., Borgstrom, S., Berardino, P., Del Gaudio, C., Lanari, R., Surface deformation analysis in the Ischia island (Italy) based on spaceborne radar interferometry. J. Volcanol. Geotherm. Res. 151, 399–416 (2006). https://doi.org/10.1016/j.jvolgeores.2005.09.010 Marotta, E., Peluso, R., Avino, R., Avvisati, G., Bellucci Sessa, E., Belviso, P., Caputo, T., Carandente, A., Cirillo, F., Pescione, R. A., Clusterisation and temporal trends of heat flux by UAS thermal camera. Remote Sens. 16, 1102 (2024). https://doi.org/10.3390/rs16061102 Marotta E., R. Peluso, R. Avino, P. Belviso, S. Caliro, A. Carandente, G. Chiodini, G. Macedonio, G. Avvisati, B. Marfè, Thermal energy release measurement with thermal camera: The case of La Solfatara Volcano (Italy). Remote Sens. 11, 167 (2019). https://doi.org/10.3390/rs11020167 Massonnet, D., Rossi, M., Carmona, C., Ardagna, F., Peltzer, G., Feigl, K., Rabaute, T., The displacement field of the Landers earthquake mapped by radar interferometry. Nature 364, 138–142 (1993). https://doi.org/10.1038/364138a0 Mercogliano, F., Barone, A., D’Auria, L., Castaldo, R., Silvestri, M., Bellucci Sessa, E., Caputo, T., Stroppiana, D., Caliro, S., Minopoli, C., Avino, R., Tizzani, P., Thermal patterns at the Campi Flegrei caldera inferred from satellite data and independent component analysis. Remote Sens. 16, 4615 (2024). https://doi.org/10.3390/rs16234615 Newhall CG, Pallister JS. Introducing the Volcanic Unrest Index (VUI): a tool to quantify and communicate the intensity of volcanic unrest. Bull Volcanol . 2015;77:77. https://doi.org/10.1007/s00445-015-0952-y Oppenheim, A. V., Schafer, R. W., Discrete-Time Signal Processing, Upper Saddle River, NJ, USA: Prentice-Hall, 1999. Orsi, G., D’Antonio, M., Civetta, L., Volcanic and deformation history of the Campi Flegrei volcanic field, Italy. in Campi Flegrei, Eds. (Springer, Berlin, 2022), pp. 1–28. https://doi.org/10.1007/978-3-642-37060-1_1 Papale P. Rational volcanic hazard forecasts and the use of volcanic alert levels. J Appl Volcanol . 2017;6(1):4. https://doi.org/10.1186/s13617-017-0064-7 Peluso, R., Il database sismico SERENADE: un sistema REST per la gestione delle localizzazioni sismiche. Miscellanea INGV 57, 7 (2020). https://doi.org/10.13127/misc/57/7 Peluso, R., SEDAN – Seismic Event Dispatcher And Notifier (version 0.11.1). Zenodo (2019). https://doi.org/10.5281/zenodo.14845439 Peluso, R., SERENADE – SEismic Restful ENAbled DatabasE (version 0.6.1). Zenodo (2018a). https://doi.org/10.5281/zenodo.14732232 Peluso, R., WESSEL – WEb interface for Seismic Signal Elaboration and Location (version 1.12.1). Zenodo (2018b). https://doi.org/10.5281/zenodo.14727127 Pepe, A., Solaro, G., Calò, F., Dema, C., A minimum acceleration approach for the retrieval of multiplatform InSAR deformation time series. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 9, 3883–3898 (2016). https://doi.org/10.1109/JSTARS.2016.2577878 Pepe, S., De Siena, L., Barone, A., Castaldo, R., D’Auria, L., Manzo, M., Casu, F., Fedi, M., Lanari, R., Bianco, F., Tizzani, P., Volcanic structures investigation through SAR and seismic interferometric methods: The 2011–2013 Campi Flegrei unrest episode. Remote Sens. Environ. 234, 111440 (2019). https://doi.org/10.1016/j.rse.2019.111440 Petrillo, Z., Chiodini, G., Mangiacapra, A., Caliro, S., Capuano, P., Russo, G., Cardellini, C., Avino, R., Defining a 3D physical model for the hydrothermal circulation at Campi Flegrei caldera (Italy). J. Volcanol. Geotherm. Res. 264, 172–183 (2013). https://doi.org/10.1016/j.jvolgeores.2013.08.008 Petrosino S, De Lauro E, Aquino I, Del Pezzo E. Evolution in unrest processes at Campi Flegrei inferred from multiparametric lags. J Volcanol Geotherm Res . 2023;431:107784. https://doi.org/10.1016/j.jvolgeores.2023.107784 Potter SH, Becker JS, Johnston DM, Rossiter KP. Communicating the status of volcanic activity: revising New Zealand’s VAL system. J Appl Volcanol . 2014;3:13. https://doi.org/10.1186/s13617-014-0013-6 Rabinowitz, P. H., Critical point theory and its applications to differential equations: A survey. in Proc. Symp. Pure Math. 45, 357–386 (1986). R Core Team: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Avaiable on: https://www.R-project.org/ (2014) Ricciolino, P., Lo Bascio, D., Esposito, R., GOSSIP – Database sismologico Pubblico INGV-Osservatorio Vesuviano. Istituto Nazionale di Geofisica e Vulcanologia (INGV) (2024). https://doi.org/10.13127/gossip Rosen, P. A., Hensley, S., Joughin, I. R., Li, F. K., Madsen, S. N., Rodriguez, E., Goldstein, R., Synthetic aperture radar interferometry. Proc. IEEE 88, 333–382 (2000). https://doi.org/10.1109/5.838084 Royston, P., Sauerbrei, W., Multivariable Model-Building: A Pragmatic Approach to Regression Analysis Based on Fractional Polynomials for Modelling Continuous Variables (Wiley, Chichester, 2008). Scotto di Uccio, F., Lomax, A., Natale, J., Muzellec, T., Festa, G., Nazeri, S., Convertito, V., Bobbio, A., Strumia, C., Zollo, A., Delineation and fine-scale structure of fault zones activated during the 2014–2024 unrest at the Campi Flegrei caldera (Southern Italy) from high-precision earthquake locations. Geophys. Res. Lett. 51, e2023GL107680 (2024). https://doi.org/10.1029/2023GL107680 Scheffer, M., Bascompte, J., Brock, W. A., Carpenter, S. R., Dakos, V., van Nes, E. H., and Sugihara, G. (2009). Early-warning signals for critical transitions. Nature, 461(7260), 53-59. https://doi.org/10.1038/nature08227 Shorrocks A. F., Decomposition procedures for distributional analysis: A unified framework based on the Shapley value. J. Econ. Inequal. 11, 99–126 (2013). https://doi.org/10.1007/s10888-011-9184-0 Shumway, R. H., Stoffer, D. S., Time Series Analysis and Its Applications: With R Examples, 4th ed. (Springer, New York, 2017). https://doi.org/10.1007/978-3-319-52452-8 Segall P. Noise sources and detection limits in volcano geodesy. J Volcanol Geotherm Res . 2019;381:106624. https://doi.org/10.1016/j.jvolgeores.2019.106624 Siniscalchi, A., Tripaldi, S., Romano, G., Chiodini, G., Improta, L., Petrillo, Z., D’Auria, L., Caliro, S., Avino, R., Reservoir structure and hydraulic properties of the Campi Flegrei geothermal system inferred by audiomagnetotelluric. J. Geophys. Res. Solid Earth 124, 5336–5356 (2019). https://doi.org/10.1029/2018JB016514 Smith, V. C., Isaia, R., Pearce, N. J. G., Tephrostratigraphy and glass compositions of post-15 kyr Campi Flegrei eruptions: Implications for eruption history and chronostratigraphic markers. Quat. Sci. Rev. 30, 3638–3660 (2011). https://doi.org/10.1016/j.quascirev.2011.07.012 Steinmann, L., Spiess, V., Sacchi, M., The Campi Flegrei caldera (Italy): Formation and evolution in interplay with sea-level variations since the Campanian Ignimbrite eruption at 39ka. J. Volcanol. Geotherm. Res. 327, 361–374 (2016). https://doi.org/10.1016/j.jvolgeores.2016.09.001 Suleiman, A. A., Analysis of multicollinearity in multiple regressions. Int. J. Adv. Technol. Eng. Sci. 3, 571–578 (2015). Sze, S. M., Physics of Semiconductor Devices, 2nd ed., New York, NY, USA: Wiley, 1981. Theußl, S., Schwendinger, F., Hornik, K., ROI: An extensible R optimization infrastructure. J. Stat. Softw. 94, 1–64 (2020). https://doi.org/10.18637/jss.v094.i15 Tizzani, P., Fernández, J., Vitale, A., Escayo, J., Barone, A., Castaldo, R., Pepe, S., De Novellis, V., Solaro, G., Pepe, A., Tramelli, A., Hu, Z., Samsonov, S. V., Vigo, I., Tiampo, K. F., Camacho, A. G., 4D imaging of the volcano feeding system beneath the urban area of the Campi Flegrei caldera. Remote Sens. Environ. 315, 114480 (2024). https://doi.org/10.1016/j.rse.2024.114480 Tramelli, A., Convertito, V., Godano, C., B value enlightens different rheological behaviour in Campi Flegrei caldera. Commun. Earth Environ. 5, 1447 (2024). https://doi.org/10.1038/s43247-024-01447-y Troiano, A., Di Giuseppe, M. G., Isaia, R., 3D structure of the Campi Flegrei caldera central sector reconstructed through short-period magnetotelluric imaging. Sci. Rep. 12, 20802 (2022). https://doi.org/10.1038/s41598-022-24998-6 Troiano, A., Isaia, R., Di Giuseppe, M. G., Tramparulo, F., Mormone, F., Vitale, V., D’Auria, L., Ricciardi, G. P., Deep electrical resistivity tomography for a 3D picture of the most active sector of Campi Flegrei caldera. Sci. Rep. 9, 15124 (2019). https://doi.org/10.1038/s41598-019-51568-0 Vanorio, A., Geremia, G., De Landro, M., Guo, Z., The recurrence of geophysical manifestations at the Campi Flegrei caldera. Sci. Adv. 11, eadt2067 (2025). https://doi.org/10.1126/sciadv.adt2067 Wickham, H., ggplot2: Elegant Graphics for Data Analysis (Springer, New York, 2009). https://doi.org/10.1007/978-0-387-98141-3 Additional Declarations No competing interests reported. Supplementary Files SupplementaryMaterial.docx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-7290139","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":500334765,"identity":"bed08e41-f1c0-45de-ac18-b1690c1978b9","order_by":0,"name":"Andrea Vitale","email":"","orcid":"","institution":"Institute for Agricultural and Forest Systems in the Mediterranean","correspondingAuthor":false,"prefix":"","firstName":"Andrea","middleName":"","lastName":"Vitale","suffix":""},{"id":500334766,"identity":"5202301e-1470-4cf2-a89a-e4d122c47849","order_by":1,"name":"Andrea Barone","email":"","orcid":"","institution":"Istituto per il Rilevamento Elettromagnetico dell'Ambiente","correspondingAuthor":false,"prefix":"","firstName":"Andrea","middleName":"","lastName":"Barone","suffix":""},{"id":500334767,"identity":"a82b67c0-6e22-45aa-b4ac-8f2f54983418","order_by":2,"name":"Enrica Marotta","email":"","orcid":"","institution":"National Institute of Geophysics and Volcanology","correspondingAuthor":false,"prefix":"","firstName":"Enrica","middleName":"","lastName":"Marotta","suffix":""},{"id":500334768,"identity":"7ea07499-4421-4a0e-8bdd-b783097f8cf1","order_by":3,"name":"Dino Franco Vitale","email":"","orcid":"","institution":"Casa di Cura San Michele","correspondingAuthor":false,"prefix":"","firstName":"Dino","middleName":"Franco","lastName":"Vitale","suffix":""},{"id":500334769,"identity":"3849d4be-cf14-4cad-a488-d2dbe795afb6","order_by":4,"name":"Susi Pepe","email":"","orcid":"","institution":"Istituto per il Rilevamento Elettromagnetico dell'Ambiente","correspondingAuthor":false,"prefix":"","firstName":"Susi","middleName":"","lastName":"Pepe","suffix":""},{"id":500334770,"identity":"8e10cd8b-e89d-41ae-8d50-95ba4f3db567","order_by":5,"name":"Rosario Peluso","email":"","orcid":"","institution":"National Institute of Geophysics and Volcanology","correspondingAuthor":false,"prefix":"","firstName":"Rosario","middleName":"","lastName":"Peluso","suffix":""},{"id":500334771,"identity":"fd7e7a8b-d488-4881-a38c-230bc962b780","order_by":6,"name":"Raffaele Castaldo","email":"","orcid":"","institution":"Istituto per il Rilevamento Elettromagnetico dell'Ambiente","correspondingAuthor":false,"prefix":"","firstName":"Raffaele","middleName":"","lastName":"Castaldo","suffix":""},{"id":500334772,"identity":"88595b1c-27b9-4466-80b6-61f504db31cf","order_by":7,"name":"Rosario Avino","email":"","orcid":"","institution":"National Institute of Geophysics and Volcanology","correspondingAuthor":false,"prefix":"","firstName":"Rosario","middleName":"","lastName":"Avino","suffix":""},{"id":500334773,"identity":"b174e5ff-7d7d-42d6-8d23-9344e42150de","order_by":8,"name":"Francesco Mercogliano","email":"","orcid":"","institution":"Parthenope University of Naples","correspondingAuthor":false,"prefix":"","firstName":"Francesco","middleName":"","lastName":"Mercogliano","suffix":""},{"id":500334774,"identity":"2077f487-60cb-4da2-84c1-2223460c12da","order_by":9,"name":"Antonio Pepe","email":"","orcid":"","institution":"Istituto per il Rilevamento Elettromagnetico dell'Ambiente","correspondingAuthor":false,"prefix":"","firstName":"Antonio","middleName":"","lastName":"Pepe","suffix":""},{"id":500334775,"identity":"afc5a21e-8dad-470b-91b6-4c2406d1262f","order_by":10,"name":"Filippo Accomando","email":"","orcid":"","institution":"Istituto per il Rilevamento Elettromagnetico dell'Ambiente","correspondingAuthor":false,"prefix":"","firstName":"Filippo","middleName":"","lastName":"Accomando","suffix":""},{"id":500334776,"identity":"d589860c-5cae-40ec-b1b8-efb400b66a2f","order_by":11,"name":"Gala Avvisati","email":"","orcid":"","institution":"National Institute of Geophysics and Volcanology","correspondingAuthor":false,"prefix":"","firstName":"Gala","middleName":"","lastName":"Avvisati","suffix":""},{"id":500334777,"identity":"608165c9-7d0d-48d9-82f9-a9304b20ba03","order_by":12,"name":"Pasquale Belviso","email":"","orcid":"","institution":"National Institute of Geophysics and Volcanology","correspondingAuthor":false,"prefix":"","firstName":"Pasquale","middleName":"","lastName":"Belviso","suffix":""},{"id":500334778,"identity":"15e526dd-309b-439f-aa3b-f58cc1a297c7","order_by":13,"name":"Eliana Bellucci Sessa","email":"","orcid":"","institution":"National Institute of Geophysics and Volcanology","correspondingAuthor":false,"prefix":"","firstName":"Eliana","middleName":"Bellucci","lastName":"Sessa","suffix":""},{"id":500334779,"identity":"a13d33cf-f649-41d0-9857-f9e43532ec7c","order_by":14,"name":"Antonio Carandente","email":"","orcid":"","institution":"National Institute of Geophysics and Volcanology","correspondingAuthor":false,"prefix":"","firstName":"Antonio","middleName":"","lastName":"Carandente","suffix":""},{"id":500334780,"identity":"56fb9976-4b3c-42db-b527-3889744dee69","order_by":15,"name":"Maddalena Perrini","email":"","orcid":"","institution":"Istituto per il Rilevamento Elettromagnetico dell'Ambiente","correspondingAuthor":false,"prefix":"","firstName":"Maddalena","middleName":"","lastName":"Perrini","suffix":""},{"id":500334781,"identity":"3390d93f-f369-45f3-a2de-e3d6e624d6e0","order_by":16,"name":"Fabio Sansivero","email":"","orcid":"","institution":"National Institute of Geophysics and Volcanology","correspondingAuthor":false,"prefix":"","firstName":"Fabio","middleName":"","lastName":"Sansivero","suffix":""},{"id":500334782,"identity":"f3deb7c6-9895-49ec-bc4f-4427c3f81a8a","order_by":17,"name":"Pietro Tizzani","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA3ElEQVRIiWNgGAWjYDACZijNBiYNbMDUAQJaGBuAlARUSxpEC349UC1QzmEGgtaYszMff1xRwVDHx3/84eeCgvPy5tIHGA9/wKPFspktsfHMGZDDDiRLzzC4bbizLwG/wwwO8xg2NrYBtTA2HJDmMbidYHCGgF8MDvN/hGhhZmz+zWNwjhgtPIwQLWzMbEBbDhDWAvSL4cyGMxKSbTxsbNY8BsmGG84AnXgGjxZz/sMPPjZU2PDL9x9/fJvnj528wRnmwx8q8DkMQkkgi4EjiqCWUTAKRsEoGAV4AAB0kEcjgUeN8QAAAABJRU5ErkJggg==","orcid":"","institution":"Istituto per il Rilevamento Elettromagnetico dell'Ambiente","correspondingAuthor":true,"prefix":"","firstName":"Pietro","middleName":"","lastName":"Tizzani","suffix":""}],"badges":[],"createdAt":"2025-08-04 10:38:33","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-7290139/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7290139/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":89345880,"identity":"e205da6e-f4e9-434c-a6dd-f072ed36cc03","added_by":"auto","created_at":"2025-08-19 04:41:45","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":3647734,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSolfatara-Pisciarelli volcanic area and geochemical monitoring network.\u003c/strong\u003e (A)provides an aerial view of the Campi Flegrei caldera and the main structural features of the study area from \u003cstrong\u003eBevilacqua et al. (2020)\u003c/strong\u003e; the red box indicates the enlarged study area in (B). (B) Outlines of the study area: Solfatara crater on the left and the Pisciarelli fumarolic field on the right.\u003c/p\u003e","description":"","filename":"Figure1.tif.jpg","url":"https://assets-eu.researchsquare.com/files/rs-7290139/v1/ecbc79eefe340c3c2ecc4b0d.jpg"},{"id":89346089,"identity":"e0982864-27b2-4884-bfc1-66025516520a","added_by":"auto","created_at":"2025-08-19 04:49:45","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":2946371,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eTemporal evolution of key geochemical parameters and CO₂ flux in the Solfatara-Pisciarelli system from 2018 to 2024.\u003c/strong\u003e (A)CO concentration time series. (B) Equilibrium temperature (T(CO-CO₂)) time series. (C) CO/CO₂ ratio. (D) Equilibrium pressure (P(CO-CO₂-H₂-H₂O)) time series indicating pressurization within the hydrothermal system. (E) CO₂/H₂O ratio reflecting variations in gas-water interaction and subsurface conditions. (F) CO₂/CH₄ ratio. (G) CO₂ flux time series (black) and its semiannual moving average (blue).\u003c/p\u003e","description":"","filename":"Figure2.tif.jpg","url":"https://assets-eu.researchsquare.com/files/rs-7290139/v1/98649e043db105660deec507.jpg"},{"id":89346453,"identity":"5845b00b-94ac-4bd5-ae65-824108579bf9","added_by":"auto","created_at":"2025-08-19 04:57:45","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":2210958,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eThermal anomalies and heat flow dynamics within the Pisciarelli area from UAV-based thermal mapping and time series analysis (2019-2024).\u003c/strong\u003e (A) Mean and annual thermal map from 2020 to 2024. (B): LST time series plot illustrating measured temperatures (black), de-seasonalized temperatures (red), and seasonal trends (green). (C): de-seasonalized temperature time series (red) and related six-month moving average (black). (D):Heat flux time series (blue) and heat release rates (red).\u003c/p\u003e","description":"","filename":"Figure3.tif.jpg","url":"https://assets-eu.researchsquare.com/files/rs-7290139/v1/04e20b0be76c183fb5159c72.jpg"},{"id":89346096,"identity":"b0a70c63-d07e-4b71-bf2d-ff9830b0601d","added_by":"auto","created_at":"2025-08-19 04:49:45","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":8937172,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSeismic clustering and temporal evolution of seismic activity within the CFc (2018-2024).\u003c/strong\u003e (A) Spatial distribution of seismic clusters within the CFc (Figs. S1-S6), mapped over the local topography. The different colors represent distinct clusters of seismic events. (B) North-South vertical cross-section showing the depth distribution of seismic clusters. (C) East-West cross-section of the seismic clusters. (D) Time series of the normalized number of seismic events per month.\u003c/p\u003e","description":"","filename":"Figure4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-7290139/v1/22c7c02e2130789436f8accb.jpg"},{"id":89346091,"identity":"c9fc6393-0fb2-462f-be8d-8910ccaeb054","added_by":"auto","created_at":"2025-08-19 04:49:45","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":2005474,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eVertical ground deformation in the Solfatara-Pisciarelli area based on InSAR measurements (2018-2024).\u003c/strong\u003e (A) InSAR mean vertical velocity map showing deformation rates across the CFc. The map highlights considered SAR pixels centered near the Pisciarelli spring, with deformation rates exceeding 12 cm/year (red zones) and a radial decrease in deformation intensity toward the outer regions. (B) Cumulative vertical deformation time series derived from the average of coherent pixels within the Pisciarelli area. (C) Time series of vertical deformation rate (black) and the semiannual moving average (blue).\u003c/p\u003e","description":"","filename":"Figure5.tif.jpg","url":"https://assets-eu.researchsquare.com/files/rs-7290139/v1/62af647a60cd6924c97a437c.jpg"},{"id":89345881,"identity":"73c495d6-1406-419c-a5f1-bd9a47642e5a","added_by":"auto","created_at":"2025-08-19 04:41:45","extension":"jpg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":739537,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eLAG and NO-LAG models comparative residuals analysis.\u003c/strong\u003e In (A) LAG Model and in (B) NO-LAG Model. The blue dots indicate deformation data, red points model fitting and green ones the residuals.\u003c/p\u003e","description":"","filename":"Figure6.tif.jpg","url":"https://assets-eu.researchsquare.com/files/rs-7290139/v1/df4d57b1896db651bc85b8f2.jpg"},{"id":89345887,"identity":"a568fb73-7f0a-4275-bd36-4d1d0b4e1042","added_by":"auto","created_at":"2025-08-19 04:41:45","extension":"jpg","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":1349467,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCovariate-adjusted deformation models illustrating the relationship between deformation and the others physical parameters.\u003c/strong\u003e (A) Deformation as a function of its first lag. (B) Positive correlation between seismic activity and deformation. (C) Negative correlation between heat flow and deformation. (D) Nonlinear relationship between CO₂flow and deformation. The shaded areas represent confidence intervals for the regression fits.\u003c/p\u003e","description":"","filename":"Figure7.tif.jpg","url":"https://assets-eu.researchsquare.com/files/rs-7290139/v1/d239b5312a7adab68f85b67a.jpg"},{"id":89346789,"identity":"093693b6-6a55-4310-850a-0e4b16cab8d3","added_by":"auto","created_at":"2025-08-19 05:05:45","extension":"jpg","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":1384718,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCovariate-adjusted deformation models illustrating the influence of the considered physical and geochemical parameters on ground deformation.\u003c/strong\u003e (A) Nonlinear positive relationship between seismic activity and deformation. (B) Negative correlation between temperature and deformation. (C) Positive correlation between heat flow and deformation. (D) Positive relationship between CO₂ flow and deformation. (E): Negative correlation between equilibrium pressure P(CO-CO₂-H₂-H₂O) and ground deformation. The shaded areas represent confidence intervals for the regression models.\u003c/p\u003e","description":"","filename":"Figure8.tif.jpg","url":"https://assets-eu.researchsquare.com/files/rs-7290139/v1/d8b2548870c4d7e8161ccdb1.jpg"},{"id":89345892,"identity":"258a3900-a249-49bb-97c2-9d0141eabf3d","added_by":"auto","created_at":"2025-08-19 04:41:45","extension":"jpg","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":3202947,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eIndividual time series for seismic activity, deformation, temperature, heat flow, CO₂ flux, and P(CO-CO₂-H₂-H₂O) and critical transitions.\u003c/strong\u003e (A): 2018-2022 Interval and its critical transition on the 30\u003csup\u003eth\u003c/sup\u003e of November 2020. (B) 2018-2024 Interval and its critical transition on the 1\u003csup\u003est\u003c/sup\u003e of April 2023.\u003c/p\u003e","description":"","filename":"Figure9.tif.jpg","url":"https://assets-eu.researchsquare.com/files/rs-7290139/v1/2c1b601f51722ee266713a89.jpg"},{"id":102332140,"identity":"837f4ec8-d8ba-473a-8c30-903e42bb3286","added_by":"auto","created_at":"2026-02-10 15:26:49","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":28161168,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7290139/v1/a782f4c1-ae95-44d0-9233-249b155c62e6.pdf"},{"id":89346093,"identity":"7dc80701-b736-4f0b-9ea8-036012c512db","added_by":"auto","created_at":"2025-08-19 04:49:45","extension":"docx","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":6446751,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryMaterial.docx","url":"https://assets-eu.researchsquare.com/files/rs-7290139/v1/920b6efcc177b10dcdabadbb.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Critical Transitions at Campi Flegrei resurgent caldera: A Novel Approach to Detect Early Warning Signals","fulltext":[{"header":"1. Introduction","content":"\u003cp\u003eUnderstanding how and when critical transitions occur in volcanic systems is a fundamental challenge in Earth Sciences. These systems are shaped by nonlinear and multivariable interactions between geological, geophysical, and geochemical processes that evolve dynamically over time. Anticipating changes in their behavior requires both continuous multidisciplinary monitoring and analytical frameworks capable of identifying subtle interdependencies and early warning signals. The Solfatara-Pisciarelli (SP) fumarolic area, located within the Campi Flegrei caldera (CFc) in Southern Italy, provides an exceptional natural laboratory to investigate these dynamics (Caliro et al., \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2025\u003c/span\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eThe CFc is a large resurgent caldera formed by two major explosive eruptions: the Campanian Ignimbrite (~\u0026thinsp;39 ka) and the Neapolitan Yellow Tuff (~\u0026thinsp;15 ka) (Steinmann et al., \u003cspan citationid=\"CR81\" class=\"CitationRef\"\u003e2016\u003c/span\u003e; Orsi et al., \u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2022\u003c/span\u003e) and has experienced repeated unrest episodes marked by ground deformation, seismicity, and gas anomalies. The SP area, northeast of the town of Pozzuoli, is among the most active sectors of the caldera (Smith et al., \u003cspan citationid=\"CR80\" class=\"CitationRef\"\u003e2011\u003c/span\u003e). It is characterized by intense hydrothermal activity, persistent degassing, and episodic seismic swarms, reflecting the complex coupling between deep magmatic sources and a shallow hydrothermal system (Barone et al., \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2025\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eGeophysical and geological studies have revealed that the SP region overlies a diatreme-like structure extending\u0026thinsp;~\u0026thinsp;3 km below the surface, composed of highly fractured volcanic rocks and fault systems that channel deep fluids toward the surface (Isaia et al., \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Troiano et al., \u003cspan citationid=\"CR87\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Barone et al, \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2025\u003c/span\u003e). Electrical resistivity tomography, magnetotelluric imaging, and 4D seismic tomography have confirmed the presence of fluid-saturated, high-permeability zones and documented dynamic changes in the caldera\u0026rsquo;s plumbing system during unrest phases (Siniscalchi et al., \u003cspan citationid=\"CR79\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Giacomuzzi et al., \u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e2024\u003c/span\u003e; Troiano et al., \u003cspan citationid=\"CR88\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). These structures play a critical role in controlling fluid migration and pressure redistribution, making the area particularly sensitive to magmatic or hydrothermal perturbations.\u003c/p\u003e\u003cp\u003eThe interaction between magmatic degassing and hydrothermal circulation also manifests in geochemical variations. Parameters such as CO₂ flux, CO/CO₂ and CO₂/H₂O ratios, and equilibrium temperatures and pressures of the hydrothermal system have shown temporal correlations with deformation and seismicity, offering key insights into the system\u0026rsquo;s evolution (Chiodini et al., \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e2003\u003c/span\u003e; \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). Notably, an increase in fumarolic sulfur emissions since 2018 further suggests a reorganization of fluid pathways and magmatic inputs (\u003cb\u003eCaliro et al., 2011\u003c/b\u003e). Gas fluxes reaching tens of thousands of tons per day, and the coexistence of superheated steam and liquid water in the shallow reservoir, underscore the intense energy exchange occurring at depth (Chiodini et al., \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2011\u003c/span\u003e; Bevilacqua et al., \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2024\u003c/span\u003e; Giudicepietro et al., \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2025\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eBetween 2018 and 2024, a comprehensive monitoring effort has captured high-resolution time series of key parameters, including vertical ground deformation (Tizzani et al., \u003cspan citationid=\"CR85\" class=\"CitationRef\"\u003e2024\u003c/span\u003e; Castaldo et al., \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Pepe et al., \u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; D\u0026rsquo;Auria et al., \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2015\u003c/span\u003e), seismicity (Tramelli et al., \u003cspan citationid=\"CR86\" class=\"CitationRef\"\u003e2024\u003c/span\u003e; \u003cb\u003eScotto di Uccio et al., 2024\u003c/b\u003e), surface temperature from UAV and satellite platforms (Marotta et al., \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2024\u003c/span\u003e; Mercogliano et al., \u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e2024\u003c/span\u003e), and gas compositions and fluxes. This rich multiparametric dataset enables an unique opportunity to assess how physical and chemical signals co-evolve and potentially anticipate transitions in the system.\u003c/p\u003e\u003cp\u003eIn this study, we address the need for an integrated methodology to identify and quantify critical transitions in volcanic systems. All analyses performed in the present study were carried out using several datasets acquired from this area of intense volcanic activity, as monitored by the multiparametric surveillance system of the INGV \u0026ndash; Osservatorio Vesuviano. We employ a novel framework that integrates Multivariable Fractional Polynomial Analysis (MFPA) (Royston and Sauerbrei, \u003cspan citationid=\"CR73\" class=\"CitationRef\"\u003e2008\u003c/span\u003e) and Global Critical Point Analysis (GCPA) (Rabinowitz, \u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e1986\u003c/span\u003e). MFPA enables the modeling of nonlinear relationships among variables such as ground deformation, seismicity, heat flow, and CO₂ flux, while incorporating time lags that capture delayed cause-effect interactions. Building on the MFPA outputs, GCPA applies quadratic optimization techniques to identify temporal breakpoints across multiple variables simultaneously, without relying on thresholds or the prediction of specific events.\u003c/p\u003e\u003cp\u003eOur aim is to identify systemic changes in the SP hydrothermal-magmatic complex by detecting associations between key parameters and their temporal structure. By uncovering lagged, nonlinear relationships and pinpointing critical transitions, this approach enhances our capacity to interpret volcanic unrest. Moreover, the method is scalable and adaptable to other high-risk volcanic areas, offering a promising tool for improving early warning capabilities and risk mitigation strategies.\u003c/p\u003e"},{"header":"2. Material and Methods","content":"\u003cp\u003eThe dataset analyzed in this study (Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e) spans the period from 2018 to 2024 and includes high-resolution time series of key geophysical and geochemical parameters essential for monitoring the volcanic activity of the Solfatara-Pisciarelli (SP) system. These parameters comprise vertical ground deformation (InSAR), seismicity (from the SERENADE-GOSSIP catalog), fumarolic CO₂ flux, and thermal data acquired through UAV-based infrared imaging. Additional geochemical parameters\u0026mdash;such as CO concentration, CO/CO₂ and CO₂/H₂O ratios, as well as equilibrium pressure and temperature\u0026mdash;were obtained from a combination of automatic monitoring stations and monthly laboratory analyses conducted by INGV\u0026ndash;Osservatorio Vesuviano.\u003c/p\u003e\u003cp\u003eThe monitoring focuses on the main active degassing vents: Bocca Grande (BG) and Bocca Nuova (BN) at Solfatara, and the Pisciarelli mud pool fumaroles, located less than 300 meters away on the external eastern flank of the Solfatara crater (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eB). Due to the close spatial proximity and the consistent geochemical and physical behavior observed across these sites\u0026mdash;such as similar gas compositions, temperature trends, and fumarolic dynamics\u0026mdash;these areas are treated as a unified hydrothermal-magmatic system. This interpretation is supported by previous studies (Siniscalchi et al., \u003cspan citationid=\"CR79\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Cardellini et al., \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Petrillo et al., \u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e2013\u003c/span\u003e), which highlight the common origin and interconnected behavior of the two fumarolic fields. Accordingly, in this study, the Solfatara and Pisciarelli fields are collectively referred to as the SP system.\u003c/p\u003e\u003cp\u003eTo improve the spatial characterization of the seismicity and its temporal evolution, a clustering analysis was applied to the SERENADE-GOSSIP earthquake catalog using the DBSCAN algorithm (Ester et al., \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e1996\u003c/span\u003e; Cirillo et al., \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). This density-based method identifies spatially coherent groups of seismic events and allows for the discrimination between background activity and structurally significant clusters. The analysis revealed a dominant vertical structure (Cluster 1), located beneath the Solfatara-Pisciarelli area, which is interpreted as a persistent fluid and/or fracture conduit related to the active hydrothermal system. Additional minor clusters were detected across the wider caldera, associated with transient or localized swarm-like activity. The identification of Cluster 1 was instrumental in isolating the subset of events most relevant to the SP system and in enhancing the temporal resolution of the multiparametric analysis.\u003c/p\u003e\u003cp\u003eAll time series were homogenized and normalized to facilitate cross-comparison and integration into the multivariate statistical models presented in the following sections. This section is organized in two parts: the first describes the sampling procedures and acquisition protocols for each dataset\u0026mdash;including sensor specifications, thermal clustering methods, DBSCAN parameters, and data preprocessing steps\u0026mdash;while the second presents the resulting time series for each parameter. Additional technical details, including spatial cross-sections of each seismic cluster, are provided in the Supplementary Materials (Figs. S1\u0026ndash;S6).\u003c/p\u003e\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\u003ch2\u003e2.1 Data collection\u003c/h2\u003e\u003cp\u003e\u003cem\u003eFumarolic Gas Composition (Bocca Grande - BG, Solfatara)\u003c/em\u003e\u003c/p\u003e\u003cp\u003eThe compositional data used in this study were collected monthly from Bocca Grande (BG), the main fumarolic vent at La Solfatara, using an automatic sampling station. Gas samples were analyzed for chemical composition at the Geochemistry Laboratory of INGV-Osservatorio Vesuviano. Samples were collected in under-vacuum flasks containing 4N NaOH solution (Giggenbach, \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e1975\u003c/span\u003e; Giggenbach and Gouguel, \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e1989\u003c/span\u003e) and analyzed for major gas species. Analytical procedures followed Cioni and Corazza (\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e1981\u003c/span\u003e\u003cb\u003e)\u003c/b\u003e, using gas chromatography with a single injection on dual molecular sieve columns (MS 5 \u0026Aring; capillary, 30 m \u0026times; 0.53 mm \u0026times; 50 \u0026micro;m) with TCD detectors and He/Ar as carrier gases. CO₂ absorbed in the alkaline solution was oxidized with H₂O₂ and quantified by acid\u0026ndash;base titration. CO was measured in dry gas samples using a water-cooled condenser (20\u0026ndash;30\u0026deg;C), chromatographic separation on a MS 5 \u0026Aring; 1/8 \u0026times; 50 in column (He carrier), and detection with a high-sensitivity Reduced Gas Detector (HgO).\u003c/p\u003e\u003cp\u003e\u003cem\u003eCO₂ Flux and Thermal Anomalies (Pisciarelli Site)\u003c/em\u003e\u003c/p\u003e\u003cp\u003eCO₂ emissions at the Pisciarelli fumarolic field were continuously monitored from 1 December 2018 to 30 May 2024, with daily resolution and point-based spatial data.\u003c/p\u003e\u003cp\u003eThermal surveys were carried out monthly from September 2019 to October 2023, with additional surveys up to May 2024. UAV-based radiometric thermal imaging was conducted during nighttime or evening hours to minimize solar interference. About 50 thermal mosaics were generated from flights at altitudes between 55 and 70 m a.g.l. (Marotta et al., \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2024\u003c/span\u003e), with some mosaics combining data from different flight plans. Thermal mosaics were processed using the method of Marotta et al. (\u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e2019\u003c/span\u003e) to delineate surface thermal anomalies and quantify emitted thermal energy. Clustering algorithms (DBSCAN and k-means; Ester et al., \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e1996\u003c/span\u003e; Cirillo et al., \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e2022\u003c/span\u003e) were applied to segment images into thermally homogeneous areas. Each cluster was annotated with temporal metadata and thermal attributes (e.g., average temperature, thermal flux, size in pixels and m\u0026sup2;), then reanalyzed with DBSCAN to track temporal evolution. Full details of data acquisition and processing are reported in Marotta et al. (\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2024\u003c/span\u003e)\u003c/p\u003e\u003cp\u003e\u003cem\u003eSeismicity (GOSSIP/SERENADE Catalog)\u003c/em\u003e\u003c/p\u003e\u003cp\u003eSeismic data were sourced from the SERENADE database (Peluso, \u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e2018a\u003c/span\u003e; \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2020\u003c/span\u003e) via the GOSSIP portal (Ricciolini et al., 2024), which provides hypocentral parameters, magnitude (Md), and station arrival times. Hypocenter locations were calculated using the Hypo71 code (Lee and Lahr, \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e1975\u003c/span\u003e) with local 1D velocity models for Vesuvius, Campi Flegrei, Ischia, and surrounding regions. The database integrates automatic Earthworm-based solutions (Johnson et al., \u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e1995\u003c/span\u003e) and manual revisions by the INGV-OV Monitoring Room and Seismic Laboratory. Static event pages are generated with the Serewrap software (Di Filippo and Peluso, \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2019\u003c/span\u003e), triggered with each new event insertion via the WESSEL portal (Peluso, \u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e2018b\u003c/span\u003e; \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2020\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eSeismic events from 2018 to May 2024 with Md\u0026thinsp;\u0026ge;\u0026thinsp;0.2 were extracted and spatially filtered within a 1000 m buffer around the target area.\u003c/p\u003e\u003cp\u003e\u003cem\u003eGround Deformation (Sentinel-1 InSAR Analysis)\u003c/em\u003e\u003c/p\u003e\u003cp\u003eGround deformation data were derived from Sentinel-1A (S1-A) C-band SAR images (λ\u0026thinsp;=\u0026thinsp;5.54 cm) acquired in Interferometric Wide mode (IW) from January 2018 to April 2024, along ascending (Path 44) and descending (Path 22) orbits, each with 191 images. The Small BAseline Subset (SBAS) multi-temporal DInSAR method (Berardino et al., \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2002\u003c/span\u003e; Jolivet et al., \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2011\u003c/span\u003e; Rosen et al., \u003cspan citationid=\"CR72\" class=\"CitationRef\"\u003e2000\u003c/span\u003e) was applied independently to both tracks. A total of 2101 interferograms were generated per track, with a maximum temporal baseline of 144 days. Interferograms were flattened using precise orbits and the NASA SRTM 1-arc-second DEM (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://dwtkns.com/srtm30m\u003c/span\u003e\u003cspan address=\"https://dwtkns.com/srtm30m\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). Ascending and descending LOS displacements were combined (Manzo et al., \u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e2006\u003c/span\u003e; Pepe et al., \u003cspan citationid=\"CR64\" class=\"CitationRef\"\u003e2016\u003c/span\u003e) to obtain vertical and horizontal deformation maps, geocoded and referenced to a stable point. The final dataset achieved cm-to-mm accuracy (Tizzani et al., \u003cspan citationid=\"CR85\" class=\"CitationRef\"\u003e2024\u003c/span\u003e), with a spatial resolution of ~\u0026thinsp;90 m and a\u0026thinsp;~\u0026thinsp;12-day temporal resolution, totaling 191 layers. Mean vertical deformation was computed from the time series of four pixels within 100 m of the target area.\u003c/p\u003e\u003cp\u003e\u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab1\" border=\"1\"\u003e\u003ccaption language=\"En\"\u003e\u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e\u003cdiv class=\"CaptionContent\"\u003e\u003cp\u003e\u003cb\u003eA summary table of the considered dataset is provided.\u003c/b\u003e The observation time interval and frequency of each individual measurement are reported.\u003c/p\u003e\u003c/div\u003e\u003c/caption\u003e\u003ccolgroup cols=\"3\"\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e\u003cthead\u003e\u003ctr\u003e\u003cth align=\"left\" colname=\"c1\"\u003e\u003cp\u003eData\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c2\"\u003e\u003cp\u003eTime-Interval\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colname=\"c3\"\u003e\u003cp\u003ePeriod\u003c/p\u003e\u003c/th\u003e\u003c/tr\u003e\u003c/thead\u003e\u003ctbody\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eGround Deformation\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e2018\u0026ndash;2024\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e12 Days\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eSeismic Data\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e2018\u0026ndash;2024\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eContinous\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eThermal Survey\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e2019\u0026ndash;2024\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eMontly\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eCO\u003csup\u003e2\u003c/sup\u003e Emission\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e2019\u0026ndash;2024\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eDaily\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eFumarolic Gas\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e2016\u0026ndash;2024\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eMontly\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003c/tbody\u003e\u003c/colgroup\u003e\u003c/table\u003e\u003c/div\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec4\" class=\"Section2\"\u003e\u003ch2\u003e2.2 Data description\u003c/h2\u003e\u003cp\u003eThe analysis of geochemical, thermal, seismic, and ground deformation data collected between 2018 and 2024 at the Solfatara-Pisciarelli system highlights a consistent multi-parameter response to evolving subsurface dynamics. The following sections briefly describe and interpret the observed trends.\u003c/p\u003e\u003cp\u003e\u003cem\u003eTemporal trends in fumarolic gas composition\u003c/em\u003e\u003c/p\u003e\u003cp\u003eThe SP hydrothermal system exhibits complex geochemical dynamics. The time series analysis of various geochemical parameters provided valuable insights into the system's temporal variations and underlying processes. Several key ratios were monitored, including T(CO-CO₂) [\u0026deg;C] which represents the equilibrium temperature for the hydrothermal system, CO [\u0026micro;mol/mol], CO/CO₂, CO₂/H₂O, CO₂/CH₄ ratio, and P(CO-CO₂-H₂-H₂O) [bar], which indicate the equilibrium pressure.\u003c/p\u003e\u003cp\u003eThese parameters provided insight into gas composition shifts, which are indicative of changes within the magmatic system (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eCO measurements (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eA): the compositional data used in this work were collected from Bocca Grande (BG), the main fumarolic vent of La Solfatara, using an automatic station. From 2018 to the beginning of May 2020, there was a notable increase in CO concentration. Post-May 2020, the concentration decreases rapidly until stabilizing around 2022. After 2022, the CO concentration exhibited limited fluctuations around a stable mean value. Post-2022 oscillations appear periodic and regular, while fluctuations between 2018 and 2021 were more irregular and pronounced. These observations align with studies indicating that high CO₂/CH₄ ratio in fumaroles can signal magmatic components influencing the hydrothermal system.\u003c/p\u003e\u003cp\u003eEquilibrium Temperature of the CO-CO₂ System (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eB): from 2018 to May 2020, there was a significant increase in equilibrium temperature. Post- May 2020, the temperature decreased rapidly, stabilizing around 2022. After 2022, temperatures exhibited limited fluctuations around a stable mean. Post-2022 oscillations appear periodic and regular, whereas fluctuations between 2018 and 2021 were less regular and more pronounced. The peak observed around 2021 likely represents a critical juncture for the system. Variations in equilibrium temperatures suggest underlying chemical or physical dynamics.\u003c/p\u003e\u003cp\u003eCO/CO₂ Ratio (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eC): from 2018 to May 2020, there was a significant increase in the CO/CO₂ ratio. Post-2021, the ratio decreased rapidly, stabilizing around 2022. After 2022, the CO/CO₂ ratio exhibited limited fluctuations around a stable mean, with oscillations appearing periodic and regular. Fluctuations between 2018 and 2021 were less regular and more pronounced. The peak around June 2020 likely represents a critical point for the system. Variations in the CO/CO₂ ratio suggest underlying geochemical and geophysical dynamics.\u003c/p\u003e\u003cp\u003eEquilibrium Pressure of the Hydrothermal System P(CO-CO\u003csub\u003e2\u003c/sub\u003e-H\u003csub\u003e2\u003c/sub\u003e-H\u003csub\u003e2\u003c/sub\u003eO) (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eD): from 2018 to 2021, there was a notable increase in equilibrium pressure. Post-2021, the pressure decreased rapidly, stabilizing around 2022. After 2022, equilibrium pressure showed limited fluctuations around a stable mean, with oscillations appearing periodic and regular. Fluctuations between 2018 and 2021 were less regular and more pronounced. The peak around 2021 likely signifies a critical point for the system. Variations in equilibrium pressure suggest underlying geochemical and geophysical dynamics.\u003c/p\u003e\u003cp\u003eCO₂/H₂O Ratio (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eE): from 2018 to 2021, there was a significant increase in the CO₂/H₂O ratio. Post-2021, the ratio decreased rapidly, stabilizing around 2022. After 2022, the CO₂/H₂O ratio exhibited limited fluctuations around a stable mean, with oscillations appearing periodic and regular. Fluctuations between 2018 and 2021 were less regular and more pronounced. The peak around 2021 likely represents a critical point for the system. Variations in the CO₂/H₂O ratio suggest underlying geochemical and geophysical dynamics.\u003c/p\u003e\u003cp\u003eCO₂/CH₄ Ratio (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eF): from 2018 to April 2020, there was a significant increase in the CO₂/CH₄ ratio. Post-2021, the ratio increased rapidly, stabilizing around 2022. After 2022, the CO₂/CH₄ ratio exhibited limited fluctuations around a stable mean, with oscillations appearing periodic and regular.\u003c/p\u003e\u003cp\u003e\u003cem\u003eVariations in CO2 fluxes at Pisciarelli and their implications\u003c/em\u003e\u003c/p\u003e\u003cp\u003eCO₂ Flux measurements (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eG): the Pisciarelli fumarolic site has been continuously monitored for CO₂ emissions over a data interval spanning from 1 December 2018 to 30 May 2024. Figure\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eG presents the time series analyzed as a monthly moving average, while the semiannual average (blue curve) highlights long-term trends and patterns. The time series reveals multiple peaks, most notably around mid-2020, late 2021, and early 2023, where the CO₂ flux exceeded 3.5 \u0026times; 10⁴ g/m\u0026sup2;/day. These peaks coincide with known periods of increased volcanic activity in the CFc. Significant reductions in CO₂ emissions were observed, with the lowest flow values occurring around late 2020 and late 2022, suggesting potential periods of reduced degassing or transient system stability. The semiannual average smooths out daily fluctuations revealing recurring cyclic patterns which may indicate a potential link between seasonal changes and degassing activity. Notably, periods of heightened activity are interspersed with phases of reduced emissions, underscoring the cyclic nature of the hydrothermal system's behavior. From 2020 onward, the amplitude of fluctuations increased, reflecting a potential intensification of subsurface processes, such as magmatic degassing or variations in hydrothermal reservoir pressurization. Over the monitoring period, the CO₂ flux shows a gradual increase in baseline values, which may signal ongoing pressurization within the magmatic-hydrothermal system or evolving permeability pathways. The detailed analysis of the CO₂ flux at Pisciarelli highlights both short-term variability and long-term trends. Peaks in emissions align with known episodes of unrest, while the semiannual average reveals broader patterns that provide crucial insights into the behavior of the volcanic system. In summary, the SP hydrothermal system has exhibited significant variations across multiple geochemical parameters between 2018 and 2024. The period from 2018 to 2021 is characterized by increasing trends, reaching peaks around 2021, followed by rapid declines and stabilization around 2022. Post-2022, the system shows regular, periodic fluctuations, suggesting the establishment of a new equilibrium state may be influenced by cyclic external factors.\u003c/p\u003e\u003cp\u003e\u003cem\u003eSpatial and temporal evolution of surface thermal anomalies\u003c/em\u003e\u003c/p\u003e\u003cp\u003eThe thermal data used in this study extend beyond those presented in Marotta et al. (\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2024\u003c/span\u003e\u003cb\u003e)\u003c/b\u003e. Figure\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e integrates spatial and temporal analyses, illustrating the heating dynamics through land surface temperature (LST) and statistical trends. Specifically, Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eA shows mean and annual thermal maps from 2020 to 2024. In Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eB, the LST time-series plot illustrates measured temperatures (black), de-seasonalized temperatures (red), and seasonal trends (green) are reported. The data reveal a pronounced seasonal cycle (green) superimposed on broader temperature variations. De-seasonalized temperatures (red) isolate long-term thermal anomalies, providing insights into subsurface processes. Figure\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eC focuses on de-seasonalized temperatures (red) alongside their six-month moving average (black). The smoothed trendline offers a clearer visualization of gradual temperature changes, indicating potential long-term geothermal or hydrothermal dynamics in the study area. Figure\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD shows the heat flux time series (blue) and the heat release rates (red). The consistent hotspots in mean map (red/yellow zones in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eA) suggest persistent thermal anomaly activity, providing a comprehensive view of the thermal behavior over the analyzed period. Specifically, the mean annual map highlights two dominant patterns of thermal activity within the studied area. The first pattern, located on the right side of the maps, corresponds to the Pisciarelli mud pool area, exhibiting the highest thermal anomalies, with intense red zones persisting across all years, likely due to concentrated hydrothermal venting and localized geothermal heat flux. The second notable anomaly pattern is observed along the external flank of the Solfatara volcano, specifically near Via Scarfoglio area on the left side of the maps. This region consistently displays medium-to-low thermal anomalies (green-to-yellow zones), indicative of diffuse heat flow and potential subsurface thermal activity extending from the Solfatara hydrothermal system. Summarizing the observations: (i) Pisciarelli anomaly is a persistent and intense hotspot with slight interannual variability in intensity and spatial extent, representing the core of hydrothermal venting and geothermal heat emission. (ii) Via Scarfoglio Pattern consists of a diffuse and moderate heat flow along the external flank of the Solfatara volcano, showing gradual spatial expansion over time, potentially linked to subsurface dynamics. (iii) The temporal evolution of the LST field emphasizes that year 2022 is remarkable for an increased thermal activity, particularly in the Pisciarelli zone (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eC) even if the overall spatial distribution of thermal anomalies remains consistent. Accordingly, a computed heat flow time series derived from the UAV-based thermal dataset from September 2019 to October 2023 is reported in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD, with a variable temporal resolution (minimum value 1 month). The estimate considers between 5 to 20 points per temporal instance, with a six-month moving average applied. Analysis of the heat flux and its variations: Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD presents the time series of heat flux (blue) and heat flux variations (red), providing insights into the geothermal dynamics of the system. The heat flux series exhibits a relatively smooth trend interspersed with significant fluctuations. From 2020 to early 2022, a gradual increase is observed, indicating a progressive intensification of geothermal activity. This trend stabilizes around mid-2022, followed by a slight decline extending through late 2023. The six-month moving average highlights long-term patterns by filtering short-term noise and emphasizing overall heat flux trajectory in the region. A comparative analysis between heat flux and its variations (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD) reveals important correlations. The most prominent peak in mid-2022, is clearly evident in both heat flux (blue) and its variations (red), suggesting a system-wide thermal perturbation. However, while the heat flux series displays a smoother evolution, the variations series exhibits pronounced spikes, indicative of episodic fluctuations in the geothermal regime. Heat flux variations, represented on the right axis of Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD, highlight relative intensity changes over time. The series of variations maintains a generally oscillatory pattern, with fluctuations contained within a predictable range. This consistency suggesting a dominant seasonal component likely influenced by atmospheric conditions and external factors affecting surface thermal measurements. In contrast, significant spikes in heat flux, particularly in mid-2020 and mid-2022, indicate abrupt events potentially linked to fluid injections, pressure changes, or localized thermal surges. These events temporarily disrupt the otherwise stable trajectory of the system. The synchronization of anomalies in both series, particularly around mid-2022, underscores a system-wide geothermal perturbation, likely driven by subsurface dynamics. By examining heat flux series as the central trend and analyzing deviations through the variations series, it becomes evident that the Pisciarelli geothermal system exhibits both stable long-term trends and episodic fluctuations.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003cem\u003eSeismic activity and clustering beneath the hydrothermal area\u003c/em\u003e\u003c/p\u003e\u003cp\u003eIn the present study, the CF seismic catalog was analyzed using the DBSCAN algorithm to identify and characterize spatial clusters of seismic activity (Supplementary Text). This method, particularly suited for datasets with varying density distributions, enables the detection of regions where seismic events are densely concentrated while distinguishing isolated events considered as noise. The results of this analysis provide a refined picture of the seismic dynamics within the caldera, offering critical insights into the spatial organization of earthquakes and their potential connection to underlying geological structures. The outcomes of the clustering analysis are illustrated in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e (see also Figs. S1-S6), where the spatial distribution of seismic clusters is represented through contour lines depicting event density. The plan view (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eA) offers an overview of how seismicity is distributed across the Campi Flegrei area, with different clusters identified by distinct colors. The contour lines effectively highlight the areas with the highest concentration of earthquakes, delineating regions of intensified seismic activity. To further investigate the depth distribution of these clusters, two vertical cross-sections were considered. The North-South cross-section (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB) provides a view of how seismicity extends in depth along this axis, revealing well-defined zones of high activity. Similarly, the East-West cross-section (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eC) confirms the presence of deep-seated clusters and highlights their spatial correlation with shallow seismicity, suggesting a structured arrangement of events that may be influenced by fault networks and magmatic or hydrothermal processes. Among the various clusters identified, particular attention is given to Cluster 1, located directly beneath the SP area. The seismicity associated with this cluster presents a highly concentrated pattern, both laterally and in depth, and is likely influenced by subsurface fluid movements and stress accumulation within the volcanic system. Given its strategic location and relevance, this study will focus specifically on the seismic events belonging to Cluster 1, aiming to investigate their temporal evolution and potential implications for the dynamic behavior of the SP hydrothermal system. The spatial clustering results obtained through DBSCAN not only confirm the presence of localized zones of high seismicity but also provide essential clues regarding the interaction between tectonic and magmatic processes within Campi Flegrei. The density contours reveal that earthquake distribution is not random but follows specific structural trends, reinforcing the hypothesis that seismicity in the area is closely linked to the caldera\u0026rsquo;s ongoing geodynamic evolution. The identification and analysis of these clusters contribute to a deeper understanding of the mechanisms driving volcanic unrest, such as fluid migration and/or fault activities.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eFigure \u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eD provides a representative time series of earthquakes per month for Cluster 1 (red in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eA-C) over the period from 2018 to 2024, with normalized density values on the \u003cem\u003ey\u003c/em\u003e-axis and time in years on the \u003cem\u003ex\u003c/em\u003e-axis.\u003c/p\u003e\u003cp\u003eFrom 2018 to the end of 2020, seismicity density remained relatively low and stable, indicating minimal seismic activity in the region. A gradual increase became evident starting in 2020, suggesting a progressive intensification of the underlying geological processes. This trend continues, culminating in a significant rise between 2022 and 2023, marking a period of pronounced seismic activity. In 2023, the seismicity density reached its peak, followed by fluctuations in early 2024. These variations may be attributed to changes in underground pressure, fluid migration, or small-scale eruptive dynamics. Toward the latest part of the dataset in 2024, a slight decline in seismicity density is observed, potentially indicating a temporary stabilization or a reduction in subsurface activity.\u003c/p\u003e\u003cp\u003eOverall, the data highlight a phase of increased seismic activity in SP, which may be linked to ongoing volcanic or geothermal processes. A more detailed analysis integrating geophysical and geochemical parameters would be necessary to further elucidate the mechanisms driving these observed seismic trends.\u003c/p\u003e\u003cp\u003e\u003cem\u003eGround deformation as an indicator of subsurface processes\u003c/em\u003e\u003c/p\u003e\u003cp\u003eThe vertical ground deformation at the Pisciarelli area, derived from InSAR data and based on the average of three coherent SAR pixels within the region of interest, exhibits notable spatial and temporal characteristics. Spatially, the mean ground velocity map (Fig.\u0026nbsp;5A) highlights a well-defined uplift area, with maximum velocities exceeding 12 cm/year. The deformation pattern decreases radially from the central uplift, forming a gradient consistent with subsurface processes such as hydrothermal activity. The vertical component was then used to compute the monthly deformation rate, which was smoothed using a six-month moving average. Temporally, the cumulative deformation time series (Fig.\u0026nbsp;5B) shows a continuous and steady increase from 2018 to 2024, reaching a total displacement of approximately 45 cm by the end of the observation period (for the pixels in Pisciarelli area). The near-linear trend in cumulative deformation indicates persistent driving forces in the system, likely related to pressurization or volumetric expansion within the subsurface structure. The time series of vertical deformation rates (Fig.\u0026nbsp;5C) provides a more detailed view of the temporal variability. The rates exhibit oscillatory behavior, with periods of acceleration and deceleration. Peaks in deformation rate occur around 2019, 2021, and late 2023, suggesting episodic intensifications in the processes driving deformation. The semiannual moving average captures overall trends while smoothing short-term fluctuations, showing a general increase in deformation activity punctuated by temporary slowdowns. The dataset, collected with high temporal resolution (~\u0026thinsp;12 days) and spatial resolution (~\u0026thinsp;90 m), effectively captures both the localized nature and dynamic evolution of the deformation. The observed patterns suggest a strong connection between subsurface processes and surface deformation, emphasizing the importance of continuous InSAR monitoring for understanding the ongoing dynamics of the Pisciarelli area.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec5\" class=\"Section2\"\u003e\u003ch2\u003e2.3 Data Processing and methods\u003c/h2\u003e\u003cp\u003eThe above-mentioned dataset were initially processed to allow comparisons across varying units and scales. In particular, normalization and temporal smoothing techniques (e.g., moving averages), have been applied to highlight significant trends.\u003c/p\u003e\u003cdiv id=\"Sec6\" class=\"Section3\"\u003e\u003ch2\u003e2.3.1 Multivariable Fractional Polynomial Analysis (MFPA)\u003c/h2\u003e\u003cp\u003eThe primary objective of this study was to investigate the association between the temporal behavior of deformation, seismic activity, and specific physical and chemical factors. This goal, along with the characteristics of the available dataset, guided the selection of specific statistical procedures and analytical approaches.\u003c/p\u003e\u003cp\u003eA multivariable regression model was employed to assess the association between our main variable of interest (vertical ground deformation) and the set of measured physical parameters. The aim is to identify statistically significant correlation for each variable, independent of other factors.\u003c/p\u003e\u003cp\u003eSince the data are structured as a single-point time series across multiple variables, an autoregressive model was adopted to mitigate statistical bias arising from the dependent variable\u0026rsquo;s non-independence, which can occur in regression models for this type of data. To address this, LAG value of the dependent factor was included as an additional independent variable in the model. The time lag order, i.e. the number of periods going back in time, was set to 1 unit (i.e., 1/10 of a year) based on the following considerations: (i) while increasing the LAG order may improve forecasting accuracy, it also increases model complexity, reduces interpretability, and raises the risk of multicollinearity. Since the goal of this work is not to develop a forecasting model, we keep this parameter low. (ii) A first-LAG model sufficiently captures the temporal dependence of the response on the dependent variable, as each value directly influences the subsequent one (Finkel, \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2007\u003c/span\u003e), thereby addressing the temporal dependence bias. Additionally, by acting as a proxy for unknown factors that affect the dependent variable, the first LAG absorbs unmeasured effects, helping to stabilize estimates of the independent effects of other covariates and enabling more accurate evaluation. To assess the coherence of the relationships identified by each model, comparisons were made between models with and without the autoregressive (LAG) variable.\u003c/p\u003e\u003cp\u003eThe identification of the best final model followed the pragmatic approach suggested by Royston and Sauerbrei (\u003cspan citationid=\"CR73\" class=\"CitationRef\"\u003e2008\u003c/span\u003e\u003cb\u003e)\u003c/b\u003e, summarized as follows:\u003c/p\u003e\u003cp\u003e\u003col\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eSelection of the final model by including factors significantly associated with the dependent variable, while assessing the functional form (linearity/non-linearity) of continuous factors using the multivariable fractional polynomial (MFP) algorithm.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eEvaluation of the overall and individual contribution of the factors selected in the final model to the dependent variable by calculating both the global variance explained (R\u0026sup2;) and the individual variance partition for each significant factor (R\u0026sup2;p). The Shapley-Owen decomposition algorithm (Shorrocks, \u003cspan citationid=\"CR76\" class=\"CitationRef\"\u003e2013\u003c/span\u003e) was used to partition global R\u0026sup2;, quantifying each variable's impact independently of others, i.e., the impact with all other factors held constant.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003c/ol\u003e\u003c/p\u003e\u003cp\u003eAll models were checked for multicollinearity by calculating the Variance Inflation Factor (VIF) for each variable in the final model. A VIF value greater than 5 was considered indicative of multicollinearity (Kutner et al., \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2005\u003c/span\u003e; Suleiman, \u003cspan citationid=\"CR82\" class=\"CitationRef\"\u003e2015\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eTo visualize the relationship between a dependent variable and a specific independent factor controlling for other significant variables in the final model, a scatterplot was created. In this plot, the dependent variable values were adjusted with respect to the mean of other significant variables considered in the model. Specifically, the relationship between \u003cem\u003ey\u003c/em\u003e and one of the \u003cem\u003ex\u003c/em\u003e variables was illustrated by adjusting observed \u003cem\u003ey\u003c/em\u003e values for each item according to deviation of other observed \u003cem\u003ex\u003c/em\u003e variables from their respective means, using estimated coefficients. Data were analyzed using \u003cem\u003eStata version 18.0\u003c/em\u003e (StataCorp LP, College Station, TX, USA). Statistical significance was set at p\u0026thinsp;\u0026le;\u0026thinsp;0.05 for univariate comparisons, variable selection, and testing between variable transformations within the MFP procedure.\u003c/p\u003e\u003cp\u003eThe factors (i.e., physical parameters) significantly associated in the performed model are used as input for a Global Critical Point Analysis to identify meaningful points in the time domain following described. The MFPA model converged in fewer than 30 iterations per run, with an average computation time below 1 minute on standard desktop hardware (Intel i7 processor, 16GB RAM). This performance confirms the suitability of the approach for potential near-real-time implementations in operational monitoring environments.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec7\" class=\"Section3\"\u003e\u003ch2\u003e2.3.2 Global Critical Point Analysis (GCPA)\u003c/h2\u003e\u003cp\u003eGCPA method is designed to identify significant temporal transitions in multiparametric historical datasets. It leverages quadratic programming (QP) (\u003cb\u003e57\u003c/b\u003e) to optimize the detection of a critical point, that represents the time when the combined dynamics of multiple variables exhibit significant changes. This critical point may correspond to a phase where the system transitions from stability to chaos or shifts between different equilibrium states. In the context of the SP system, GCPA serves as a valuable tool for analyzing variations in geophysical and geochemical parameters, helping to detect potential systemic changes without aiming to predict specific events, such as an eruption. Instead, the methodology focuses on understanding the system's evolution and identifying precursors to critical transitions, essential in volcanic systems where gradual or subtle changes may indicate an approaching tipping point (global critical point). The GCPA statistical technique identifies the temporal moment when significant changes occur across multiple variables simultaneously. These variables include physical or chemical parameters recorded over time, such as ground deformation, seismicity, temperature, heat flow, and CO₂ flux. By pinpointing a single temporal moment representing overall change across all analyzed time series, GCPA simplifies the complexity of multiparametric analysis and identifies a common breaking point explaining observed variations. Moreover, the application of quadratic optimization using the Karmarkar method (Karmarkar, \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e1984\u003c/span\u003e) has improved the efficiency of global critical point detection, allowing for a more precise assessment of systemic changes. This optimization framework enhances the reliability of transition detection by providing a mathematical basis for identifying system shifts beyond simple threshold-based analyses. The analysis was performed using R software (R Core Team, \u003cspan citationid=\"CR70\" class=\"CitationRef\"\u003e2014\u003c/span\u003e) and involved five main steps:\u003c/p\u003e\u003cp\u003e\u003col\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eData Preprocessing and Normalization: input time series are normalized to standardize datasets with different units or scales, (e.g., temperature, seismicity, gas flux), ensuring meaningful and robust comparison of variables while reducing the risk of bias in subsequent analysis.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eQuadratic Programming for Critical Point Identification: the problem is modeled using a quadratic objective function that maximizes deviations from expected or average behavior. The identity matrix Q represents the quadratic structure of the optimization problem, while the vector c captures normalized absolute deviations. Temporal constraints are applied to ensure that the identified critical point lies within the observed data range, preventing the selection of spurious points.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003c/ol\u003e\u003cdiv class=\"BlockQuote\"\u003e\u003cp\u003eTo formally identify the time of critical transition, we structured the analysis as a quadratic programming (QP) optimization problem. Specifically, the critical point was defined as the time index minimizing the following objective function:\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$\\:\\left(\\raisebox{1ex}{$1$}\\!\\left/\\:\\!\\raisebox{-1ex}{$2$}\\right.\\right)\\bullet\\:{\\varvec{x}}^{\\varvec{T}}\\bullet\\:\\varvec{Q}\\bullet\\:\\varvec{x}+{\\varvec{c}}^{\\varvec{T}}\\bullet\\:\\varvec{x}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003e\u003cdiv class=\"BlockQuote\"\u003e\u003cp\u003ewhich represents a quadratic function in matrix form, commonly used in quadratic optimization (Quadratic Programming, QP); where \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\varvec{x}\\in\\:{\\varvec{R}}^{\\varvec{n}}\\)\u003c/span\u003e\u003c/span\u003e is a binary indicator vector (with 1 at the candidate critical time and 0 elsewhere), \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\varvec{Q}\\in\\:{\\varvec{R}}^{\\varvec{n}\\times\\:\\varvec{n}}\\)\u003c/span\u003e\u003c/span\u003e is a symmetric matrix representing the covariance structure between the standardized input variables (often also defined as positive semi-definite, to ensure convexity), and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\varvec{c}\\in\\:\\varvec{R}\\)\u003c/span\u003e\u003c/span\u003e is a column vector of linear coefficients collecting the normalized mean deviation of each parameter over time. More specifically, the function:\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Equb\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equb\" name=\"EquationSource\"\u003e\n$$\\:\\varvec{f}\\left(\\varvec{x}\\right)=\\left(\\frac{1}{2}\\right)\\bullet\\:{\\varvec{x}}^{\\varvec{T}}\\bullet\\:\\varvec{Q}\\bullet\\:\\varvec{x}+{\\varvec{c}}^{\\varvec{T}}\\bullet\\:\\varvec{x}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003e\u003cdiv class=\"BlockQuote\"\u003e\u003cp\u003eis a typical objective function in quadratic optimization problems, where:\u003c/p\u003e\u003cp\u003e\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\left(\\frac{1}{2}\\right)\\bullet\\:{\\varvec{x}}^{\\varvec{T}}\\bullet\\:\\varvec{Q}\\bullet\\:\\varvec{x}\\)\u003c/span\u003e\u003c/span\u003e, is the quadratic term, and represents a paraboloid if QQ is positive definite;\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Equc\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equc\" name=\"EquationSource\"\u003e\n$$\\:{\\varvec{c}}^{\\varvec{T}}\\bullet\\:\\varvec{x},\\:\\text{i}\\text{s}\\:\\text{t}\\text{h}\\text{e}\\:\\text{l}\\text{i}\\text{n}\\text{e}\\text{a}\\text{r}\\:\\text{t}\\text{e}\\text{r}\\text{m}.$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003e\u003cdiv class=\"BlockQuote\"\u003e\u003cp\u003eNote that the factor \u003cem\u003e1/2\u003c/em\u003e is conventional and useful because the derivative of the quadratic term is \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\varvec{Q}\\bullet\\:\\varvec{x}\\)\u003c/span\u003e\u003c/span\u003e which simplifies the gradient computation:\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Equd\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equd\" name=\"EquationSource\"\u003e\n$$\\:\\nabla\\:\\varvec{f}\\left(\\varvec{x}\\right)=\\varvec{Q}\\bullet\\:\\varvec{x}+\\varvec{c}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003e\u003cdiv class=\"BlockQuote\"\u003e\u003cp\u003eFinally, we remark that the matrix \u003cb\u003eQ\u003c/b\u003e encodes the co-evolution and coupling between variables (e.g., deformation, seismicity, gas), while \u003cb\u003ec\u003c/b\u003e provides a temporal profile of their joint anomalies. This formulation allows us to detect the time step where the system collectively exhibits the greatest departure from background behavior.\u003c/p\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003e\u003col start=\"3\"\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eOptimization Techniques: the R Optimization Infrastructure (ROI) package (Theu\u0026szlig;l et al., \u003cspan citationid=\"CR84\" class=\"CitationRef\"\u003e2020\u003c/span\u003e), combined with the \u003cem\u003eROI.plugin.quadprog plugin\u003c/em\u003e, is used to solve the quadratic programming problem. This setup efficiently handles large datasets with multiple parameters and integrates advanced optimization techniques, including:\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003c/ol\u003e\u003c/p\u003e\u003cp\u003e\u003cul\u003e\u003cli\u003e\u003cp\u003e\u003cem\u003eInterior-Point Methods\u003c/em\u003e: ensure convergence by iteratively refining the solution within the feasible region.\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003e\u003cem\u003eConjugate Gradient Methods\u003c/em\u003e: particularly suitable for large and sparse datasets.\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003e\u003cem\u003ePenalty and Barrier Methods\u003c/em\u003e: facilitate constraint handling and enhance computational efficiency.\u003c/p\u003e\u003c/li\u003e\u003c/ul\u003e\u003c/p\u003e\u003cp\u003e\u003col start=\"4\"\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eCritical Point Detection: the index corresponding to the maximum deviation within the normalized time series is identified as the global critical point, representing the moment when significant transitions occur across the analyzed parameters (Boyd and Vandenberghe, \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2004\u003c/span\u003e).\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eVisualization and Interpretation: the identified critical point is visualized within the context of normalized and smoothed time series using the \u003cem\u003eggplot2 package\u003c/em\u003e (Wickham, \u003cspan citationid=\"CR90\" class=\"CitationRef\"\u003e2009\u003c/span\u003e) for trend enhancement. Additional visualizations compare raw, normalized, and smoothed data to provide a comprehensive view of the dataset's temporal evolution (Shumway and Stoffer, \u003cspan citationid=\"CR77\" class=\"CitationRef\"\u003e2017\u003c/span\u003e).\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003c/ol\u003e\u003c/p\u003e\u003cp\u003eThe GCPA method is particularly suited to complex systems with dynamic and interdependent parameters. By integrating significant geophysical and geochemical variables, such as ground deformation, seismicity, temperature, heat flow, and CO₂ flux, the algorithm identifies moments of significant systemic change. These insights are crucial for understanding volcanic behavior, assessing hazards, and refining monitoring strategies. The GCPA framework is adaptable to large datasets and diverse scientific applications, allowing the incorporation of various constraints and objective functions tailored to specific research needs. For our analysis, we use the time-series of physical parameters retrieved by the MFPA. To evaluate the method's response and its robustness to variations in multiparametric time series, we performed two analyses by considering different observation intervals for the parameters: the first (2018\u0026ndash;2023) that represents a subset of the second (2018\u0026ndash;2024).\u003c/p\u003e\u003c/div\u003e\u003c/div\u003e"},{"header":"3. Results","content":"\u003cdiv id=\"Sec9\" class=\"Section2\"\u003e\u003ch2\u003e3.1 Best-fitting multivariable model for ground deformation using time-lagged analysis\u003c/h2\u003e\u003cp\u003eWe applied MFPA to determine the best combination of geophysical and geochemical parameters that explain ground deformation at the SP site, considered as the dependent variable. This choice is motivated by the prominent role that vertical ground deformation plays in signaling pressurization processes at the Campi Flegrei caldera. Long-term geodetic observations have consistently linked episodes of uplift to periods of increased hydrothermal or magmatic activity (Del Gaudio et al., \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2010\u003c/span\u003e; Tizzani et al., \u003cspan citationid=\"CR85\" class=\"CitationRef\"\u003e2024\u003c/span\u003e; Vanorio et al., \u003cspan citationid=\"CR89\" class=\"CitationRef\"\u003e2025\u003c/span\u003e). In the Solfatara-Pisciarelli area, deformation is often the first measurable indicator of subsurface changes, and it integrates complex interactions such as gas accumulation, heat transfer, and fluid migration. Therefore, it provides an effective dependent variable for exploring multiparametric cause-effect relationships in the SP system. In addition, deformation is measured with high spatial and temporal resolution via InSAR techniques, ensuring the reliability and continuity required for regression-based analyses. This makes it a valuable proxy for unrest modeling, particularly when combined with complementary parameters such as CO₂ flux, seismicity, and heat flow.\u003c/p\u003e\u003cp\u003eThe goal was to capture significant nonlinear associations and delayed effects that influence vertical deformation dynamics, offering a robust interpretation of the physical mechanisms behind observed surface displacements. Following the assessement of multicollinearity among all parameters (geochemical parameters: T(CO-CO₂) [\u0026deg;C], CO [\u0026micro;mol/mol], CO/CO₂, CO₂/H₂O, CO₂/CH₄, P(CO-CO₂-H₂-H₂O) [bar] and CO₂ flux; land surface temperature; heat flux and seismicity), only those showing no collinearity were included in the MFPA. Two models were compared: one that incorporates a 30-day time lag (LAG model), and one that does not (NO-LAG model).\u003c/p\u003e\u003cp\u003eTo enhance the interpretability of temporal dependencies among monitored parameters, we developed two regression models: one incorporating a time-lag structure (LAG model) and one without (NO-LAG model). In the LAG model, the dependent variable (vertical ground deformation) is regressed against its own past value (i.e., the value observed 30 days earlier), in addition to other explanatory variables. This autoregressive term captures delayed cause-effect dynamics that are common in volcanic systems, where processes such as pressurization, degassing, and heat transfer propagate over time through subsurface pathways.\u003c/p\u003e\u003cp\u003eThe rationale behind the 30-day lag was based on:\u003c/p\u003e\u003cp\u003e\u003col style=\"list-style-type:lower-roman;\"\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003ethe average response time of deformation to known geochemical variations at CFc as inferred from empirical monitoring (Chiodini et al., \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2017\u003c/span\u003e);\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003ethe need to introduce temporal structure into the model without overfitting or increasing multicollinearity;\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003ethe goal of capturing the immediate memory of the system while ensuring robustness.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003c/ol\u003e\u003c/p\u003e\u003cp\u003eThis approach aligns with previous studies modeling dynamic interactions in volcanic environments, where LAGs in the order of weeks to months are commonly used to represent internal feedbacks and delayed surface responses.\u003c/p\u003e\u003cp\u003eFrom a statistical perspective, the first-LAG structure offers a balance between capturing temporal autocorrelation and avoiding overfitting or collinearity issues associated with higher-order lags. Testing multiple lag values would increase model complexity and reduce interpretability, without providing substantial gain in explained variance for the current scope of identifying critical transitions (rather than forecasting). Nonetheless, future developments will explore multi-lag configurations and rolling-window approaches to further generalize the framework.\u003c/p\u003e\u003cp\u003e\u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab2\" border=\"1\"\u003e\u003ccaption language=\"En\"\u003e\u003cdiv class=\"CaptionNumber\"\u003eTable 2\u003c/div\u003e\u003cdiv class=\"CaptionContent\"\u003e\u003cp\u003esummarizes the coefficients, explained variance, and model selection metrics for both approaches. In the LAG model (Panel A), we found that four predictors significantly contribute to deformation: the first lag of deformation, seismic activity, heat flow, and CO₂ flux. This model explained 99.9% of the total variance (global R\u0026sup2; = 0.999), demonstrating that it provides the optimal multivariable configuration for this dataset. The first-lag term alone explained 59% of the variance, followed by seismic activity (33.1%), heat flow (4.9%), and CO₂ flux (3.0%). The Akaike Information Criterion (AIC = -12.3) and Bayesian Information Criterion (BIC = -4.9) further supported the model\u0026rsquo;s.\u003c/p\u003e\u003c/div\u003e\u003c/caption\u003e\u003ccolgroup cols=\"11\"\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c5\" colnum=\"5\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c6\" colnum=\"6\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c7\" colnum=\"7\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c8\" colnum=\"8\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c9\" colnum=\"9\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c10\" colnum=\"10\"\u003e\u003c/div\u003e\u003cdiv align=\"left\" class=\"colspec\" colname=\"c11\" colnum=\"11\"\u003e\u003c/div\u003e\u003cthead\u003e\u003ctr\u003e\u003cth align=\"left\" colname=\"c1\"\u003e\u0026nbsp;\u003c/th\u003e\u003cth align=\"left\" colspan=\"5\" nameend=\"c6\" namest=\"c2\"\u003e\u003cp\u003ePanel A\u003c/p\u003e\u003c/th\u003e\u003cth align=\"left\" colspan=\"5\" nameend=\"c11\" namest=\"c7\"\u003e\u003cp\u003ePanel B\u003c/p\u003e\u003c/th\u003e\u003c/tr\u003e\u003c/thead\u003e\u003ctbody\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colspan=\"5\" nameend=\"c6\" namest=\"c2\"\u003e\u003cp\u003eGlobal R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.999\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colspan=\"5\" nameend=\"c11\" namest=\"c7\"\u003e\u003cp\u003eGlobal R\u003csup\u003e2\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.97\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u0026nbsp;\u003c/td\u003e\u003ctd align=\"left\" colspan=\"5\" nameend=\"c6\" namest=\"c2\"\u003e\u003cp\u003eAIC= -12.3; BIC=-4.9\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colspan=\"5\" nameend=\"c11\" namest=\"c7\"\u003e\u003cp\u003eAIC\u0026thinsp;=\u0026thinsp;114.9; BIC\u0026thinsp;=\u0026thinsp;123.8\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eParameter\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eFunct. Form\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003eCoeff.\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003ep\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003eR\u003csup\u003e2\u003c/sup\u003ep (%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003eBIF (%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003eFunct. Form\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003eCoeff.\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003ep\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c10\"\u003e\u003cp\u003eR\u003csup\u003e2\u003c/sup\u003ep (%)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c11\"\u003e\u003cp\u003eBIF (%)\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cem\u003eFirst Lag\u003c/em\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e\u003cb\u003eLin.\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e\u003cb\u003e0,98\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e\u003cb\u003e\u0026le;\u0026thinsp;0.001\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e\u003cb\u003e59\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e\u003cb\u003e100\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e\u003cb\u003e-\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e\u003cb\u003e-\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e\u003cb\u003e-\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c10\"\u003e\u003cp\u003e\u003cb\u003e-\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c11\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eSeismicity\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e\u003cb\u003eLin.\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e\u003cb\u003e9,3\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e\u003cb\u003e\u0026le;\u0026thinsp;0.001\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e\u003cb\u003e33,1\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e\u003cb\u003e97,4\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e\u003cb\u003e1/x\u003c/b\u003e\u003csup\u003e\u003cb\u003e2\u003c/b\u003e\u003c/sup\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e\u003cb\u003e-0,02\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e\u003cb\u003e\u0026le;\u0026thinsp;0.001\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c10\"\u003e\u003cp\u003e\u003cb\u003e41,8\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c11\"\u003e\u003cp\u003e\u003cb\u003e99,7\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eTemperature\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eLin.\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e0,005\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e=\u0026thinsp;0.58\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e19,2\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e\u003cb\u003eLin.\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e\u003cb\u003e-0,2\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e\u003cb\u003e\u0026le;\u0026thinsp;0.001\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c10\"\u003e\u003cp\u003e\u003cb\u003e1,0\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c11\"\u003e\u003cp\u003e\u003cb\u003e84,7\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eHeat Flow\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e\u003cb\u003eLin.\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e\u003cb\u003e0,00023\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e\u003cb\u003e=\u0026thinsp;0.003\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e\u003cb\u003e4,9\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e\u003cb\u003e24,4\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e\u003cb\u003ex\u003c/b\u003e\u003csup\u003e\u003cb\u003e2\u003c/b\u003e\u003c/sup\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e\u003cb\u003e2,7\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e\u003cb\u003e\u0026le;\u0026thinsp;0.001\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c10\"\u003e\u003cp\u003e\u003cb\u003e12,7\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c11\"\u003e\u003cp\u003e\u003cb\u003e99,7\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eCO\u003csub\u003e2\u003c/sub\u003e Flow\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e\u003cb\u003ex\u003c/b\u003e\u003csup\u003e\u003cb\u003e3\u003c/b\u003e\u003c/sup\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e\u003cb\u003e-0,01\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e\u003cb\u003e\u0026le;\u0026thinsp;0.001\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e\u003cb\u003e3,0\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e\u003cb\u003e65,9\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e\u003cb\u003eLin.\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e\u003cb\u003e0,0002\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e\u003cb\u003e\u0026le;\u0026thinsp;0.001\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c10\"\u003e\u003cp\u003e\u003cb\u003e9,8\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c11\"\u003e\u003cp\u003e\u003cb\u003e97,0\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eP(CO-CO\u003csub\u003e2\u003c/sub\u003e-H\u003csub\u003e2\u003c/sub\u003e-H\u003csub\u003e2\u003c/sub\u003eO)\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eLin.\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e-0,03\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e=\u0026thinsp;0.36\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e67,9\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e\u003cb\u003eLin.\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e\u003cb\u003e-0,7\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e\u003cb\u003e\u0026le;\u0026thinsp;0.001\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c10\"\u003e\u003cp\u003e\u003cb\u003e34,7\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c11\"\u003e\u003cp\u003e\u003cb\u003e99,5\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003eCO\u003csub\u003e2\u003c/sub\u003e-H\u003csub\u003e2\u003c/sub\u003eO\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003eLin.\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e-9,5\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e=\u0026thinsp;0.15\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e24,1\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003eLin.\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e-27,2\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e=\u0026thinsp;0.44\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c10\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c11\"\u003e\u003cp\u003e18,1\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003ctr\u003e\u003ctd align=\"left\" colname=\"c1\"\u003e\u003cp\u003e\u003cem\u003eConstant\u003c/em\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c2\"\u003e\u003cp\u003e\u003cem\u003e-\u003c/em\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c3\"\u003e\u003cp\u003e18\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c4\"\u003e\u003cp\u003e\u0026le;\u0026thinsp;0.001\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c5\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c6\"\u003e\u003cp\u003e\u003cb\u003e-\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c7\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c8\"\u003e\u003cp\u003e18,8\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c9\"\u003e\u003cp\u003e\u0026le;\u0026thinsp;0.001\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c10\"\u003e\u003cp\u003e-\u003c/p\u003e\u003c/td\u003e\u003ctd align=\"left\" colname=\"c11\"\u003e\u003cp\u003e\u003cb\u003e-\u003c/b\u003e\u003c/p\u003e\u003c/td\u003e\u003c/tr\u003e\u003c/tbody\u003e\u003c/colgroup\u003e\u003c/table\u003e\u003c/div\u003e\u003c/p\u003e\u003cp\u003eTable\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e. \u003cb\u003eDeformation regression models with (PANEL A) and without (PANEL B) LAG autoregressive term.\u003c/b\u003e Bold face marks the significant factors of each final model. AIC\u0026thinsp;=\u0026thinsp;Akaike information criterion. BIC\u0026thinsp;=\u0026thinsp;Bayesian Information Criterion. R\u003csup\u003e2\u003c/sup\u003ep\u0026thinsp;=\u0026thinsp;Fraction of explained variance. Funct. Form\u0026thinsp;=\u0026thinsp;Functional form. Coeff. = Coefficient. Lin. = Linear. p\u0026thinsp;=\u0026thinsp;observed alfa error. BIF\u0026thinsp;=\u0026thinsp;Bootstrap inclusion frequency.\u003c/p\u003e\u003cp\u003eFigure \u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eA displays the observed deformation time series (blue dots) overlaid with the LAG model prediction (red line), showing excellent agreement. The residuals (green dots) are minimal and randomly distributed. In contrast, the NO-LAG model (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eB) resulted in increased residuals, lower explanatory power (R\u0026sup2; = 0.97), and significantly higher AIC (114.9) and BIC (123.8) values. Covariate-adjusted partial relationships (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e7\u003c/span\u003e) revealed a strong positive association for the first-lag deformation (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e7\u003c/span\u003eA) and seismic activity (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e7\u003c/span\u003eB). Heat flow (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e7\u003c/span\u003eC) showed a negative linear relationship, while CO₂ flux (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e7\u003c/span\u003eD) displayed a nonlinear response. The bootstrap inclusion frequency (BIF) analysis confirmed the LAG model's robustness: the LAG term and seismicity exceeded a 95% inclusion rate, while CO₂ flux and P(CO\u0026ndash;CO₂\u0026ndash;H₂\u0026ndash;H₂O) surpassed the 50% threshold. Conversely, parameters such as temperature and the CO₂/H₂O ratio displayed low inclusion frequencies and limited statistical significance. In the NO-LAG model, seismicity and equilibrium pressure of hydrothermal system emerged as dominant predictors, explaining 41.8% and 34.7% of the variance, respectively. Figure\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e8\u003c/span\u003e illustrates nonlinear relationships and inverted trends among the five most relevant predictors, including negative associations for temperature and pressure, and positive trends for heat flow. Despite strong individual contributions, the NO-LAG model underperformed relative to the LAG model across all evaluation metrics.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec10\" class=\"Section2\"\u003e\u003ch2\u003e3.2 Identification of critical transitions using GCPA\u003c/h2\u003e\u003cp\u003eThe analysis was performed over two overlapping time windows: 2018\u0026ndash;2023 and 2018\u0026ndash;2024. The first interval was selected to identify transitional behaviors prior to the onset of the 2023 critical point (CP2), while the extended window allows for the detection of the additional system-wide reorganization that characterizes CP2. This approach supports a robust validation of the Global Critical Point Analysis (GCPA) framework across different temporal baselines.\u003c/p\u003e\u003cp\u003eWe used the outputs of the MFPA as input for the Global Critical Point Analysis (GCPA), focusing on six variables: ground deformation, seismicity, temperature, heat flow, CO₂ flux, and equilibrium pressure of hydrothermal system P(CO-CO₂-H₂-H₂O). We ran GCPA over two time intervals: 2018\u0026ndash;2023 and 2018\u0026ndash;2024. The analysis revealed two distinct critical transitions:\u003c/p\u003e\u003cp\u003e\u003cul\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003eCP1 \u0026ndash; November 30th, 2020 (\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e9\u003c/span\u003eA\u003cb\u003e)\u003c/b\u003e: This transition was characterized by synchronized peaks in CO₂ flux (\u0026gt;\u0026thinsp;3.5 \u0026times; 10⁴ g/m\u0026sup2;/day), equilibrium pressure of hydrothermal system. Seismic activity showed a sustained increase, and the rate of ground deformation rose sharply, indicating intensified magmatic degassing and increased pressurization within the system.\u003c/p\u003e\u003c/li\u003e\u003c/ul\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003cul\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003eCP2 \u0026ndash; April 1st, 2023 (\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e9\u003c/span\u003eB\u003cb\u003e)\u003c/b\u003e: This event coincided with a new phase of unrest. Seismicity peaked, ground deformation accelerated, and multiple geochemical indicators as CO₂ Flux and equilibrium pressure, reached critical values.\u003c/p\u003e\u003c/li\u003e\u003c/ul\u003e\u003c/p\u003e\u003cp\u003eFinally, we highlight that the integration of individual time series into composite metrics, as shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e9\u003c/span\u003e, enabled the identification of significant changes and evolutions within the volcanic system. Figure\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e9\u003c/span\u003eA highlights the synchronization of variables leading to CP1, while Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e9\u003c/span\u003eB illustrates the multivariable interaction underlying CP2. These panels effectively demonstrate how different parameters interact and highlight critical transitions.\u003c/p\u003e\u003cp\u003eIn this context, the interpretation of the acquired data in light of the performed analysis results highlights that:\u003c/p\u003e\u003cp\u003e\u003col style=\"list-style-type:lower-roman;\"\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eOn the temporal evolution of monitored parameters: Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e9\u003c/span\u003eA shows that, approaching CP1, CO₂ flux, equilibrium temperature, ground deformation, seismicity, land surface temperature, and heat flux all exhibited a coordinated increase. This synchronized behavior suggests a coherent systemic shift leading to the first critical point.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eOn the behavior of variables near CP2: Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e9\u003c/span\u003eB indicates that, prior to CP2, trends across the same set of variables became more differentiated. While CO₂ flux, seismicity, and ground deformation peaked distinctly, other parameters such as heat flux and land surface temperature showed more irregular evolution. This divergence suggests a more complex, possibly fragmented transition.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eOn the contribution of the composite curve: The composite curves in both panels effectively integrate the normalized trajectories of all included parameters. These aggregated trends clearly delineate the timing and intensity of CP1 and CP2, serving as robust indicators of multiparametric transitions.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003cspan\u003e\u003cli\u003e\u003cp\u003eOn the identification of critical transitions: The temporal co-evolution of CO₂ flux, equilibrium temperature, ground deformation, seismicity, land surface temperature, and heat flux confirms the occurrence of two major system reorganizations CP1 and CP2 marking distinct phases in the recent dynamics of the Solfatara-Pisciarelli volcanic system.\u003c/p\u003e\u003c/li\u003e\u003c/span\u003e\u003c/ol\u003e\u003c/p\u003e\u003c/div\u003e"},{"header":"4. Discussion","content":"\u003cdiv id=\"Sec12\" class=\"Section2\"\u003e\u003ch2\u003e4.1 The role of time-lagged regression in ground deformation modeling\u003c/h2\u003e\u003cp\u003eThe proposed methodology allows us to reveal a strong overall association between the variables of interest and the set of factors (i.e., physical parameters) tested in the models. Specifically, regarding ground deformation, seismic activity plays a significant role, accounting for one-third of the overall variance even when a lag term is included in the model. In other words, it explains a substantial portion of the variance in a model while also accounting for unmeasured associative factors. This evidence is further supported by the high stability of seismic activity in the final models, as indicated by its consistently high frequency of statistical significance (BIF above 97%, Table\u0026nbsp;\u003cspan refid=\"Tab2\" class=\"InternalRef\"\u003e2\u003c/span\u003e) in the bootstrapped datasets. In models without the LAG term, other factors gain prominence in the association, attempting to compensate for the 'unmeasured' influences. The high BIF values observed for almost all considered variables in these models reflect this compensatory effect. In particular, for a substantial number of factors (except seismic activity), we observed non-significant levels paired with high BIF values, and vice versa. This suggests that none of these factors, aside from seismic activity, establish a stable association.\u003c/p\u003e\u003cp\u003eThe LAG model effectively accounts for the delays between independent variables in the historical time series. It demonstrates a strong correspondence between observed data (\"Raw data\") and the model fit (\"Fit\"), with residuals (\"Residual\") remaining close to zero (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eA). This indicates that the model successfully captures the temporal delays between independent variables (seismicity, heat flux, and CO₂ flux) and ground deformation.\u003c/p\u003e\u003cp\u003eBy including a LAG term, the model improves predictive accuracy by accounting for the dynamics of the cause-effect relationships between these independent parameters and the dependent variable. In contrast, in the NO-LAG model, the residuals exhibit greater variability and a systematic temporal oscillation, suggesting that the model fails to adequately explain the temporal delays between the independent variables and ground deformation.\u003c/p\u003e\u003cp\u003eIn complex systems like the CF volcanic field, processes such as fluid transport or gas release require time to propagate, leading to delays between the cause (variations in independent variables) and the effect (ground deformation). The inclusion of this time interval in the lag model allows for a more accurate representation of these dynamic processes, resulting in a closer fit between observed data and the model. The near-stationary residuals close to zero in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eA further confirm the model's effectiveness in capturing temporal delays.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec13\" class=\"Section2\"\u003e\u003ch2\u003e4.2 Identification of critical transitions through GCPA\u003c/h2\u003e\u003cp\u003eFurthermore, the application of the GCPA has been fundamental in identifying the timing and characteristics of critical transitions within the volcanic system. The GCPA method enabled us to detect two major critical points (CP1: November 30th 2020, CP2: April 1st 2023), which correspond to significant variations in multiple monitored parameters, including deformation rates, CO₂ flux, and seismicity. These transitions suggest that the SP system underwent two significant evolutionary phases, likely reflecting pressure accumulation and subsequent adjustments within the magmatic-hydrothermal environment.\u003c/p\u003e\u003cp\u003eThe construction of composite metrics by integrating multiple normalized time series allowed us to observe the synchronized evolution of key geophysical and geochemical parameters. This aggregation process revealed collective behaviors and system-wide responses, such as those leading to CP1 and CP2, that would have been difficult to identify through the analysis of individual variables alone. By emphasizing moments of convergence among signals, composite metrics provide a valuable tool for detecting emerging instability in complex volcanic systems.\u003c/p\u003e\u003cp\u003eThe first transition (CP1) marked the onset of an instability phase, characterized by an abrupt increase in gas fluxes and localized ground deformation acceleration. A subsequent stabilization phase followed, leading up to the second transition (CP2). It was associated with an additional shift in seismic and geochemical parameters, indicating renewed activity in the SP system. Identifying such transitions is crucial for volcanic hazard assessment, as it provides insights into the system\u0026rsquo;s evolution over time and highlights periods of potential increased unrest.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec14\" class=\"Section2\"\u003e\u003ch2\u003e4.3 Methodological advantages over classical and non-parametric approaches\u003c/h2\u003e\u003cp\u003eOne of the major strengths of this study is the integration of advanced statistical modeling (MFPA) with critical transition detection (GCPA), providing a novel framework for revealing instability phase in a volcanic system. This approach offers three main advantages:\u003c/p\u003e\u003cp\u003e\u003cul\u003e\u003cli\u003e\u003cp\u003erefined distinction between transient fluctuations and systemic transitions. The methodology may reduce the risk of false positives in hazard assessment by distinguishing between short-term variations and significant shifts in system dynamics.\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003edetection of major critical transitions. The identification of two major critical transitions, corresponding to significant changes in deformation rates, CO₂ flux, and seismicity, highlights the importance of data-driven approaches for monitoring evolving volcanic conditions.\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003escalability to other volcanic systems. The validated methodology can be extended to other volcanic systems, enhancing global monitoring strategies by applying similar multiparametric analyses in high-risk calderas.\u003c/p\u003e\u003c/li\u003e\u003c/ul\u003e\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec15\" class=\"Section2\"\u003e\u003ch2\u003e4.4 Link between critical transitions and overpressure source\u003c/h2\u003e\u003cp\u003eOur results allow to make a connection between the identified critical transition points (CP1 and CP2) and the evolution of the shallow overpressure source presented by (20) achieved using 4D geodetic tomographic inversion techniques. The GCPA results highlight a distinct critical point corresponding to a major shift in system dynamics, indicating a transition from stable conditions to an accelerated phase of pressure accumulation and deformation. The temporal coincidence between the GCPA-derived critical point and the inferred migration of the overpressure source suggests a strong linkage. This is supported by the observed increase in CO₂ flux and seismic activity recorded prior to the modeled overpressure source migration, as well as by the acceleration in ground deformation patterns detected through InSAR data. Specifically, the increase in CO₂ flux, deformation rates, and seismicity, key indicators in the GCPA model, precedes the modeled ascent of pressurized fluids. This supports the hypothesis that the critical point represents the onset of enhanced permeability and fluid mobilization, facilitating the upward migration of the over-pressured magmatic-hydrothermal mixture.\u003c/p\u003e\u003c/div\u003e\u003cdiv id=\"Sec16\" class=\"Section2\"\u003e\u003ch2\u003e4.5 Applicability and perspectives for development\u003c/h2\u003e\u003cp\u003eThe integrated MFPA\u0026ndash;GCPA framework presented here is not restricted to the Solfatara-Pisciarelli area, but can be generalized to other volcanic systems with dense multiparametric monitoring. The core advantage lies in the model's ability to capture nonlinear and time-delayed interactions between diverse geophysical and geochemical signals, without relying on pre-set thresholds or specific event forecasts. For example, volcanoes such as Yellowstone, Long Valley, or Taal \u0026mdash; where persistent deformation, degassing, and shallow seismicity coexist \u0026mdash; could benefit from this approach, particularly in identifying early signs of systemic reorganization or slow-onset unrest. Similarly, densely monitored stratovolcanoes like Etna or Stromboli, where cyclic eruptive behaviors are observed, may provide ideal testbeds for implementing time-lagged multivariate models.\u003c/p\u003e\u003cp\u003eThe case of Campi Flegrei is particularly suited to demonstrate this methodology, due to the exceptional continuity of time series, the known recurrence of hydrothermal/magmatic transitions, and the availability of high-quality deformation, gas and seismic datasets. In this context, the proposed framework enhances our ability to recognize critical transitions that may remain hidden in traditional analyses.\u003c/p\u003e\u003cp\u003eWe anticipate that applying this approach across different volcanic settings will improve comparative analyses, foster cross-caldera learning, and contribute to harmonized early warning protocols at global scale.\u003c/p\u003e\u003cp\u003eThe accuracy of the identified critical point depends on the resolution and consistency of input datasets, and uncertainties in model assumptions regarding fluid migration dynamics may affect interpretations. Further validation through independent geophysical and geochemical observations is necessary to refine these findings for identifying and interpreting critical transitions in volcanic systems.\u003c/p\u003e\u003cp\u003eThis study demonstrates that the combination of MFPA and GCPA provides an effective tool for understanding the dynamics of complex volcanic systems. The identification of critical transitions in the SP area underscores the potential of this methodology in early warning systems and risk assessment frameworks.\u003c/p\u003e\u003cp\u003eFurthermore, the introduction of a time-lagged model significantly improves the predictive power of ground deformation analysis, capturing key interactions between geophysical and geochemical parameters. These insights are particularly relevant for optimizing volcanic monitoring networks and anticipating hazardous events with greater accuracy.\u003c/p\u003e\u003cp\u003eCompared to non-parametric approaches such as decision trees, random forests, or k-nearest neighbors (k-NN), the Multivariable Fractional Polynomial Analysis (MFPA) offers several advantages. First, it allows the identification of interpretable functional relationships between variables, including the assessment of non-linear trends within a parametric framework. This aspect is crucial when the physical meaning of coefficients and lag structures is necessary to relate model output to volcanic processes.\u003c/p\u003e\u003cp\u003eSecond, the MFPA enables the partitioning of explained variance and the decomposition of predictor contributions (via R\u0026sup2;p and BIF), providing quantitative insights into which parameters are most strongly associated with observed dynamics. In contrast, non-parametric methods are often \u0026ldquo;black-box\u0026rdquo; in nature, lacking direct interpretability unless paired with post-hoc explainability tools\u003c/p\u003e\u003cp\u003eThird, parametric models are computationally efficient and robust when applied to moderate-sized time series, as in our case, and they require fewer tuning parameters or training iterations. This makes them particularly suitable for real-time applications and for systems where interpretability and parsimony are as critical as predictive power.\u003c/p\u003e\u003cp\u003eNonetheless, future developments may incorporate ensemble strategies where MFPA is complemented by non-parametric classifiers or anomaly detectors to enhance forecasting capabilities and sensitivity to rare events.\u003c/p\u003e\u003c/div\u003e"},{"header":"5. Conclusions","content":"\u003cp\u003eThe joint application of GCPA and MFPA represents significant methodological advancement, as it enables both the identification of non-linear interactions between variables and the detection of critical transitions in volcanic activity. Modeling temporal delays and systemic tipping points is crucial for understanding the evolution of volcanic unrest and supporting decision-making processes in hazard mitigation. To further develop this approach, future research should focus on: (i) expanding the dataset to include real-time monitoring data and high-resolution satellite imagery for enhanced early warning signals detection; (ii) applying this methodology to other active volcanic systems to validate its robustness and scalability across different geological contexts; (iii) developing hybrid models that integrate AI-driven pattern recognition with our current framework to enhance forecasting capabilities.\u003c/p\u003e\u003cp\u003eThis study utilized a uniquely detailed dataset for the SP area, enabling a critical step in understanding volcanic system interactions through a specifically calibrated one-dimensional model. This approach paves the way for future multidimensional models and broader regional applications, potentially refining volcanic monitoring and hazard assessment by integrating geophysical and geochemical parameters.\u003c/p\u003e\u003cp\u003eIn conclusion, this research represents a step forward in volcanic hazard assessment, demonstrating the feasibility of a data-driven, multiparametric approach to detecting critical transitions in active volcanic systems. The integration of advanced statistical methods with real-world monitoring data can substantially improve early warning systems, ultimately enhancing the safety and resilience of communities living in volcanic regions.\u003c/p\u003e\u003cp\u003eIn this context, in recent years, the identification of early warning signals (EWS) for critical transitions in volcanic systems has become a central goal for improving monitoring capabilities and response preparedness. Several studies have explored this challenge through approaches ranging from tipping point analysis in complex systems to stochastic modeling and machine learning techniques. Early warning indicators based on statistical metrics such as increasing variance, lag-1 autocorrelation (critical slowing down), and kurtosis have shown promising results in specific case studies (Scheffer et al., \u003cspan citationid=\"CR75\" class=\"CitationRef\"\u003e2009\u003c/span\u003e; Dakos et al., \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2012\u003c/span\u003e; Lenton, \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e2013\u003c/span\u003e). However, their application in multiparametric systems is often hindered by high noise sensitivity and limited robustness under complex dynamics (Bernuzzi \u0026amp; Boers \u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Gualandi and Liu, \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Segall, \u003cspan citationid=\"CR78\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Boers, \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). In comparison, the integrated MFPA\u0026ndash;GCPA framework presents three key advantages:\u003c/p\u003e\u003cp\u003e\u003cul\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003e(i)\u003c/b\u003e It explicitly models nonlinear and time-delayed relationships between geophysical and geochemical parameters\u0026mdash;an aspect rarely addressed in conventional EWS models, which often assume Markovian and instantaneous dynamics (K\u0026eacute;fi et al., \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2014\u003c/span\u003e).\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003e(ii)\u003c/b\u003e It avoids arbitrary thresholds: critical points are identified through mathematical optimization over normalized signals, thereby reducing the risk of false positives due to random fluctuations (\u003cb\u003eBoettiger \u0026amp; Hastings, 2013\u003c/b\u003e).\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003e(iii)\u003c/b\u003e It enables the quantitative and temporal integration of multiple signals into a single composite metric, rather than treating them separately, thus facilitating the detection of systemic transitions rather than isolated local events\u0026mdash;an essential feature in highly interdependent systems such as volcanic ones (Newhall \u0026amp; Pallister \u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Behr et al., \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2024\u003c/span\u003e).\u003c/p\u003e\u003c/li\u003e\u003c/ul\u003e\u003c/p\u003e\u003cp\u003eThe MFPA\u0026ndash;GCPA framework thus represents a hybrid approach: more flexible than classical descriptive models, yet more interpretable and statistically robust than many artificial intelligence-based methods (Bountos et al., \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Behr et al., \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). This combination of interpretability and methodological rigor makes it especially well-suited for real-time volcanic surveillance, where both reliability and transparency of signals are critical.\u003c/p\u003e\u003cp\u003eFrom an operational perspective, the application of this framework offers significant potential as a decision-support tool for civil protection authorities and monitoring agencies, such as INGV or global volcanic alert centers. Specifically:\u003c/p\u003e\u003cp\u003e\u003cul\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003eThe inclusion of time lags\u003c/b\u003e enables earlier detection of transitional phases compared to conventional signals, thereby improving the responsiveness of early warning systems (Papale, \u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Petrosino et al., \u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e2023\u003c/span\u003e).\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003eThe identification of systemic transitions\u003c/b\u003e (CP1 and CP2) provides a temporal map of structural changes within the system, which can inform dynamic updates to hazard levels (Petrosino et al., \u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e2023\u003c/span\u003e).\u003c/p\u003e\u003c/li\u003e\u003cli\u003e\u003cp\u003e\u003cb\u003eThe construction of normalized composite metrics\u003c/b\u003e allows for seamless integration into monitoring dashboards, enhancing the interpretability of signals even for non-technical users (GVMID working group, \u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e2024\u003c/span\u003e; Potter et al., \u003cspan citationid=\"CR68\" class=\"CitationRef\"\u003e2014\u003c/span\u003e).\u003c/p\u003e\u003c/li\u003e\u003c/ul\u003e\u003c/p\u003e\u003cp\u003eThe proposed approach could be implemented in near-real time by integrating existing data streams (e.g., InSAR, seismicity, gas emissions) with an automated MFPA\u0026ndash;GCPA analysis module. This would enhance the predictive capabilities of surveillance networks, providing a robust mechanism for detecting critical phases, even in the absence of clear eruptive precursors (Ebmeier et al., \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2018\u003c/span\u003e; Bevilacqua et al., \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2020\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eFinally, the method's adaptability to diverse geological settings makes this framework a strong candidate for the standardization of next-generation multidisciplinary monitoring strategies, particularly in high-risk caldera systems such as Campi Flegrei, Yellowstone, Taupo, and Toba (Newhall \u0026amp; Pallister, \u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; \u003cb\u003eBiggs et al., 2016\u003c/b\u003e).\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgments\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe thank the DTA.AD004.065 project, \"Development and applications of geophysical methods based on the use of environmental remote sensing data,\" and its subproject DTA.AD004.065.001, \"Geophysical analysis and modeling of environmental processes,\" for partially supported this research activity. UAVs thermal data were supported both by the DPC-INGV B2 project WP2 Task 5 (2019-2021) and by the DPC-INGV agreement Annex A. Seismic data are supported by the DPC-INGV agreement Annex A.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors received no funding for this work.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthor contributions:\u003c/strong\u003e\u0026nbsp; \u0026nbsp;\u003c/p\u003e\n\u003cp\u003eVitale A (AV), Barone A (AB), Marotta E (EM), Vitale D F (DFV), Pepe S (SP), Peluso R (RP), Castaldo R (RC), Avino R (RA), Mercogliano F (FM), Pepe A (AP), Accomando F (FA), Avvisati G (GA), Belviso P (PB), Bellucci Sessa E (EBS), Carandente A (AC), Perrini M (MP), Sansivero F (FS), Tizzani P (PT)\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;Conceptualization: PT, AV, DFV\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;Methodology: PT, AV, DFV, EM\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;Investigation: AB, AP, SP, RA, AC, PB, GA, EBS, FS, RC, FM, EM,FA, RP\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;Visualization: PT, AV, AB, MP\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;Supervision: PT, AV, EM\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;Writing original draft: PT, AV, DFV, EM, AB, RC, FS, MP, SP, RP, EBS\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; Writing review \u0026amp; editing:\u0026nbsp;\u003c/p\u003e\u003cp\u003eData and materials availability: \u0026ldquo;The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.\u0026rdquo;\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eAvino, R., Peluso, R., Marotta, E., Belviso, P., Caliro, S., Carandente, A., Chiodini, G., Macedonio, G., Avvisati, G., Marf\u0026egrave;, B., Thermal energy release measurement with thermal camera: The case of La Solfatara Volcano (Italy). Remote Sens. 11, 167 (2019). https://doi.org/10.3390/rs11020167\u003c/li\u003e\n\u003cli\u003eBarone, A., Gola, G., Caliro, S., Chiodini, G., Tizzani, P., Castaldo, R., Long-term thermo-fluid dynamic modeling of Solfatara hydrothermal system, Campi Flegrei caldera. J. Volcanol. Geotherm. Res. 459, 108277 (2025). https://doi.org/10.1016/j.jvolgeores.2025.108277\u003c/li\u003e\n\u003cli\u003eBehr Y, Cerminara M, Rinaldi AP, Azzaro R. Probabilistic, multi-sensor eruption forecasting: a Bayesian framework. \u003cem\u003eGeophys Res Lett\u003c/em\u003e. 2024;51:e2024GL112029. https://doi.org/10.1029/2024GL112029\u003c/li\u003e\n\u003cli\u003eBernuzzi S, Boers N. Critical slowing-down in dynamical systems driven by non-stationary correlated noise. \u003cem\u003ePhys Rev Res\u003c/em\u003e. 2022;4:013230. https://doi.org/10.1103/PhysRevResearch.4.013230\u003c/li\u003e\n\u003cli\u003eBerardino, P., Fornaro, G., Lanari, R., Sansosti, E., A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 40, 2375\u0026ndash;2383 (2002). https://doi.org/10.1109/TGRS.2002.803792\u003c/li\u003e\n\u003cli\u003eBevilacqua A, Neri A, Genco R, Isaia R, Marzocchi W. Failure forecast on Campi Flegrei GPS data 2011\u0026ndash;2020. \u003cem\u003earXiv preprint\u003c/em\u003e. 2020. arXiv:2007.02756. Available from: https://arxiv.org/abs/2007.02756\u003c/li\u003e\n\u003cli\u003eBevilacqua, A., Neri, A., De Martino, P., Isaia, R., Novellino, A., Tramparulo, F., Vitale, S., Radial interpolation of GPS and leveling data of ground deformation in a resurgent caldera: Application to Campi Flegrei (Italy). J. Geod. 94, 13 (2020). https://doi.org/10.1007/s00190-020-01355-x\u003c/li\u003e\n\u003cli\u003eBevilacqua, A., Neri, A., De Martino, P., Giudicepietro, F., Ricciolino, P., Patra, A., Pitman, E. B., Bursik, M., Voight, B., Flandoli, F., Macedonio, G., Accelerating upper crustal deformation and seismicity of Campi Flegrei caldera (Italy) during the 2000\u0026ndash;2023 unrest. Commun. Earth Environ. 5, 742 (2024). https://doi.org/10.1038/s43247-024-01865-y\u003c/li\u003e\n\u003cli\u003eBiggs J, Pritchard ME. Global satellite observations of magmatic and volcanic deformation: implications for volcano monitoring. In: Pyle D, Mather TA, Biggs J, editors. \u003cem\u003eGlobal Volcanic Hazards and Risk\u003c/em\u003e. Cambridge: Cambridge University Press; 2016. p. 1\u0026ndash;34.\u003c/li\u003e\n\u003cli\u003eBoers N. Observation-based early-warning signals for AMOC collapse. \u003cem\u003eNat Clim Chang\u003c/em\u003e. 2021;11:680\u0026ndash;688. https://doi.org/10.1038/s41558-021-01181-1\u003c/li\u003e\n\u003cli\u003eBoettiger C, Ross N, Hastings A. Early warning signals: The charted and uncharted territories. \u003cem\u003eTheor Ecol\u003c/em\u003e. 2013;6:255\u0026ndash;264. https://doi.org/10.1007/s12080-013-0192-6\u003c/li\u003e\n\u003cli\u003eBountos N, Andronopoulos S, Tzovaras D. Self-supervised contrastive learning for volcanic unrest detection. \u003cem\u003earXiv preprint\u003c/em\u003e. 2022. arXiv:2202.04030. Available from: https://arxiv.org/abs/2202.04030\u003c/li\u003e\n\u003cli\u003eBoyd, S., Vandenberghe, L., Convex Optimization (Cambridge Univ. Press, Cambridge, UK, 2004). https://doi.org/10.1017/CBO9780511804441\u003c/li\u003e\n\u003cli\u003eCaliro, S., Avino, R., Capecchiacci, F., Carandente, A., Chiodini, G., Cuoco, E., Minopoli, C., Rufino, F., Santi, A., Rizzo, A. L., Aiuppa, A., Allocca, V., De Vita, P., Di Vito, M. A., Chemical and isotopic characterization of groundwater and thermal waters from the Campi Flegrei caldera (southern Italy). J. Volcanol. Geotherm. Res. 460, 108280 (2025). https://doi.org/10.1016/j.jvolgeores.2025.108280\u003c/li\u003e\n\u003cli\u003eCaliro, S., Chiodini, G., Avino, R., Carandente, A., Cuoco, E., Di Vito, M. A., Minopoli, C., Rufino, F., Santi, A., Lages, J., Mangiacapra, A., Monteleone, B., Pappalardo, L., Taracs\u0026aacute;k, Z., Tramati, C., Vizzini, S., Aiuppa, A., Escalation of caldera unrest indicated by increasing emission of isotopically light sulfur. Nat. Geosci. 18, 123\u0026ndash;132 (2025). https://doi.org/10.1038/s41561-024-01632-w\u003c/li\u003e\n\u003cli\u003eCardellini, C., Chiodini, G., Frondini, F., Avino, R., Bagnato, E., Caliro, S., Lelli, M., Rosiello, A., Monitoring diffuse volcanic degassing during volcanic unrests: The case of Campi Flegrei (Italy). Sci. Rep. 7, 6757 (2017). https://doi.org/10.1038/s41598-017-06169-5\u003c/li\u003e\n\u003cli\u003eCastaldo, R., Tizzani, P., Solaro, G., Inflating source imaging and stress/strain field analysis at Campi Flegrei Caldera: The 2009\u0026ndash;2013 unrest episode. Remote Sens. 13, 2298 (2021). https://doi.org/10.3390/rs13122298\u003c/li\u003e\n\u003cli\u003eChiodini, G., Vandemeulebrouck, J., Caliro, S., D\u0026rsquo;Auria, L., De Martino, P., Mangiacapra, A., Petrillo, Z., Evidence of thermal driven processes triggering the 2005\u0026ndash;2014 unrest at Campi Flegrei caldera. Earth Planet. Sci. Lett. 414, 58\u0026ndash;67 (2015). https://doi.org/10.1016/j.epsl.2015.01.012\u003c/li\u003e\n\u003cli\u003eChiodini, G., Giudicepietro, F., Vandemeulebrouck, J., Aiuppa, A., Caliro, S., De Cesare, W., Tamburello, G., Avino, R., Orazi, M., D\u0026rsquo;Auria, L., Fumarolic tremor and geochemical signals during a volcanic unrest. Geology 45, 1131\u0026ndash;1134 (2017). https://doi.org/10.1130/G39447.1\u003c/li\u003e\n\u003cli\u003eChiodini, G., Todesco, M., Caliro, S., Del Gaudio, C., Macedonio, G., Russo, M., Magma degassing as a trigger of bradyseismic events: The case of Phlegrean Fields (Italy). Geophys. Res. Lett. 30, 1434 (2003). https://doi.org/10.1029/2002GL016790\u003c/li\u003e\n\u003cli\u003eChiodini, G., Avino, R., Caliro, S., Minopoli, C., Temperature and pressure gas geoindicators at the Solfatara fumaroles (Campi Flegrei). Ann. Geophys. 54, 2 (2011). https://doi.org/10.4401/ag-5002\u003c/li\u003e\n\u003cli\u003eChiodini, G., Caliro, S., De Martino, P., Avino, R., Gherardi, F., Early signals of new volcanic unrest at Campi Flegrei caldera? Insights from geochemical data and physical simulations. Geology 40, 943\u0026ndash;946 (2012). https://doi.org/10.1130/G33251.1\u003c/li\u003e\n\u003cli\u003eChiodini, G., Caliro, S., Avino, R., Bagnato, E., Capecchiacci, F., Carandente, A., Allocca, V., Aiuppa, A., The hydrothermal system of the Campi Flegrei Caldera, Italy. in Campi Flegrei: A Restless Caldera in a Densely Populated Area, G. Orsi, M. D\u0026rsquo;Antonio, L. Civetta, Eds. (Springer, Berlin, 2022), pp. 239\u0026ndash;255. https://doi.org/10.1007/978-3-642-37060-1_9\u003c/li\u003e\n\u003cli\u003eCioni, R., Corazza, E., Medium temperature fumarolic gas sampling. Bull. Volcanol. 44, 23\u0026ndash;29 (1981).\u003c/li\u003e\n\u003cli\u003eCirillo, F., Avvisati, G., Belviso, P., Marotta, E., Peluso, R., Pescione, P. R., STARTED (StaTistical Analysis clusteRed ThErmal Data) (version 1.0.1). Zenodo (2022). https://doi.org/10.5281/zenodo.5886707\u003c/li\u003e\n\u003cli\u003eDakos V, Carpenter SR, Brock WA, Ellison AM, Guttal V, Ives AR, et al. Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. \u003cem\u003ePLoS One\u003c/em\u003e. 2012;7(7):e41010. https://doi.org/10.1371/journal.pone.0041010\u003c/li\u003e\n\u003cli\u003eD\u0026rsquo;Auria, L., Pepe, S., Castaldo, R., Giudicepietro, F., Macedonio, G., Ricciolino, P., Tizzani, P., Casu, F., Lanari, R., Manzo, M., Martini, M., Sansosti, E., Zinno, I., Magma injection beneath the urban area of Naples: A new mechanism for the 2012\u0026ndash;2013 volcanic unrest at Campi Flegrei caldera. Sci. Rep. 5, 13100 (2015). https://doi.org/10.1038/srep13100\u003c/li\u003e\n\u003cli\u003eDel Gaudio, C., Aquino, I., Ricciardi, G.P., Ricco, C., Scandone, R., Unrest episodes at Campi Flegrei: A reconstruction of vertical ground movements during 1905\u0026ndash;2009. J. Volcanol. Geotherm. Res. 195, 48\u0026ndash;56 (2010). https://doi.org/10.1016/j.jvolgeores.2010.05.014\u003c/li\u003e\n\u003cli\u003eDi Filippo, A., Peluso, R., SEREWrap \u0026ndash; SERENADE Wrapper (version 5.3.2). Zenodo (2019). https://doi.org/10.5281/zenodo.14844623\u003c/li\u003e\n\u003cli\u003eDirac, P. A. M., The Principles of Quantum Mechanics, 4th ed., Oxford, UK: Oxford University Press, 1958.\u003c/li\u003e\n\u003cli\u003eEbmeier SK, Elliott JR, Naismith A, Biggs J. Synthesis of global satellite observations of magmatic and volcanic deformation: implications for volcano monitoring. \u003cem\u003eGeol Soc Lond Spec Publ\u003c/em\u003e. 2018;446:305\u0026ndash;322. https://doi.org/10.1144/SP446.5\u003c/li\u003e\n\u003cli\u003eEster, M., Kriegel, H. P., Sander, J., Xu, X., A density-based algorithm for discovering clusters in large spatial databases with noise. in Proc. 2nd Int. Conf. Knowl. Discovery Data Mining (KDD-96), (1996), pp. 226\u0026ndash;231.\u003c/li\u003e\n\u003cli\u003eFinkel, S. S., Linear panel analysis. in Handbook of Longitudinal Research: Design, Measurement, and Analysis, S. Menard, Ed. (Elsevier/Academic Press, Burlington, MA, 2007), pp. 486\u0026ndash;487.\u003c/li\u003e\n\u003cli\u003eFletcher, R., Practical Methods of Optimization (Wiley, Chichester, ed. 2, 2013).\u003c/li\u003e\n\u003cli\u003eGabriel, A. K., Goldstein, R. M., Zebker, H. A., Mapping small elevation changes over large areas: Differential interferometry. J. Geophys. Res. 94, 9183\u0026ndash;9191 (1989). https://doi.org/10.1029/JB094iB07p09183\u003c/li\u003e\n\u003cli\u003eGiacomuzzi, G., Chiarabba, C., Bianco, F., De Gori, P., Piana Agostinetti, N., Tracking transient changes in the plumbing system at Campi Flegrei Caldera. Earth Planet. Sci. Lett. 637, 118744 (2024). https://doi.org/10.1016/j.epsl.2024.118744\u003c/li\u003e\n\u003cli\u003eGiggenbach, W. F., A simple method for the collection and analysis of volcanic gas samples. Bull. Volcanol. 39, 132\u0026ndash;145 (1975).\u003c/li\u003e\n\u003cli\u003eGiggenbach, W. F., Gouguel, R. L., Collection and analysis of geothermal and volcanic water and gas discharges. Chemistry Division Report, DSIR, New Zealand, No. CD 2401 (1989).\u003c/li\u003e\n\u003cli\u003eGiudicepietro, F., Avino, R., Bellucci Sessa, E., Bevilacqua, A., Bonano, M., Caliro, S., Chiodini, G., Burst-like swarms in the Campi Flegrei caldera accelerating unrest from 2021 to 2024. Nat. Commun. 16, 1548 (2025). https://doi.org/10.1038/s41467-025-39587-1\u003c/li\u003e\n\u003cli\u003eGoldberg, D. E., Genetic Algorithms in Search, Optimization, and Machine Learning, Reading, MA, USA: Addison-Wesley, 1989.\u003c/li\u003e\n\u003cli\u003eGualandi A, Liu Z. Variational Bayesian Independent Component Analysis for InSAR displacement time-series: application to the 2009 Sierra Negra unrest. \u003cem\u003eJ Geophys Res Solid Earth\u003c/em\u003e. 2021;126:e2020JB020845. https://doi.org/10.1029/2020JB020845\u003c/li\u003e\n\u003cli\u003eGVMID Working Group. The Global Volcano Monitoring Infrastructure Database: toward integrated dashboards. \u003cem\u003eFront Earth Sci\u003c/em\u003e. 2024;12:1284889. https://doi.org/10.3389/feart.2024.1284889\u003c/li\u003e\n\u003cli\u003eIsaia, R., Vitale, S., Di Giuseppe, M. G., Iannuzzi, E., Tramparulo, F. D., Troiano, A., Stratigraphy, structure, and volcanotectonic evolution of Solfatara maar-diatreme (Campi Flegrei, Italy). Geol. Soc. Am. Bull. 127, 1485\u0026ndash;1504 (2015). https://doi.org/10.1130/B31183.1\u003c/li\u003e\n\u003cli\u003eJohnson, C. E., Bittenbinder, A., Bogaert, B., Dietz, L., Kohler, W., Earthworm: A flexible approach to seismic network processing. IRIS Newsletter 14, 1\u0026ndash;4 (1995).\u003c/li\u003e\n\u003cli\u003eJolivet, R., Grandin, R., Lasserre, C., Doin, M. P., Peltzer, G., Systematic InSAR tropospheric phase delay corrections from global meteorological reanalysis data. Geophys. Res. Lett. 38, L17311 (2011). https://doi.org/10.1029/2011GL048757\u003c/li\u003e\n\u003cli\u003eKarmarkar, N., A new polynomial-time algorithm for linear programming. Combinatorica 4, 373\u0026ndash;395 (1984). https://doi.org/10.1007/BF02579150\u003c/li\u003e\n\u003cli\u003eK\u0026eacute;fi S, Guttal V, Brock WA, Carpenter SR, Ellison AM, Livina V, et al. Early warning signals of ecological transitions: Methods for spatial patterns. \u003cem\u003ePLoS One\u003c/em\u003e. 2014;9(3):e92097. https://doi.org/10.1371/journal.pone.0092097\u003c/li\u003e\n\u003cli\u003eKutner, M. H., Nachtsheim, C. J., Neter, J., Li, W., Applied Linear Statistical Models, 5th ed. (McGraw-Hill/Irwin, New York, NY, 2005), p. 409.\u003c/li\u003e\n\u003cli\u003eLee, W. H. K., Lahr, J. C., HYPO71: A computer program for determining hypocenter, magnitude, and first motion pattern of local earthquakes. U.S. Geol. Surv. Open File Rep. 75\u0026ndash;311 (1975).\u003c/li\u003e\n\u003cli\u003eLenton TM. Environmental tipping points. \u003cem\u003eAnnu Rev Environ Resour\u003c/em\u003e. 2013;38:1\u0026ndash;29. https://doi.org/10.1146/annurev-environ-102511-084654\u003c/li\u003e\n\u003cli\u003eManzo, M., Ricciardi, G. P., Casu, F., Ventura, G., Zeni, G., Borgstrom, S., Berardino, P., Del Gaudio, C., Lanari, R., Surface deformation analysis in the Ischia island (Italy) based on spaceborne radar interferometry. J. Volcanol. Geotherm. Res. 151, 399\u0026ndash;416 (2006). https://doi.org/10.1016/j.jvolgeores.2005.09.010\u003c/li\u003e\n\u003cli\u003eMarotta, E., Peluso, R., Avino, R., Avvisati, G., Bellucci Sessa, E., Belviso, P., Caputo, T., Carandente, A., Cirillo, F., Pescione, R. A., Clusterisation and temporal trends of heat flux by UAS thermal camera. Remote Sens. 16, 1102 (2024). https://doi.org/10.3390/rs16061102\u003c/li\u003e\n\u003cli\u003eMarotta E., R. Peluso, R. Avino, P. Belviso, S. Caliro, A. Carandente, G. Chiodini, G. Macedonio, G. Avvisati, B. Marf\u0026egrave;, Thermal energy release measurement with thermal camera: The case of La Solfatara Volcano (Italy). Remote Sens. 11, 167 (2019). https://doi.org/10.3390/rs11020167\u003c/li\u003e\n\u003cli\u003eMassonnet, D., Rossi, M., Carmona, C., Ardagna, F., Peltzer, G., Feigl, K., Rabaute, T., The displacement field of the Landers earthquake mapped by radar interferometry. Nature 364, 138\u0026ndash;142 (1993). https://doi.org/10.1038/364138a0\u003c/li\u003e\n\u003cli\u003eMercogliano, F., Barone, A., D\u0026rsquo;Auria, L., Castaldo, R., Silvestri, M., Bellucci Sessa, E., Caputo, T., Stroppiana, D., Caliro, S., Minopoli, C., Avino, R., Tizzani, P., Thermal patterns at the Campi Flegrei caldera inferred from satellite data and independent component analysis. Remote Sens. 16, 4615 (2024). https://doi.org/10.3390/rs16234615\u003c/li\u003e\n\u003cli\u003eNewhall CG, Pallister JS. Introducing the Volcanic Unrest Index (VUI): a tool to quantify and communicate the intensity of volcanic unrest. \u003cem\u003eBull Volcanol\u003c/em\u003e. 2015;77:77. https://doi.org/10.1007/s00445-015-0952-y\u003c/li\u003e\n\u003cli\u003eOppenheim, A. V., Schafer, R. W., Discrete-Time Signal Processing, Upper Saddle River, NJ, USA: Prentice-Hall, 1999.\u003c/li\u003e\n\u003cli\u003eOrsi, G., D\u0026rsquo;Antonio, M., Civetta, L., Volcanic and deformation history of the Campi Flegrei volcanic field, Italy. in Campi Flegrei, Eds. (Springer, Berlin, 2022), pp. 1\u0026ndash;28. https://doi.org/10.1007/978-3-642-37060-1_1\u003c/li\u003e\n\u003cli\u003ePapale P. Rational volcanic hazard forecasts and the use of volcanic alert levels. \u003cem\u003eJ Appl Volcanol\u003c/em\u003e. 2017;6(1):4. https://doi.org/10.1186/s13617-017-0064-7\u003c/li\u003e\n\u003cli\u003ePeluso, R., Il database sismico SERENADE: un sistema REST per la gestione delle localizzazioni sismiche. Miscellanea INGV 57, 7 (2020). https://doi.org/10.13127/misc/57/7\u003c/li\u003e\n\u003cli\u003ePeluso, R., SEDAN \u0026ndash; Seismic Event Dispatcher And Notifier (version 0.11.1). Zenodo (2019). https://doi.org/10.5281/zenodo.14845439\u003c/li\u003e\n\u003cli\u003ePeluso, R., SERENADE \u0026ndash; SEismic Restful ENAbled DatabasE (version 0.6.1). Zenodo (2018a). https://doi.org/10.5281/zenodo.14732232\u003c/li\u003e\n\u003cli\u003ePeluso, R., WESSEL \u0026ndash; WEb interface for Seismic Signal Elaboration and Location (version 1.12.1). Zenodo (2018b). https://doi.org/10.5281/zenodo.14727127\u003c/li\u003e\n\u003cli\u003ePepe, A., Solaro, G., Cal\u0026ograve;, F., Dema, C., A minimum acceleration approach for the retrieval of multiplatform InSAR deformation time series. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 9, 3883\u0026ndash;3898 (2016). https://doi.org/10.1109/JSTARS.2016.2577878\u003c/li\u003e\n\u003cli\u003ePepe, S., De Siena, L., Barone, A., Castaldo, R., D\u0026rsquo;Auria, L., Manzo, M., Casu, F., Fedi, M., Lanari, R., Bianco, F., Tizzani, P., Volcanic structures investigation through SAR and seismic interferometric methods: The 2011\u0026ndash;2013 Campi Flegrei unrest episode. Remote Sens. Environ. 234, 111440 (2019). https://doi.org/10.1016/j.rse.2019.111440\u003c/li\u003e\n\u003cli\u003ePetrillo, Z., Chiodini, G., Mangiacapra, A., Caliro, S., Capuano, P., Russo, G., Cardellini, C., Avino, R., Defining a 3D physical model for the hydrothermal circulation at Campi Flegrei caldera (Italy). J. Volcanol. Geotherm. Res. 264, 172\u0026ndash;183 (2013). https://doi.org/10.1016/j.jvolgeores.2013.08.008\u003c/li\u003e\n\u003cli\u003ePetrosino S, De Lauro E, Aquino I, Del Pezzo E. Evolution in unrest processes at Campi Flegrei inferred from multiparametric lags. \u003cem\u003eJ Volcanol Geotherm Res\u003c/em\u003e. 2023;431:107784. https://doi.org/10.1016/j.jvolgeores.2023.107784\u003c/li\u003e\n\u003cli\u003ePotter SH, Becker JS, Johnston DM, Rossiter KP. Communicating the status of volcanic activity: revising New Zealand\u0026rsquo;s VAL system. \u003cem\u003eJ Appl Volcanol\u003c/em\u003e. 2014;3:13. https://doi.org/10.1186/s13617-014-0013-6\u003c/li\u003e\n\u003cli\u003eRabinowitz, P. H., Critical point theory and its applications to differential equations: A survey. in Proc. Symp. Pure Math. 45, 357\u0026ndash;386 (1986).\u003c/li\u003e\n\u003cli\u003eR Core Team: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Avaiable on: https://www.R-project.org/ (2014)\u003c/li\u003e\n\u003cli\u003eRicciolino, P., Lo Bascio, D., Esposito, R., GOSSIP \u0026ndash; Database sismologico Pubblico INGV-Osservatorio Vesuviano. Istituto Nazionale di Geofisica e Vulcanologia (INGV) (2024). https://doi.org/10.13127/gossip\u003c/li\u003e\n\u003cli\u003eRosen, P. A., Hensley, S., Joughin, I. R., Li, F. K., Madsen, S. N., Rodriguez, E., Goldstein, R., Synthetic aperture radar interferometry. Proc. IEEE 88, 333\u0026ndash;382 (2000). https://doi.org/10.1109/5.838084\u003c/li\u003e\n\u003cli\u003eRoyston, P., Sauerbrei, W., Multivariable Model-Building: A Pragmatic Approach to Regression Analysis Based on Fractional Polynomials for Modelling Continuous Variables (Wiley, Chichester, 2008).\u003c/li\u003e\n\u003cli\u003eScotto di Uccio, F., Lomax, A., Natale, J., Muzellec, T., Festa, G., Nazeri, S., Convertito, V., Bobbio, A., Strumia, C., Zollo, A., Delineation and fine-scale structure of fault zones activated during the 2014\u0026ndash;2024 unrest at the Campi Flegrei caldera (Southern Italy) from high-precision earthquake locations. Geophys. Res. Lett. 51, e2023GL107680 (2024). https://doi.org/10.1029/2023GL107680\u003c/li\u003e\n\u003cli\u003eScheffer, M., Bascompte, J., Brock, W. A., Carpenter, S. R., Dakos, V., van Nes, E. H., and Sugihara, G. (2009). Early-warning signals for critical transitions. Nature, 461(7260), 53-59. https://doi.org/10.1038/nature08227\u003c/li\u003e\n\u003cli\u003eShorrocks A. F., Decomposition procedures for distributional analysis: A unified framework based on the Shapley value. J. Econ. Inequal. 11, 99\u0026ndash;126 (2013). https://doi.org/10.1007/s10888-011-9184-0\u003c/li\u003e\n\u003cli\u003eShumway, R. H., Stoffer, D. S., Time Series Analysis and Its Applications: With R Examples, 4th ed. (Springer, New York, 2017). https://doi.org/10.1007/978-3-319-52452-8\u003c/li\u003e\n\u003cli\u003eSegall P. Noise sources and detection limits in volcano geodesy. \u003cem\u003eJ Volcanol Geotherm Res\u003c/em\u003e. 2019;381:106624. https://doi.org/10.1016/j.jvolgeores.2019.106624\u003c/li\u003e\n\u003cli\u003eSiniscalchi, A., Tripaldi, S., Romano, G., Chiodini, G., Improta, L., Petrillo, Z., D\u0026rsquo;Auria, L., Caliro, S., Avino, R., Reservoir structure and hydraulic properties of the Campi Flegrei geothermal system inferred by audiomagnetotelluric. J. Geophys. Res. Solid Earth 124, 5336\u0026ndash;5356 (2019). https://doi.org/10.1029/2018JB016514\u003c/li\u003e\n\u003cli\u003eSmith, V. C., Isaia, R., Pearce, N. J. G., Tephrostratigraphy and glass compositions of post-15 kyr Campi Flegrei eruptions: Implications for eruption history and chronostratigraphic markers. Quat. Sci. Rev. 30, 3638\u0026ndash;3660 (2011). https://doi.org/10.1016/j.quascirev.2011.07.012\u003c/li\u003e\n\u003cli\u003eSteinmann, L., Spiess, V., Sacchi, M., The Campi Flegrei caldera (Italy): Formation and evolution in interplay with sea-level variations since the Campanian Ignimbrite eruption at 39ka. J. Volcanol. Geotherm. Res. 327, 361\u0026ndash;374 (2016). https://doi.org/10.1016/j.jvolgeores.2016.09.001\u003c/li\u003e\n\u003cli\u003eSuleiman, A. A., Analysis of multicollinearity in multiple regressions. Int. J. Adv. Technol. Eng. Sci. 3, 571\u0026ndash;578 (2015).\u003c/li\u003e\n\u003cli\u003eSze, S. M., Physics of Semiconductor Devices, 2nd ed., New York, NY, USA: Wiley, 1981.\u003c/li\u003e\n\u003cli\u003eTheu\u0026szlig;l, S., Schwendinger, F., Hornik, K., ROI: An extensible R optimization infrastructure. J. Stat. Softw. 94, 1\u0026ndash;64 (2020). https://doi.org/10.18637/jss.v094.i15\u003c/li\u003e\n\u003cli\u003eTizzani, P., Fern\u0026aacute;ndez, J., Vitale, A., Escayo, J., Barone, A., Castaldo, R., Pepe, S., De Novellis, V., Solaro, G., Pepe, A., Tramelli, A., Hu, Z., Samsonov, S. V., Vigo, I., Tiampo, K. F., Camacho, A. G., 4D imaging of the volcano feeding system beneath the urban area of the Campi Flegrei caldera. Remote Sens. Environ. 315, 114480 (2024). https://doi.org/10.1016/j.rse.2024.114480\u003c/li\u003e\n\u003cli\u003eTramelli, A., Convertito, V., Godano, C., B value enlightens different rheological behaviour in Campi Flegrei caldera. Commun. Earth Environ. 5, 1447 (2024). https://doi.org/10.1038/s43247-024-01447-y\u003c/li\u003e\n\u003cli\u003eTroiano, A., Di Giuseppe, M. G., Isaia, R., 3D structure of the Campi Flegrei caldera central sector reconstructed through short-period magnetotelluric imaging. Sci. Rep. 12, 20802 (2022). https://doi.org/10.1038/s41598-022-24998-6\u003c/li\u003e\n\u003cli\u003eTroiano, A., Isaia, R., Di Giuseppe, M. G., Tramparulo, F., Mormone, F., Vitale, V., D\u0026rsquo;Auria, L., Ricciardi, G. P., Deep electrical resistivity tomography for a 3D picture of the most active sector of Campi Flegrei caldera. Sci. Rep. 9, 15124 (2019). https://doi.org/10.1038/s41598-019-51568-0\u003c/li\u003e\n\u003cli\u003eVanorio, A., Geremia, G., De Landro, M., Guo, Z., The recurrence of geophysical manifestations at the Campi Flegrei caldera. Sci. Adv. 11, eadt2067 (2025). https://doi.org/10.1126/sciadv.adt2067\u003c/li\u003e\n\u003cli\u003eWickham, H., ggplot2: Elegant Graphics for Data Analysis (Springer, New York, 2009). https://doi.org/10.1007/978-0-387-98141-3\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Volcanology, Ground deformation, Seismicity, Gas emissions, Critical transitions, Early warning","lastPublishedDoi":"10.21203/rs.3.rs-7290139/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7290139/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eUnderstanding how volcanic systems evolve over time is a major challenge due to their complex behavior and constantly changing conditions. This study explores a new way to detect early warning signals of volcanic unrest by analyzing how different types of data, such as ground deformation, gas emissions, temperature, and earthquakes, interact with each other. Focusing on the Solfatara-Pisciarelli volcano system, that is a more active area in Campi Flegrei Caldera (Southern Italy), we used two advanced methods to identify critical transitions in the system: one to model the nonlinear relationships between variables, and the other to detect key moments when the system's behavior shifts. By including time delays between signals (LAG) we found that our model becomes much more accurate in identifying these changes. In contrast, models that ignore time lags show higher uncertainty. These results highlight the importance of considering how different warning signs may not appear all at once, and they offer a promising tool for improving the monitoring of active volcanoes and supporting early risk management.\u003c/p\u003e","manuscriptTitle":"Critical Transitions at Campi Flegrei resurgent caldera: A Novel Approach to Detect Early Warning Signals","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-08-19 04:41:40","doi":"10.21203/rs.3.rs-7290139/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"52622bb4-e372-4572-8bcf-fb1a843cebdd","owner":[],"postedDate":"August 19th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[{"id":53148222,"name":"Earth and environmental sciences/Natural hazards"},{"id":53148223,"name":"Earth and environmental sciences/Solid earth sciences"}],"tags":[],"updatedAt":"2026-02-10T15:25:17+00:00","versionOfRecord":[],"versionCreatedAt":"2025-08-19 04:41:40","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-7290139","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7290139","identity":"rs-7290139","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}