1. Introduction: “Campi Flegrei” Brief History and Current Scenario
The Campi Flegrei, Pozzuoli, Italy, is a caldera volcanic field that has been active for at least 80 kyr. Campi Flegrei is the largest active caldera in Europe (∼12 km diameter) with a population of more than 360,000 people. From the suburbs of Naples, it extends westward to the Tyrrhenian Sea and is partially submerged below the Gulf of Pozzuoli (
Figure 1). The oldest rock outcroppings, in fact, date back to ca. 80 kyr.; however, much older explosive activity is reported in the literature [
1]. Campi Flegrei eruptive history has been characterized by at least two large-magnitude events that generated the 12 km wide caldera: the eruption of the Campanian Ignimbrite, which occurred about 40,000 years ago, and the Neapolitan Yellow Tuff eruption [
2], aged 15 ka [
3], which occurred 39 ka ago [
4,
5]. Recent studies suggest the existence of a third high energy event that occurred around 29,000 kyr ago, that is, the eruption of the Masseria del Monte tuff [
6].
The volcanic activity preceding the eruption of the Campanian Ignimbrite is evidenced by pyroclastic deposits and lava emerging only along the slopes bordering the Campi Flegrei and in excavations and drillings in the urban area of Naples [
7,
8]. They represent the products of at least 20 explosive and effusive eruptions generated by eruptive centers that are rarely visible and, in most cases, have been buried under more recent products or hidden due to the intense urbanization of the area.
The eruptive activity of the last 15,000 years, therefore following the eruption of the Neapolitan Yellow Tuff, was mainly explosive, predominantly phreatomagmatic, and produced fall deposits with interspersed pyroclastic current deposits with the construction of cinder cones and tuff cones. Effusive eruptions, generating lava domes in the central–eastern area of the caldera, and low-energy Strombolian eruptions also occurred, though less frequently.
This activity was interrupted by two long periods of quiescence that allowed it to be divided into three different eras [
9,
10]: first epoch of activity (15,000–10,600 years); second epoch of activity (9600–9200 years); and third epoch of activity (5500–3800 years).
In total, approximately 70 eruptions were identified during the three epochs, mainly explosive, with eruptive centers largely located in the central–eastern sector of the caldera. The eruptions of Pomici (12,300 years ago) and Agnano-Monte Spina (4550 years ago) are among the highest energy eruptions.
The last eruption, which occurred in 1538, produced the Monte Nuovo tuff cone in the western sector of the caldera, just west of the town of Pozzuoli. In the years preceding this eruption, the caldera underwent a significant ground uplift displaying a bell-shaped pattern. After the eruption, the same area was subject to constant subsidence until the mid-20th century [
11]. Starting from the second half of the 20th century, a period of uplift episodes began, which led to two bradyseismic crises in the years 1950–52, 1969–72, and 1982–84, the latter characterized by intense seismicity. The term “bradyseismic” indicates that phenomenon by which the soil is subject to slow vertical movements of subsidence or uplift.
As highlighted in
Section 3, currently, the ground has been constantly and gradually rising since 2005.
Despite considerable progress in monitoring and modeling volcanic activity at Campi Flegrei, significant challenges persist in comprehensively capturing the spatial and temporal patterns of bradyseismic events. Numerous studies have contributed to advancing knowledge in this field, employing a range of geophysical, geochemical, and remote sensing methodologies. This article aims to build on this foundation by providing a comprehensive review of relevant research, with a focus on the methodologies that have been applied in the Campi Flegrei region and their broader implications for understanding bradyseismic dynamics.
The scientific exploration of bradyseismic activity in Campi Flegrei dates back to the mid-20th century, with early studies primarily relying on ground-based measurements such as leveling surveys and borehole tiltmeters. Key findings from these foundational investigations established the cyclic nature of ground deformation and identified significant episodes of uplift, notably in 1969–1972 and 1982–1984. These events were linked to magmatic intrusions and hydrothermal fluid migrations beneath the caldera. For instance, Barberi et al. [
12] and Dvorak and Mastrolorenzo [
13] provided critical insights into the interplay between magma chamber pressurization and fault reactivation during uplift phases.
Many studies have focused on the study of geophysical spatial dynamics, while few studies have employed proper geographic methods of spatial analysis to complement geophysical analyses.
The advent of satellite-based remote sensing technologies revolutionized the study of bradyseismic phenomena. Synthetic Aperture Radar Interferometry (InSAR) has emerged as one of the most powerful tools for detecting and quantifying ground deformation over large spatial extents. Some studies [
14,
15] employed InSAR data to map ground deformation patterns in unprecedented detail, revealing spatial asymmetries in uplift and subsidence across the caldera. Global Navigation Satellite System (GNSS) networks have also played a pivotal role in tracking the temporal evolution of ground movements. GNSS data have provided high-resolution time series of deformation, enabling researchers to correlate uplift episodes with seismic swarms and changes in fumarolic activity. More recently, integrated approaches combining InSAR and GNSS data have improved the accuracy and reliability of deformation models, as demonstrated in [
16].
Geophysical modeling has been instrumental in interpreting the observed deformation patterns and constraining the underlying sources of bradyseismic activity. Analytical and numerical models have been developed to simulate the pressurization of magmatic and hydrothermal reservoirs and their interaction with caldera structures. One widely cited model [
14] employed a finite-element approach to simulate stress distribution and fault slip during uplift episodes, highlighting the role of structural heterogeneities in modulating deformation dynamics.
Thermo–hydro–mechanical–chemical (THMC) models have further advanced our understanding of the coupled processes driving bradyseismic activity. These models incorporate the effects of heat transfer, fluid flow, and chemical reactions on deformation patterns, providing a more comprehensive representation of subsurface processes. For example, Tramelli et al. [
17] demonstrated how the interplay between magmatic degassing and hydrothermal circulation could explain the observed uplift and seismicity at Campi Flegrei.
Similar geographical approaches could be recognized in recent years in studies regarding the application of GIS, machine learning, and big data analytics to bradyseismic research. Machine learning algorithms have been used to analyze large volumes of geophysical data, identify patterns, and forecast deformation trends. Works such as those by Bellucci Sessa et al. [
18], Felpeto et al. [
19], Pignone et al. [
20] have highlighted the potential of these geography-based approaches to improve hazard assessment by integrating multi-source data, including seismic, geodetic, and geochemical observations.
The multiparameter approach to earthquake forecasting based on complementary geographical approaches integrates geospatial and geophysical data across multiple disciplines to improve earthquake prediction. This strategy leverages the spatial variability of tectonic, environmental, and societal factors to identify patterns and anomalies associated with seismic events [
21]. Similar research using this methodology can be recognized in Pulinets et al. [
22].
The approach makes use of geospatial analysis, through Geographic Information Systems (GISs), to map and analyze spatial correlations between seismic activity and potential precursors.
The data entry level is obtained through the integration of spatial data, including fault maps, seismic records, crustal deformation patterns, and land use information, acquired by a distributed monitoring network, particularly through the deployment of dense networks of seismic, GPS, and strain monitoring stations tailored to regional conditions. Regional variability in precursor phenomena, such as gas emissions or electromagnetic anomalies, is in fact used to refine forecasts.
In this sense, the geographical complementary approach based on the spatial integration of parameters could facilitate the analysis of how different precursors manifest across geographical settings. Furthermore, seismic clustering could be identified through mapping spatial patterns of foreshocks and aftershocks to identify high-risk zones. Combining geology, geophysics, and geography allows for more holistic models of earthquake risk, and in terms of scalability, the geospatial tools enable multi-scale analysis, from local fault zones to regional tectonic systems.
The methodologies employed in Campi Flegrei have demonstrated broad applicability to other volcanic regions experiencing episodic deformation. While, in fact, Campi Flegrei remains one of the most studied calderas for bradyseismic activity, analogous research has been conducted in other volcanic regions, such as Yellowstone (USA), Long Valley (USA), and Rabaul (Papua New Guinea). These studies have provided valuable comparative insights into the mechanisms of caldera deformation and their relationship to magmatic and hydrothermal processes. For instance, the deformation episodes at Rabaul caldera, analyzed by Saunders et al. [
23], exhibit striking similarities to those at Campi Flegrei, underscoring the importance of the regional context in interpreting bradyseismic phenomena.
2. Geographical Perspective as an Added Value
Considering the high probability, highlighted by seismologists, of the repetition of other events in the near future, it is necessary to prepare and make effective a Preventive Integrated Synergistic Strategy (SSIP) that is not dictated by emergency triggers, but rather represents a policy characterized by different monitoring stages, in a range of variation that oscillates between normality and severity, capable of dealing with a phenomenon, such as that at the Campi Flegrei, of a long duration, with critical phases of worsening of intensity, and multifactorial nature. Therefore, a similar strategy can only be developed through the interoperability of scientific areas that are very different but certainly complementary to each other.
The UNISDR (United Nations Office for Disaster Reduction) expresses itself in relation to reducing the risk of dangerous events by highlighting the need for systematic efforts, therefore, codified and organized preventive measures, aimed at predicting, managing, and containing adverse risk factors, mitigating the vulnerability of things and people [
24]. According to this approach, both a horizontal-type coordination between different disciplines and a vertical-type one involving many governments, institutions, NGOs, and communities are necessary, the effectiveness of which is reflected in the systematic application of protection actions as opposed to extemporaneous mechanisms [
25].
The role of geography in the dialogue with other disciplines and its supporting role in multidisciplinary cooperative activities has often been reiterated. In a work from a few years ago, several geographers under the collective acronym of GECO, of the AGeI group “Geography, development cooperation and local development”, contributed to highlighting, for example, the relationship of the discipline with development cooperation [
14]. More recently, in November 2019, in Pisa, in the meeting called “Interpreting the fourth industrial revolution: geography in dialogue with other disciplines”, organized as part of the Excellence Project of the Department of Civilization and Forms of Knowledge with the collaboration of the Association of Italian Geographers and the Society of Geographical Studies, we wanted to reflect on the usefulness of exchanging ideas, models, and methodologies between different disciplines.
In relation to the earthquake situation in the Campi Flegrei area, the common denominator when reflecting on the synergy of disciplines and, in particular, from the geographical perspective refers to its relationship with risk. In this way, it is necessary to keep in mind that knowledge of the risk levels associated with certain geographical areas and that can be deduced from thematic maps created before a possible threat constitutes a fundamental prerequisite for the effectiveness of the efforts made in disaster mitigation.
Concerning the preventive knowledge of the geographical context for the purposes of reducing natural threats, two perspectives can be recognized, i.e., one endogenous and one exogenous. The first relates to the case in which an area is home to endemic risky phenomena. Typical cases are represented by areas exposed to volcanic and seismic risk, in which there is a systematic situation with the possible triggering of some known and predictable risk factors. In the case of such events, evacuation plans are drawn up, and progressive bands of increasing risk are defined. The second concerns external dangers that can be channeled into certain areas, e.g., in the circumstance that the hazard is not specific to a specific territory, and therefore its peculiar virulence of impact must be assessed in relation to the territorial characteristics that can hinder or favor its spread and effects. Precisely, with reference to diffusion phenomena, it is possible to use a model, such as geographical scattering [
26], that allows for obtaining a spatial reading of the seismic phenomenon.
It is therefore necessary to model a multi-variable study in which to consider morphological, economic, and social variables inherent in the degree of technology, mobility, and infrastructure available and cross-reference them with the characteristics of the hazard. At a preventative level, some threats, including seismic hazard, can certainly be modelled.
On the basis of this model, starting from the real data collected by the INGV, an analysis was carried out to verify the spatial patterns that highlight the possible presence of scattering, i.e., areas with a greater concentration of high-magnitude seismic events, thus identifying a sort of spatial take off. As illustrated in the following paragraph, this analysis was developed through concentration maps, spatial ellipses of standard deviation, Ripley’s K function, and a search for spatial autocorrelations based on the Moran index.
3. Data and Methods
Low seismicity occurred west of Solfatara (as seen in
Figure 1) in 2005–2023. At depths of 1.5–3.0 km between Pozzuoli and Solfatara, only isolated and low-energy events have occurred since 2005, in contrast to the pronounced seismicity in 1982–84 that included duration magnitudes greater than 3. This difference has been accentuated by the cluster of events developing at depths of 2–3 km beneath Solfatara since 2005. This cluster appears to be slightly shifted to the NE compared with its position in 1982–84 [
27].
This analysis was based on the comparison between a historical series of data collected by stations (
Figure 2) of the INGV—National Institute of Geophysics and Volcanology—Vesuvius Observatory section, from 2005 to 2022 (
Figure 3a), a data set for the period 2022–2024 (excluding the month of May 2024) (
Figure 3b), in which the bradyseismic phenomenon restarted, and a further data set relating to the month of May 2024 alone (
Figure 3c), very important for the frequency and intensity of the seismic events (the INGV—National Institute of Geophysics and Volcanology, Vesuvius Observatory section, shares for free its monitoring data on its official website:
https://terremoti.ov.ingv.it/gossip/flegrei/index.html (accessed on 15 October 2024).
Spatial Pattern Real Data Analysis
The heatmap analysis already highlighted an initial dichotomous evidence of concentration/dispersion in relation to the comparison between the values of the month of May 2024 and those of the 2005–2022 historical series (
Figure 4a,b). In fact, the historical data presented a much wider dispersion.
We then proceeded with the SDE—Standard Deviation Ellipse—analytical–spatial methodology in order to verify the evidence from the heatmaps. SDE is a spatial analysis tool used to describe the geographic distribution of a set of points [
28]. This methodology allows us to understand the orientation, dispersion, and main direction of the spatial distribution of data. The main components of the SDE are the following:
- -
The centroid: The central point of the data distribution. It is calculated as the average of the x and y coordinates of all points;
- -
The principal axes: The major and minor axes of the ellipse represent the direction of maximum and minimum dispersion of the data. The major axis is the direction along which the data are most dispersed, while the minor axis represents the direction of least dispersion;
- -
The standard deviation: It indicates the dispersion of the points around the centroid. The standard deviation along the major and minor axes is used to define the shape of the ellipse.
The calculation of the centroid, the standard deviations along the major and minor axes, and the orientation angle is performed according to the following:
where
Cov(
x,
y) is the covariance between the
x and
y coordinates.
The results of the SDE (
Figure 5) clearly showed a confirmation of the spatial dynamics of the events, with a greater dispersion of the data of the 2005–2024 historical time series compared to those of the last two years 2023 and 2024, when recovery of the bradyseismic phenomenon occurred. In particular, then, the month of May 2024 revealed a further concentration of the data.
To complete the SDE analysis, Ripley’s K was evaluated in the compared scenarios.
Ripley’s K is a geostatistical tool used to analyze the spatial distribution of points in an area and understand whether they are randomly distributed, clustered, or dispersed. Specifically, Ripley’s K function compares the observed distribution of points to a random distribution to determine whether there are significant patterns of clustering or dispersion.
The calculation of the function
K is based on the following expression:
where
λ is the average dot density (number of dots divided by the total area);
N is the total number of points;
dij is the distance between i and j;
I(dij < r) is a function equal to 1 if dij < r or 0 otherwise.
To facilitate the interpretation, a normalized version of the
K function called Ripley’s
L function is often used:
where
L(r) = 0 indicates a random distribution;
L(r) > 0 indicates aggregation;
L(r) < 0 indicates dispersion.
In the specific case, the comparison between the image in
Figure 5 relating to the month of May 2024 and that relating to the 2005–2024 historical series showed a non-random concentration, which, however, confirmed the greater dispersion highlighted by the SDE in the case of the historical series (
Figure 6a,b).
Finally, we proceeded with a research analysis of any spatial autocorrelation. Spatial autocorrelation is a measure that evaluates the degree of similarity between values of a variable at spatially close points [
29]. In other words, it measures how similar the values of a variable observed in one location are to those observed in nearby locations. Spatial autocorrelation is critical for understanding spatial patterns in geographic data and for validating spatial models. One of the most common methods for measuring spatial autocorrelation is the Moran index [
30].
The Moran index (
I) is a statistic that quantifies spatial autocorrelation in a geographic data set. It is defined as follows [
31]:
where
N is the total number of spatial units;
xi is the value of the variable of interest in the georeferenced position i;
is the average of the values regarding the variable of interest;
Wij is an element of the contiguity matrix W that defines the proximity between the locations i and j;
W it is the sum of all elements.
Regarding the interpretation of the Moran index:
I > 0 indicates a positive spatial autocorrelation (a tendency of similar values to be close to each other);
I < 0 indicates a negative spatial autocorrelation (a tendency of dissimilar values to be close to each other);
I = 0 indicates absent spatial autocorrelation (the values randomly distributed in space).
In the search for spatial autocorrelation, it is useful to use the LISA (Local Indicators of Spatial Association) method, as it allows the identification of spatial clusters and local outliers [
32]. Unlike the global Moran index, which provides a single measure of spatial correlation for an entire dataset, the LISA method allows one to examine local variations in spatial correlation, offering a more detailed view of spatial patterns. Through the LISA method, the following elements are identified:
- -
Clusters: areas where similar values are close to each other (high local positive autocorrelation).
- -
Outliers: areas where dissimilar values are close to each other (high local negative autocorrelation).
- -
Significance: statistical evaluation of the significance of the identified clusters and outliers; in other words, significance tests (z-score tests or Monte Carlo simulations) are performed through appropriate software to determine whether the index values are significantly different from those that would be expected for a random distribution.
In details, the spatial weights matrix (WW) used in our analysis was based on inverse distance weights. Specifically, we defined , where dij represents the Euclidean distance between the locations ii and jj, with a distance threshold of 1 km to limit the interactions to geographically proximate locations.
The weights were row-standardized to ensure that the sum of the weights for each row equaled 1, which is a common practice to account for differences in the number of neighbors.
The significance of the Moran index was tested using 999 permutations to compare the observed value of the index against a reference distribution under the null hypothesis of spatial randomness.
The results obtained from the analysis of real data showed, through the Moran index, a strong spatial correlation both for the time series data and, especially, for the data of the month of May 2024 (
Figure 7a,b).
Similarly, the LISA method highlighted a strong clustering in both data sets, characterized only by homogeneous conditions of the high–high or low–low type and not by crossed conditions of the high–low type or the opposite type. In particular, the data of the month of May 2024 clearly identified some restricted geographical areas of the high–high type that require monitoring and attention from the relevant bodies and can be considered as additional data in the ongoing analyses of the domain related to other disciplines.
This scenario is visible in the comparison between the clustering histograms (
Figure 8a,b), in the comparison between the connectivity graphs (
Figure 9a,b), and finally, in the comparison between the cartographic representations of the LISA method (
Figure 10a,b).
4. Conclusions
The bradyseismic phenomenon currently underway in the Campi Flegrei area constitutes a long-term dynamic; however, in the current resumption of its activity, it manifests a specific spatial concentration of high-magnitude values, which contrasts with the more dispersed distribution of the past.
The proposed research aimed to represent a geographical reading/perspective, through spatial analysis tools, of the ongoing phenomenon in order to make complementary information available to scientists and experts capable of making the seismic picture in their possession clearer, based on geophysical and geological studies.
The results of these approaches provided a spatial framework in which some areas can be better monitored, as they are characterized by high-magnitude events and strong spatial autocorrelation. This evidence suggests integrating these conclusions with the plans for defining the risk areas and with the evidence provided by the subjects responsible for surveying the state of infrastructures and buildings in some particular areas (such as, for example, those in
Figure 11). The strong impact on the socio-economic activity of the area, in terms of the abandonment of various homes by the residents and the closure of economic, commercial, and entertainment activities, requires the use of interdisciplinary approaches that can broaden the analysis framework to facilitate the choices of political decision makers. The proposal of this geographical approach, complementary to existing approaches, can increase the accuracy of investigations, improve the investigation resolution of the phenomenon, and limit critical areas to more restricted geographical areas.
4.1. Implications for Policy Makers and Practitioners
The geographical approach could be seen as a method of thinking and problem solving that integrates and organizes all relevant information in the crucial context of a location. Policy makers can use this type of approach, or directly make use of its results, to identify patterns and trends, model scenarios and solutions, and finally make well-founded strategic decisions.
The idiographic nature of the geographical approach, which adopts the methods of field research, presents the risk of producing models with high uncertainty due to the fact of neglecting elements of interdependence between the causal variables of the seismic phenomenon. However, in its declared ancillary purpose related to the earth sciences, in this specific field regarding seismic dynamics, the feared risk does not constitute a critical factor.
In this sense, the added value of spatial tools, such as the one proposed, must be framed in the broader framework of system integrator actions that must be carried out by policy makers to quickly address such impactful emergencies.
4.2. Implications for Civil Protection
Civil protection agencies, which are at the forefront of disaster response and risk mitigation, require innovative tools and methodologies to anticipate and manage the impacts of seismic events. One such innovation is the complementary geographical approach to earthquake monitoring, which integrates multiple spatial and geophysical parameters to enhance the accuracy and reliability of seismic risk assessments.
Civil protection agencies are responsible for safeguarding communities from natural disasters. Their work includes some different activities such as monitoring seismic activity and identifying high-risk areas, developing early warning systems to reduce the loss of lives, designing mitigation strategies to protect critical infrastructures, and coordinating emergency response efforts during and after seismic events.
To perform these tasks effectively, agencies rely on robust monitoring systems. However, the traditional earthquake monitoring methods often face limitations, such as gaps in data coverage, inadequate integration of diverse datasets, and challenges in capturing the complexity of local and regional tectonic behaviors.
This proposed complementary geographical approach seeks to address these limitations by combining diverse datasets and techniques within a unified framework. Particularly, by adopting this approach, the civil protection agencies could optimize the integration of multidimensional data, combining seismic records, geophysical measurements, satellite imagery, and socio-geographic data, to create a holistic understanding of earthquake dynamics. Furthermore, they could emphasize a localized focus, tailoring the monitoring systems to reflect the unique geographical and tectonic characteristics of each monitored area. For example, urban neighborhoods may require highly localized monitoring to detect micro-scale precursors of seismic activity. Finally, the protection agencies, through complementary geographical spatial and temporal analyses, based on tools like GISs and spatial autocorrelation indices (e.g., Moran index), could identify patterns of clustering or anomalies over time and space. Therefore, by adopting a complementary geographical approach, civil protection agencies can improve risk assessment with enhanced spatial precision, enabling the identification of micro-zones with heightened seismic vulnerability, particularly in densely populated urban areas.
Support resource allocation through geographical analysis helps prioritize resources for high-risk areas, ensuring that preparedness and mitigation efforts are targeted and cost-effective, with less impacts on the citizens of the monitored areas. Moreover, a data-driven understanding of seismic risks will empower agencies to communicate effectively with communities, fostering trust and encouraging preparedness measures.