Next Article in Journal
Arctic Sea Ice Albedo Estimation from Fengyun-3C/Visible and Infra-Red Radiometer
Next Article in Special Issue
Calibration of SAR Polarimetric Images by Covariance Matching Estimation Technique with Initial Search
Previous Article in Journal
A Gaussian Function Model of Mesoscale Eddy Temperature Anomalies and Research of Spatial Distribution Characteristics
Previous Article in Special Issue
A Hierarchical Fusion SAR Image Change-Detection Method Based on HF-CRF Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On the Initial Phase of the Ongoing Unrest at Campi Flegrei and Its Relation with Subsidence at Vesuvio (Italy)

by
Antonella Amoruso
1,*,
Adriano Gualandi
2,3 and
Luca Crescentini
1
1
Department of Physics, University of Salerno, 84084 Fisciano, Italy
2
Department of Earth Sciences, University of Cambridge, Cambridge CB2 3EQ, UK
3
Osservatorio Nazionale Terremoti, Istituto Nazionale di Geofisica e Vulcanologia, 00143 Rome, Italy
*
Author to whom correspondence should be addressed.
Remote Sens. 2024, 16(10), 1717; https://doi.org/10.3390/rs16101717
Submission received: 28 March 2024 / Revised: 5 May 2024 / Accepted: 10 May 2024 / Published: 12 May 2024

Abstract

:
The densely inhabited area of Naples (Italy), between the Campi Flegrei and Vesuvio volcanoes, is one of the most hazardous regions in the world. After two decades of sustained subsidence, Campi Flegrei has been experiencing an accelerating uplift since 2005. The uplift is currently associated with unusual seismicity and increased degassing. To try to identify the cause of the shift from subsidence to uplift and explore any connection between Campi Flegrei and Vesuvio, we analysed the ground displacement time series of the two volcanoes from 1993 to 2010, obtained from ERS/ENVISAT Synthetic Aperture Radar imagery. To distinguish between the various sources of deformation, we used simple scatter plots and a blind source separation technique called variational Bayesian independent component analysis (vbICA). We obtained consistent results using both approaches. Specifically, with vbICA, we identified two significant independent components (ICs). IC1 describes the subsidence that occurred at Campi Flegrei prior to 2000, including the mini-uplifts of 2000 and 2005, and part of the post-2005 uplift. The expansion and contraction of two volumes beneath Campi Flegrei satisfy IC1: a sill-shaped volume at a depth of approximately 3 km and a small volume at a depth of 1–2 km, respectively. The two sources of deformation reproduce the large-scale deformation in the Campi Flegrei area and the local deformation in the Solfatara area, respectively. In the Campi Flegrei area, IC2 exhibits primarily uplift, which is concentrated in the eastern part of the caldera. The deformation pattern is complex and difficult to interpret. If we model it using simple spheroidal deformation sources, the pattern suggests that two volumes at depths of approximately 9 and 8 km are experiencing opposite activity, namely contraction (beneath the southwestern part of the caldera) and expansion (beneath the central part of the caldera). In the Vesuvio area, IC2 is consistent with the deformation induced by the contraction of a volume at a depth of around 9 km. The contraction beneath Vesuvio is smaller in magnitude than the expansion/contraction beneath Campi Flegrei. The correlation observed after 2002 between uplift at Campi Flegrei and subsidence at Vesuvio suggests the transfer of magma and/or magmatic fluids between the two plumbing systems at 8–9 km depth. This implies that part of the ongoing unrest at Campi Flegrei may have been promoted by mass transfer from below Vesuvio.

Graphical Abstract

1. Introduction

The deformation style of volcanoes can be detected at the surface using geodetic techniques, primarily global navigation satellite systems (GNSSs) and satellite synthetic aperture radar (SAR) [1]. Ground deformation can provide information on the kinematics of both the aquifer and the entire plumbing system at various levels and depths, as well as regional tectonics. The study of the initial stages of subsidence or uplift, or the transition between the two, can reveal weak processes that may not be detected when deformation data are dominated by the effects of stronger or shallower processes. Focusing on this transition can also highlight the diagnostic features of the processes governing it.
The Campi Flegrei caldera and Vesuvio volcano are two active and high-risk volcanoes located ∼25 km apart from each other. They are situated on opposite sides of the densely populated city of Naples, Italy (Figure 1). Over the past ∼15 ka, several dozen eruptions have occurred at Campi Flegrei, accompanied by the resurgence of the centre of the caldera. After approximately 3 ka of quiescence and centuries of subsidence, the last eruption occurred in 1538, preceded by a period of increasing seismicity and uplift [2]. Magma was nearly erupted again between 1540 and 1582, after which post-eruptive deflation began [3].
Campi Flegrei has experienced intermittent unrest since the 1950s, with four main episodes occurring between 1950–1952, 1969–1972, 1982–1984, and 2005 to the present [4,5]. The period of unrest between 1982 and 1984 was followed by prolonged subsidence, which was interrupted by a few shorter phases of uplift (referred to as mini-uplifts) until 2004. From 2005 to the present, there has been an almost continuous uplift, with an accelerating trend. As of January 2024, the ground level at Rione Terra (R in Figure 1), which is the area of maximum uplift, is 25 cm higher than it was at its 1984 peak and 120 cm higher than it was in 2005 [6] (Figure 2).
Figure 1. Map of the Campi Flegrei and Vesuvio areas, UTM WGS84 33N coordinates. Palette colours indicate elevations [7], cyan areas are bodies of water. Black dots and labels: R, Rione Terra; S, Solfatara; P, Pisciarelli. The red rectangle encloses the area referred to in this paper as the “Campi Flegrei area”. The orange and magenta curves indicate the rims of the Campi Flegrei caldera collapses of the Campanian Ignimbrite (about 39 ka) and the Neapolitan Yellow Tuff (about 15 ka) eruptions [8]. The blue curve indicates the rim of the collapsed caldera of an older structure, Mt. Somma, with the Vesuvio cone rising from the centre.
Figure 1. Map of the Campi Flegrei and Vesuvio areas, UTM WGS84 33N coordinates. Palette colours indicate elevations [7], cyan areas are bodies of water. Black dots and labels: R, Rione Terra; S, Solfatara; P, Pisciarelli. The red rectangle encloses the area referred to in this paper as the “Campi Flegrei area”. The orange and magenta curves indicate the rims of the Campi Flegrei caldera collapses of the Campanian Ignimbrite (about 39 ka) and the Neapolitan Yellow Tuff (about 15 ka) eruptions [8]. The blue curve indicates the rim of the collapsed caldera of an older structure, Mt. Somma, with the Vesuvio cone rising from the centre.
Remotesensing 16 01717 g001
Seismicity accompanied ground deformation at Campi Flegrei during the 1982–1984 unrest. Seismic activity nearly vanished during the two decades of slow subsidence and mini-uplifts following 1984 and resumed after the onset of the current unrest. Seismic activity persisted in 2012–2014 and has since continued to increase, although the average rate remains lower than during 1982–1984 (e.g., [9,10,11]). The current distribution of earthquakes closely resembles that recorded during the uplift period of 1982–1984 and is mostly limited to the first 3 km in depth [6]. Seismicity in the caldera is attributed, at least in part, to the injection of magmatic fluids from a shallow reservoir and/or the deep magma chamber [12,13,14].
The Solfatara–Pisciarelli fumarole area (S and P in Figure 1) is currently experiencing high surface fluid release (e.g., [6,15]). Changes in fumarole gas composition and temperature, as well as surface gas release, indicate heating of the hydrothermal system and massive magma outgassing in the deep part of the plumbing system (e.g., [12,16,17]).
Figure 2. Changes in ground level at Rione Terra (R in Figure 1), which is the area of maximum vertical displacement at Campi Flegrei. (a) From 1905 to the end of 2023. (b) Period investigated in this work. Blue points, levelling data [4]; red points, vertical displacement from GNSS data [18]; green points, temporal function of the first independent component obtained from the analysis of SAR data, reversed in sign and properly scaled.
Figure 2. Changes in ground level at Rione Terra (R in Figure 1), which is the area of maximum vertical displacement at Campi Flegrei. (a) From 1905 to the end of 2023. (b) Period investigated in this work. Blue points, levelling data [4]; red points, vertical displacement from GNSS data [18]; green points, temporal function of the first independent component obtained from the analysis of SAR data, reversed in sign and properly scaled.
Remotesensing 16 01717 g002
The study of the relationship between uplift and gas outflow has extensively utilised various geochemical species and their ratios, such as the CO2/H2O ratio, which indicates the relative contributions of magmatic and hydrothermal sources [16]. The sudden increase in the CO2/H2O ratio in 2000 and the significant rise in the Pisciarelli fumarole flow rate and temperature in 2005–2006 are noteworthy for this study.
Campi Flegrei’s large-scale deformation can be explained almost completely by the expansion/contraction of a sill-like volume centred off Rione Terra at about 3–4 km depth (e.g., [19,20,21]); the location is consistent with the coda-wave high-attenuation patch at 3–4 km depth detected beneath Rione Terra [22]. The sill may receive magma and/or magmatic fluids from the partially molten volume detected beneath the caldera at depths greater than 7.5 km (e.g., [23]). The sill contracts during periods of subsidence. The residual deformation is mostly limited to the Solfatara area and can be attributed to a local source, which can be explained by hydrothermal-driven poroelastic deformation near the point where hot magmatic fluids are injected into the hydrothermal system [24]. The two paired large-scale and local deformation sources are consistent with seismic (e.g., [22]) and geochemical (e.g., [12]) data. Furthermore, the connection between sill activity and seismicity, particularly the impact of transient stress rates on fault slips, has been highlighted [25]. However, the hydrothermal system may still play a role and interact with the seismic activity (e.g., [26]).
Somma-Vesuvio (Figure 1) is a stratovolcano with a summit caldera (Mt. Somma) and a recent cone (Vesuvio) resulting from several Plinian eruptions (e.g., [27]). The most recent Plinian eruption took place in 79 AD, and was followed by sub-Plinian eruptions in 472 and 1631 AD. The latter was followed by semi-persistent activity that lasted until 1944.
Vesuvio’s eruptions in the last 2000 years were fed by a magma reservoir approximately 4 km deep (e.g., [28,29]), although there is currently no evidence of stored melt [30]. A large sill-like magma chamber is located at approximately 8 km depth [31,32]. Between 2002 and 2010, ground velocity maps indicate contraction at 7–8 km depth [33]. Seismic activity at depths greater than about 2 km significantly decreased around 2002 (e.g., [34]).
Seismic tomography [31] and reflections [23], as well as textural and geochemical investigations [35], indicate the presence of a single partially molten layer with a top at approximately 8 km depth beneath both Campi Flegrei and Vesuvio. The percentage of melt varies from around 80% below Campi Flegrei [23] to around 10% below Vesuvio [31,32]. Gravity modelling [36] suggests the presence of several magma pockets at depths greater than approximately 8 km.
There is a long-standing controversy over the nature—magmatic, hydrothermal, or hybrid—of the unrest episodes at Campi Flegrei. The determination of intrusion density is reliable only when considering the effects of heterogeneities in elastic parameters on ground displacements and gravity variations [37]. These effects also depend on the shape of the source and the physical and chemical characteristics of the intruded fluid, such as compressibility. Whatever the nature of the Campi Flegrei unrest episodes, they appear to belong to an evolutionary sequence of long-term perturbations of the magmatic system (e.g., [38]); possible future behaviours of the caldera range from a deceleration of ground motion to phreatic explosions or an upwelling of magma from the deep part of the plumbing system.
This work analyses the transition from subsidence to uplift at the Campi Flegrei caldera, which occurred between 2000 and 2005. The ground deformation at Campi Flegrei is characterised by two distinct patterns. The first is related to the 3–4 km deep sill, while the second dominates deformation at the caldera borders and is complex and difficult to interpret. Modelling the second pattern using simple spheroidal deformation sources suggests that two volumes at depths of approximately 9 and 8 km were experiencing opposite activity, namely contraction and expansion. Contraction occurred beneath the southwestern part of the caldera, while expansion occurred beneath the central part of the caldera. The transition of Campi Flegrei from subsidence to uplift coincided with the onset of subsidence at the nearby Vesuvio volcano. The Vesuvio subsidence is consistent with contraction occurring at a depth of approximately 9 km beneath Vesuvio itself.

2. Materials and Methods

2.1. SAR Data and Analysis

Line-of-sight (LOS) ground displacement time series were obtained by the Istituto per il Rilevamento Elettromagnetico dell’Ambiente del Consiglio Nazionale delle Ricerche (IREA-CNR), Naples, Italy, from ERS-ENVISAT SAR images using the DInSAR small baseline subset (SBAS) technique [39,40,41]. These time series have already been used in the past by us and other researchers to study the ground deformation at Campi Flegrei and Vesuvio (e.g., [19,42,43]); they include 138 ascending-orbit images from January 1993 to September 2010 and 155 descending-orbit images from June 1992 to September 2010 (Figure 3). The temporal interferometric coherence of the time series is greater than 0.55 for the ascending orbit and 0.5 for the descending orbit. The accuracy of the displacement time series is approximately ±0.5 cm [44].
LOS displacement data are irregularly located, spaced approximately 50 to 80 m apart in covered areas. To decrease the amount of data and noise, we estimate the time series of LOS displacements on a regular (150 m × 150 m) UTM WGS84-33N grid by computing a median value for every non-empty block in the grid. We also mask out the Vomero area within the city of Napoli and five spots along a quasi-circular “girdle” at about 10 km around Vesuvio because they are subjected to local anthropogenic subsidence [33,41].
The satellite revisit time is 35 days, and the ascending-orbit images are acquired 6 days after the descending-orbit images. The time series of ascending and descending orbits have non-coincident missing times; to have compatible time series, we use displacement data only for paired images of ascending and descending orbits, i.e., when the available ascending- and descending-orbit data are separated by 6 days (Figure 3).
To determine the direction of the LOS per pixel in the ascending and descending orbit images, we first use the European Space Agency Sentinel Application Platform (ESA SNAP) to obtain the incidence angle (IA) per pixel with respect to the reference ellipsoid ( θ a and θ d for the ascending and descending orbits, respectively). We then calculate the average angle of the horizontal projection of the LOS from grid east ( β a and β d , respectively, positive clockwise) from the coefficients of the best-fit plane to the IA. We obtain
θ a = 22.94 + 6.93 × 10 5 ( x 437481 ) + 1.39 × 10 5 ( y 4522079 ) θ d = 21.62 7.10 × 10 5 ( x 437481 ) + 1.56 × 10 5 ( y 4522079 ) β a = arctan ( 1.39 / 6.93 ) = 11.3 β d = arctan ( 1.56 / 7.10 ) = 12.4
where all angles are in degrees and ( x , y ) are the pixel coordinates.
At each pixel, a plane is defined by the ascending and descending LOS. The unit vector n ^ , which is perpendicular to this plane, is obtained by taking the cross product of the unit vectors of the two LOS, n ^ a and n ^ d . Rotating the horizontal projection of n ^ by 90 around the vertical axis gives the unit vector m ^ . This lies at the intersection of the plane of the LOS and the horizontal surface, and consequently is used to define the horizontal motion u H . The unit vector p ^ , which is perpendicular to m ^ and lies on the plane of the LOS, is obtained by taking the cross product between n ^ and m ^ . Its direction is almost vertical, and we use it to define the quasi-vertical displacement u V (Supplementary Information Figure S1, modified from [45]).
In equations, the displacements u a and u d along the two LOS are given by
u a u d = n ^ a · m ^ n ^ a · p ^ n ^ d · m ^ n ^ d · p ^ u H u V
where u H is the projection of the ground displacement along m ^ ; u V is the projection of the ground displacement along p ^ ; and “·” denotes the dot product. The horizontal displacement u H and the quasi-vertical displacement u V are obtained by inverting Equation (2). The quasi-vertical displacement u V would actually be vertical if the azimuths of the two LOS were +/− 90°.
By using the direction cosines of the two LOS relative to the east, north, and up directions and following the above calculations, when the azimuths of the horizontal lines perpendicular to the two LOS— α a and α d for the ascending- and descending-orbit, respectively—are equal in magnitude but opposite in sign, i.e., | α a |   =   | α d |   =   α , we find
u H = cos ( θ a ) u d cos ( θ d ) u a / sin 2 ( α ) sin 2 ( θ a θ d ) + cos 2 ( α ) sin 2 ( θ a + θ d ) , u V = cos 2 ( α ) sin ( θ a + θ d ) sin 2 ( α ) sin ( θ a θ d ) sin ( θ d ) u a + + cos 2 ( α ) sin ( θ a + θ d ) + sin 2 ( α ) sin ( θ a θ d ) sin ( θ a ) u d / sin 2 ( α ) sin 2 ( θ a θ d ) + cos 2 ( α ) sin 2 ( θ a + θ d ) × × 1 ( sin 2 ( α ) cos ( θ a θ d ) + cos 2 ( α ) cos ( θ a + θ d ) 2 .
Since the grid convergence (i.e., the angle between true north and grid north, positive clockwise) in the study region is about 0.47 , we approximate α a = 12 and α d = 12 . Thus, we can compute u H and u V using Equation (3) with α = 12 . As the two LOS have similar angles of incidence, u H is pointing approximately east; the angle of u V from the vertical is less than 5 .

2.2. vbICA

In principle, different sources can contribute to the observed ground deformation. For this reason, we cast the problem as a blind source separation (BSS) one, where one attempts to disentangle the different sources from the knowledge of the observed mix. Independent component analysis (ICA) is a commonly used approach [46]. A clear advantage in using a decomposition technique like ICA when studying SAR measurements stems from the fact that it offers the chance to reduce the dimensionality of the dataset. In other words, we can extract only a small number of relevant independent components (ICs) to describe the ground deformation. This helps because we can visually inspect the time evolution of a few components instead of investigating thousands of time series, one per pixel, resolved by the SAR satellites. Other approaches to reduce the dimensionality of the data exist, and a popular and intuitive one is principal component analysis (PCA). The advantage of ICA with respect to PCA is that it generalises it, searching for components that are not only statistically uncorrelated but also statistically independent. The main assumption of ICA is rooted in information theory: the best way to separate a signal made of a mix of different sources is to minimise the mutual information between the components. The standard ICA decomposition relies on two assumptions: (1) the sources are linearly mixed (i.e., the mixing is described by a matrix A); (2) the sources are not moving, so the mixing matrix A is constant. Under these assumptions, the problem consists of finding the source matrix S that generates the observations X:
X A S + N
with each source vector of S being statistically independent from the others and N some possible additional Gaussian noise.
The data matrix can be organised either in a spatial or temporal mode. With this, we mean that we can search for independent components in either the temporal or the spatial domain. If X contains M rows, with M being twice the number of pixels in the images (because we have both u H and u V ), and T columns, with T being the total number of observed epochs, we say that the decomposition is performed in temporal mode (T-mode). In fact, in this case, the sources (which are the rows of S) will be statistically independent time series. If instead we organise X as a T × M matrix, the decomposition is performed in spatial mode (S-mode). Despite both having their pros and cons, the S-mode is the most commonly used to study SAR time series because of the abundance of spatial information (thousands of pixels) with respect to the temporal one (tens or hundreds of epochs). Here, we perform the decomposition in S-mode. This means that each row of S corresponds to a spatial distribution or map, and the goal is to find a finite number of rows that are statistically independent from each other and that can describe the data if appropriately linearly mixed. The columns of the mixing matrix A represent the weight of each component in the mix at every given time step, and thus can be described as temporal functions.
The most commonly used ICA algorithm is the FastICA one [47], but it presents two major problems when dealing with geodetic data: (1) it does not allow us to deal with missing data; (2) it uses approximations of the kurtosis and skewness that assume a unimodal distribution of the sources. The variational Bayesian ICA (vbICA) is an alternative algorithm that uses a modelling approach instead of a mapping approach. It is a generative model that employs Bayesian inference to determine the best probability density function (pdf) of the sources and mixing matrix. It offers more flexibility than FastICA because it models the pdf of the sources with a mix of Gaussians, thus allowing us to describe multimodal distributions.
In practice, the vbICA model is defined by a set of parameters, W. These describe the sources S via the aforementioned mix of Gaussians and their mixing proportions, the mixing matrix A and the additional noise N. For each Gaussian in the mix, designated by an indicator variable q, we need to specify a mean ( μ ) and a precision ( β ), as well as the probability π of being selected to contribute to the explanation of the source. The independence of the sources is realised by assuming that p ( S ) = i = 1 L p ( s i ) , with L being the number of components. Other assumptions imply a Gaussian distribution of the columns of the mixing matrix A and a zero-mean Gaussian additive noise, for which the precision Λ needs to be estimated. The parameters are thus indicated by W = { S , q , μ , β , π , A , Λ } . The goal is to find the posterior pdf for each of these parameters. One possibility is to estimate this via Monte Carlo integration, but a more efficient way is to use a variational approximation. An approximating pdf for W is introduced, and the Kullback–Leibler divergence between this approximating pdf and the true posterior is minimised by maximising the model negative free energy (NFE). The choice of the approximating pdf closely follows the work of [48], with all parameters being independent except the parameters q and S, for which the joint distribution is considered. Another advantage of vbICA over more traditional approaches to the BSS problem consists of the fact that the variational Bayesian framework enables us to handle missing data. The original algorithm [48] has been adapted to the study of geodetic time series (GNSS [49] and InSAR [50]), showing better results when compared to FastICA [47].
To help find the best ICA solution, it is common to whiten the data (i.e., to attribute to each principal component the same variance) using results from a PCA (e.g., [51]) or the singular values of the data variance–covariance matrix, which can be obtained from a singular value decomposition of X. As in [50], as a first step, we apply vbICA without whitening, and select the number of components to retain using the automatic relevance determination method (ARD [52]). In practice, each column i of the mixing matrix A is described by a Gaussian random variable for which the mean prior is set to zero, and the precision is controlled by a parameter α i . The posterior of the precision will be an adaptation of the prior precision. The idea behind ARD is that large values of α i will tend to suppress the contribution of the i-th component. Thus, we monitor the ratio between the minimum and maximum variance associated with each column to estimate the number of components supported by the data. Finally, we apply vbICA again, using only the chosen number of components and performing the whitening preprocessing step. The order in which vbICA provides the independent components is not significant.
Each temporal function starts at zero, and the value corresponding to the maximum of its modulus is non-negative. The full range of variation, i.e., the difference between the maximum and minimum of the temporal function, is 1. Therefore, the spatial distribution of each IC provides the per-pixel displacement related to the period during which the maximum variation in the temporal function occurred. One can obtain the temporal evolution at each pixel by simply multiplying the spatial distribution by the corresponding temporal function.

2.3. Mathematical Modelling of Deformation Sources

Realistic, accurate modelling of ground deformation is virtually impossible due to the complexity of the physical processes and medium. Nonetheless, the spatial distribution of the two ICs (i.e., the rows of the matrix S from Equation (4)) shows some remarkable coherent structure, and their spatial patterns can be interpreted in terms of the ground displacement caused by geometrically simple pressurised cavities. Before providing some modelling details, we want to dispel any doubt about the use of the term “source”. On the one hand, when we solve the blind source separation (BSS) problem, each IC corresponds to a source that is separated from the other sources by virtue of their statistical independence. On the other hand, each IC may be the result of multiple physical processes or sources (in our case, pressurised cavities). For example, as it will be clearer later, we need multiple physical sources to explain the spatial distribution associated with IC2.
The pressurised cavities used to model the surface deformation are either very small (point-like) or finite spheroids with a vertical polar axis and a variable aspect ratio (polar-to-equatorial radius). Rock heterogeneity and topography are known to have an impact on surface deformation. At Campi Flegrei, topography does not appear to be a relevant factor. At Vesuvio, which has non-negligible relief (almost 1300 m a.s.l.), it could potentially play a more important role; however, ground displacement data are mainly limited to regions with relatively smooth topography.
The presence of soft superficial layers above a deformation source can make it appear shallower, known as the “focusing effect” [53,54]. Soft superficial layers also increase the ratio of horizontal to vertical displacement [37]. Seismic tomography indicates that a layered model is a suitable first-order approximation of the true heterogeneity of rocks, at least below Campi Flegrei (e.g., [23]). Therefore, we calculate ground displacements using the technique of [19,55], which can be applied in both homogeneous and layered half-spaces. This technique represents an improvement upon that of [56], which employs Eshelby’s elastic inclusion theory [57] and Mindlin’s half-space point force solution [58] to compute approximate ground displacements through the use of three vector dipoles, which are three pairs of equal and opposite forces with the force and arm in the same direction (Supplementary Information Figure S2). The vector dipoles are aligned with the axis of the ellipsoid and have magnitudes that are calculated in terms of elliptic integrals as a function of the ellipsoid axis ratios. In this work, the group of three vector dipoles is designated by the symbol M D .
In [19,55], the deformation caused by a pressurised ellipsoid is approximated by the deformation caused by seven point sources. The approximation involves a single central point source and six additional point sources, distributed symmetrically along the ellipsoid semi-axes a, b, and c, at distances a / 5 , b / 5 , and c / 5 from the ellipsoid centre. Each distributed source is represented by M = M D / 2 , while the central point source is represented by M 0 = 2 M D . In the case of a layered half-space, the PSGRN post-seismic deformation code [59] is employed to compute the displacement Green’s functions of a given layered half-space for M D . In this work, the technique has been enhanced by incorporating a free parameter, f, with a range of 0 f 1 , to determine the optimal distances of the six point sources from the ellipsoid centre; these distances are now f a , f b , and f c . Consequently, the distributed sources are represented by M = M D 1 / 10 f 2 , while the central point source is represented by M 0 = M D 1 3 / 5 f 2 (Supplementary Information Figure S2).
This study employs the Campi Flegrei elastic layered model of [54], which is based on the 1D Vp model of [60] and the Nafe–Drake curve [61]. Poisson’s ratio is set to 0.25 at all depths to account for the drained response of the medium (see Table 1). To validate our computations for a finite spheroidal source, we used COMSOL Multiphysics® v. 5.2 and finite element method (FEM) modelling with a 2D axisymmetric domain. The domain extends 50,000 m horizontally from the sill axis and 50,000 m below the surface. We assigned zero normal displacements at the lateral and bottom boundaries of the domain, while the upper boundary (ground surface) is stress-free. A uniform pressure is assigned to the walls of the spheroidal cavity. The properties of the elastic medium within the domain are given in Table 1.

2.4. Inverting Spatial Patterns for Deformation Source Parameters

The source parameters that best fit the data, including the three-dimensional centre position, potency, aspect ratio, and, in the case of finite spheroids, equatorial radius and f, are determined by minimising the absolute deviation of residuals using the L1-norm misfit function. This approach is preferred as it is less sensitive to outliers than chi-square minimisation. Optimisation of the misfit function is achieved using adaptive simulated annealing (ASA) [62] and the neighbourhood algorithm (NA) [63]. ASA aims to find the best global fit of a non-linear non-convex misfit function on a multidimensional parameter space; this method ensures a broad global search in the early stages and ample rapid convergence in the final phases, self-optimising programme options recursively. NA generates ensembles of model parameters that preferentially sample the good data-fitting regions of the parameter space and is also used for global optimisation. As with any other Monte Carlo method, success in finding the optimal model parameters is never guaranteed by ASA or NA. Therefore, we cross-check dozens of ASA and NA inversions to increase confidence in the results. Increasing the number of free parameters typically makes minimising the misfit more challenging.
Parameter uncertainties are estimated in the form of marginal probability density functions (MPDFs) using NA Bayes (NAB) [64]. NAB draws quantitative inferences from the entire ensemble generated by NA and enables the measurement of resolution and trade-offs between parameter pairs within a Bayesian framework. It is important to note that finding true MPDFs is never guaranteed, and therefore we average MPDFs from five NAB runs to increase confidence.
The best-fitting parameters are independent of data uncertainties, as long as they are equal for all data. However, the MPDFs are highly sensitive to uncertainties in the data. Unfortunately, the actual uncertainties are unknown. To address this, we first assume the same level of uncertainty for all data and use ASA to obtain the best-fit parameters. Then, we estimate the MPDFs with NAB by using the mean absolute deviation of the ASA best-fit model as the data uncertainty.

3. Results

We analyse horizontal and quasi-vertical ground displacements from 1993 to 2010. Horizontal displacements ( u H ) are approximately directed from west to east; the angle of quasi-vertical displacements ( u V ) from the vertical is less than 5 ° . The available time series contains breaks, including four breaks lasting several months to over a year during the Campi Flegrei subsidence period. Additionally, there is a long break from late 2000 to mid-2002, after which only two data points are available until mid-2003, when the time series becomes approximately continuous (see Section 2).

3.1. Pre-Analysis of Displacement Data Using Scatter Plots

Scatter plots are created to sample both the central and peripheral areas of Campi Flegrei. The plots compare u V and u H in selected areas of the Campi Flegrei with the quasi-vertical displacement u V in the area of maximum 1993–2000 subsidence from 1993 to 2010. In general, when using the vertical displacement observed in the area of maximum uplift or subsidence as the independent variable, a positive correlation of the horizontal displacement to the east and a negative correlation to the west are expected. In areas with easting values close to those of the area of maximum uplift or subsidence, the horizontal displacement is expected to be almost zero.
All the scatter plots (Figure 4 and Figure 5 and Supplementary Information Figures S3–S36) confirm that the deformation field during the Campi Flegrei subsidence was stationary. The characteristics of the deformation field changed dramatically around 2003, after the long data break. The scatter plots pertaining to the central part of the caldera (Figure 4e and Supplementary Information Figures S12, S13 and S15–S18) indicate that a deformation source, similar to that of the subsidence period, persisted after subsidence had ceased. This source transitioned from contraction to expansion and was the primary cause of deformation in that area. This coexisted with a second deformation pattern that prevailed at the caldera borders. The eastern part of the caldera is associated with a significant uplift, as shown in Figure 4 and Figure 5 and Supplementary Information Figures S3–S36, due to this second deformation pattern.
Based on the scatter plots and the observed changes in the deformation field around Vesuvio during the same period [33], we apply vbICA to displacement data from 1993 to 2010, first in the Campi Flegrei and Vesuvio areas separately, to confirm the coeval change in deformation style in both areas, and then in a larger area including both volcanoes, trying to identify any processes that may have affected the whole area.

3.2. vbICA Results

3.2.1. Campi Flegrei Area

The ARD criterion suggests keeping two ICs based on the trend of the ratio of minimum-to-maximum posterior mixing matrix variances in relation to the number of independent components, not shown here for the sake of brevity. The temporal function for the first component strictly follows the vertical displacement at Rione Terra from levelling and GNSS data (see Figure 2 and Figure 6a). The spatial distributions of IC1 are very similar to those of the mean velocities during subsidence (1995–2000, [19]; 1993–2000, [42]). The temporal function retrieved for the second component increases at a steady rate (Figure 6d) after the data break, which extends from December 2000 to October 2002 (see Figure 3). The spatial distribution of u V for IC2 shows a distinct spot in the Solfatara area (Figure 6e). Otherwise, both u H and u V suggest that IC2 is associated with uplift that is primarily concentrated in the eastern part of the caldera. The palettes in Figure 6b,c,e,f indicate that the ratio of horizontal to quasi-vertical displacement ( u H / u V ) of IC2 is greater than that of IC1.

3.2.2. Vesuvio Area

Once again, the ARD criterion suggests keeping two ICs. The temporal function of IC1 indicates a decrease in ground displacement speed, which is almost exponential (Figure 7a). A westward horizontal displacement from Vesuvio and an eastward displacement towards Vesuvio can be observed in the spatial distribution of u H for IC1 (Figure 7c). This finding is in line with the average horizontal velocity in the Vesuvio area reported in [33] for the period 1993–2000. The spatial distribution of u V for IC1 suggests subsidence in the east and uplift in the west; its interpretation is not straightforward (Figure 7b). The temporal function of IC2 is more erratic but appears to increase steadily after 2001 (Figure 7d). The distributions of IC2 indicate subsidence (Figure 7e,f), which is consistent with the mean velocities over the period 2002–2010 [33].

3.2.3. Campi Flegrei and Vesuvio Combined

Since the focus is on processes that involve the entire area, including Campi Flegrei and Vesuvio, this section provides more detailed information than Section 3.2.1 and Section 3.2.2. Figure 8a,b show the ARD and explained variance for the first six components. The first component accounts for over 91% of the variance. According to ARD, either two or four components should be retained. In both cases, the temporal function of IC1 follows the vertical displacement at Rione Terra from levelling and GNSS data (see Figure 9). However, the general trend is in line with the decrease in ground displacement speed in the vicinity of Vesuvio, as evidenced by the analysis of the Vesuvio area alone (Figure 7a). The temporal function of IC2 remains almost constant until the data break, which extends from December 2000 to October 2002 (see Figure 3). After this point, it increases at a steady rate (Figure 9). The spatial distributions of IC1 and IC2 also remain almost unchanged when either two or four components are retained. This shows the robustness of the results.
The spatial distributions of IC1 and IC2 (Figure 10 and Figure 11) are consistent with those obtained by analysing the Campi Flegrei and Vesuvio areas separately (Figure 6 and Figure 7). In the Vesuvio area, u V for IC1 (Figure 10b) is very small except within the rim of the Mt. Somma caldera (Figure 1). The spatial distribution of u H for IC2 (Figure 11a) provides an overview of the entire area. It reveals a westward displacement in the area west of Campi Flegrei, an eastward displacement between Campi Flegrei and Vesuvio, and a westward displacement in the area east of Vesuvio. The spatial distribution of u V indicates subsidence at Vesuvio and uplift at Campi Flegrei, which are mostly concentrated in the eastern part of the caldera. These observations are consistent with expansion occurring below Campi Flegrei and contraction occurring below Vesuvio.

3.3. Campi Flegrei and Vesuvio Combined: Modelling of the Spatial Distributions of IC1 and IC2

The combined spatial distributions and temporal functions of IC1 and IC2 indicate an uplifted area in the central and eastern parts of Campi Flegrei after approximately 2003. This observation is supported by the scatter plots of u V and u H (Figure 4 and Figure 5, and Supplementary Information Figures S3–S36). The ground displacements associated with IC2 are small and somewhat noisy. Their interpretation is not straightforward. However, we attempt to model them mathematically using simple pressurised spheroidal cavities.

3.3.1. First Independent Component

As deformation related to IC1 is mainly limited to the Campi Flegrei area, we only invert the spatial distribution of pixel displacements ( u H and u V ) within the red rectangle in Figure 1. Following [19], we investigate whether the change in volume of two pressurised spheroids, one finite and the other point-like, both with a vertical polar axis, can account for the spatial distribution of the first component. The spheroids are embedded in a layered medium (see Table 1). The model has 12 free parameters, including the three-dimensional position, potency, and polar-to-equatorial radius of the point-like source, as well as the three-dimensional position, potency, polar-to-equatorial radius, equatorial radius, and f of the finite source. Table 2 presents the model parameters obtained from ASA and NA inversions, while Figure 12 shows the MPDFs of the model parameters estimated using NAB.
These results are in line with findings in [19]. The mean absolute deviation is 0.63 cm, which is consistent with the accuracy of the displacement time series [44]. Assuming a homogeneous medium, the mean absolute deviation increases to 0.78 cm. This supports the use of a layered medium.
Figure 13 shows the spatial distribution of IC1, model predictions, planimetric positions of the sources, and model residuals. The only visible discrepancy between the spatial distribution of IC1 and model predictions is related to u H in the Solfatara area. This residual anomaly can be attributed to two main reasons: (1) the inadequate approximation of the effects of fluid injection at the base of the hydrothermal system with a pressurised spheroid embedded in an elastic medium; and (2) the ineffectiveness of IC1 in capturing the entire displacement history at Solfatara. A local anomaly is also present at Solfatara in the IC2 spatial distribution (Figure 11b), as well as in the IC3 and IC4 spatial distributions (Figure S37). The Solfatara anomaly exhibits distinct characteristics in each of the ICs.
Figure S38 compares the radial and vertical superficial displacements caused by the finite spheroid, as computed using our approximate approach and FEM. The results demonstrate the effectiveness of our approximate semi-analytical approach (Section 2).

3.3.2. Second Independent Component

The spatial distribution of IC2 in the Vesuvio area is consistent with deep depressurisation, as previously noted by [33]. However, the spatial distribution in the Campi Flegrei area is more complex and its interpretation is not straightforward. We attempt to model it using simple point-like deformation sources. As a first step, we invert the spatial distribution of u H and u V in the whole area including Campi Flegrei and Vesuvio for two point-like sources: a pressurised spheroid beneath Vesuvio and a generic moment tensor [65] beneath Campi Flegrei. It is assumed that the medium is homogeneous because using the same layered model throughout the entire study area is not feasible. The model has 14 free parameters, including the three-dimensional position, potency, and polar-to-equatorial radius of the spheroid, as well as the three-dimensional position and six elements of the (symmetric) moment tensor. The mean absolute deviation is 0.53 cm. Inside the Campi Flegrei area (red rectangle in Figure 1), the mean absolute deviation is 0.67 cm.
As a second step, we test a model that involves three pressurised point-like spheroids. Two of these spheroids are located under Campi Flegrei, and one is located under Vesuvio. The model has 15 free parameters, including the three-dimensional position, potency, and polar-to-equatorial radius of each spheroid. The mean absolute deviation is now 0.48 cm. The mean absolute deviation within the Campi Flegrei area (highlighted in the red rectangle in Figure 1) is 0.52 cm.
The decrease in the mean absolute deviation within the Campi Flegrei area when transitioning from 14 to 15 degrees of freedom indicates that the three-source model is better supported by the data than the two-source model. Table 3 presents the model parameters of the three deformation sources obtained from ASA and NA inversions, while Figure 14 shows the MPDFs of the model parameters estimated using NAB.
Two volumes at depths of approximately 9 and 8 km were experiencing opposite activity, namely contraction and expansion. Contraction occurred beneath the southwestern part of the caldera, while expansion occurred beneath the central part of the caldera. Deformation in the Vesuvio area is induced by contraction at a depth of approximately 9 km. The contraction beneath Vesuvio was smaller in magnitude than the expansion/contraction beneath Campi Flegrei.
Figure 15 displays the spatial distribution of IC2, model predictions, planimetric positions of the sources, and residuals. The model successfully accounts for the fundamental characteristics of the IC2 spatial distributions. Due to its reliance on deep sources embedded in a homogeneous half-space, the ground displacement field computed by the model is necessarily very smooth. As a result, the model is unable to account for local features.

4. Discussion

The dynamics of the Campi Flegrei caldera may be influenced by various processes that occur in different parts of the plumbing and hydrothermal systems, each with its own time scale and surface expression (e.g., [13,15,35,42,66,67,68]). The general dominance of the effects of the 3–4 km deep sill-shaped deformation source makes it challenging to evaluate the contribution of other potential sources of deformation, particularly those at greater depths. The objective of our work was to enhance the comprehension of the origins of the current unrest and to identify the specific components of the plumbing system that are involved, as well as the extent and manner of their involvement. This will increase the knowledge necessary for forecasting future evolution. The co-occurrence of the transition from subsidence to uplift at Campi Flegrei and the onset of deep contraction below Vesuvio emphasise the significance of studying surface deformation over a broad area that includes both volcanoes. With this objective in mind, we analysed the time series of ground displacements from 1993 to 2010 ERS-ENVISAT SAR images using vbICA. The images have broad spatial coverage, enabling the reliable detection of local deformations. Due to the absence of Neapolitan Volcanoes Continuous GPS (NeVoCGPS) stations at potentially interesting locations until after 2010, GNSS data were not utilised. Furthermore, the oldest stations only began operating in 2000 [69].
Consistent results were obtained using various approaches, including scatter plots, vbICA applied in the Campi Flegrei area, vbICA applied in the Vesuvio area, and vbICA applied in the entire area encompassing both Campi Flegrei and Vesuvio. For the sake of brevity, we will only discuss the ICs obtained with the last approach, when four components are retained.
We ascribe the first two independent components (IC1 and IC2) to sources of volcanic deformation. The remaining two components (IC3 and IC4, Figure S37) are likely non-volcanic in origin. IC3 shows a temporal evolution with a large displacement in the early available dates, possibly picking up data-processing noise. IC4 has no clear trend and a cyclical behaviour, which might be seasonal and indicative of hydrological processes. The spatial patterns are also more scattered with respect to IC1 and IC2, and the overall displacements are quite small (max 3 cm ). For these reasons, we do not consider them further here.
The spatial patterns of IC1 indicate that the central area of Campi Flegrei experiences much more deformation than the rest of the study area. Its temporal function strictly follows the vertical displacement at Rione Terra from levelling and GNSS data (Figure 2).
In the Campi Flegrei area, the spatial patterns of u H and u V are very similar to those of the mean velocities during subsidence from 1995 to 2000 [19] and from 1993 to 2000 [42]. The displacements are consistent with the expansion and contraction of two deformation sources. These sources include a sill-shaped volume at a depth of approximately 3 km and a small volume at a depth of 1–2 km. The former is located below the area of maximum uplift, while the latter is located below Solfatara. They serve to reproduce the large-scale deformation in the Campi Flegrei area and the local deformation in the Solfatara area.
Both sources of deformation align with previous findings when inverting average velocities during the 1995–2000 Campi Flegrei subsidence for the best-fit pressurised triaxial ellipsoid (PTE) located off Rione Terra and the point-pressurised spheroid (PS) located below Solfatara [19]. The sill-shaped deformation source also agrees with the results of other authors [20,21,67], who did not, however, distinguish the anomalous behaviour of the Solfatara area. Although the physical characteristics of the sill-shaped deformation source are unclear, the rock temperature at a depth of 3–4 km below Campi Flegrei allows for the emplacement of a thin magmatic sill [70,71].
The IC1 temporal function fully represents the activity of the sill-shaped source, which is responsible for large-scale deformation at Campi Flegrei. The potency histories of the PTE and PS, obtained by inverting the same SAR displacement time series as in the present work epoch by epoch, are similar but not identical, as reported by [19]. Furthermore, the deformation source below Solfatara has been contracting since 2015, while the sill-shaped source has been expanding [42]. Local anomalies are also present in the IC2, IC3, and IC4 spatial distributions at Solfatara, each with distinct characteristics and different temporal evolution. Therefore, we deem it implausible to attribute the deformation anomaly in the Solfatara area to a distortion of the strain field caused by the presence of stationary local rock heterogeneities.
From this point forward, we will refer to the sill-shaped source and the deformation source located below Solfatara, which were obtained in this study, as S1 and S2, respectively. This will avoid confusion between the deformation sources related to the first independent component and those related to the second one.
This study’s key findings concern the identification and analysis of IC2, which represents processes occurring at both Campi Flegrei and Vesuvio (Figure 11). The overall patterns of u H and u V are consistent with subsidence in the Vesuvio area and uplift in the eastern part of Campi Flegrei. In the Campi Flegrei area, IC2 characterises most of the deformation that cannot be attributed to the S1 and S2 sources. The related ground displacements are small, resulting in somewhat noisy spatial distributions. However, maximum values of u V are comparable to those related to IC1 during 2002–2010, at approximately 4 cm. It should be noted that the spatial pattern of u V closely resembles Figure 14d of [19], which displays the vertical displacements of March 2010 relative to February 1993, with the displacements due to the PTE and PS sources (i.e., the sources analogous to S1 and S2) removed.
The u V pattern of IC2 (Figure 11b) is characterised by an area with very small u V in the central part of the caldera. This area borders the area of maximum uplift for IC2 and is where the spatial pattern of IC1 shows its highest values. This feature of IC2 could be due to possible leakage or cross-talk between the independent components. If that were the case, and the vertical spatial pattern associated with IC2 was positive in the whole central-eastern part of Campi Flegrei, we might hypothesise the presence of a single expansive source below eastern Campi Flegrei. Otherwise, it is necessary to consider the effects of two deformation sources with opposite polarities. In both cases, to fully explain the IC2 spatial pattern in the whole study area, which encompasses Campi Flegrei and Vesuvio, it is necessary to assume the existence of an additional contracting source located beneath Vesuvio. Based on a simple mathematical model, which includes two pressurised point spheroids beneath Campi Flegrei and one beneath Vesuvio, all embedded in a homogeneous half-space, IC2 could be the result of two deformation sources located deeper than roughly 8 km below Campi Flegrei, as well as contraction at a depth of approximately 9 km beneath Vesuvio. Therefore, IC2 may reflect the processes taking place in the deep part of the plumbing systems of Campi Flegrei and Vesuvio.
The two deep sources of deformation below Campi Flegrei are located in different areas of the caldera (Figure 15c,d). One is situated in the central area, while the other is situated in the southwestern part. The planimetric location of the latter source is not well constrained due to its non-barycentric position with respect to the available data. However, it appears to be situated between the rims of the caldera’s collapses due to the Campanian Ignimbrite (about 39 ka) and the Neapolitan Yellow Tuff (about 15 ka) eruptions (Figure 1). The deformation source located in the southwestern part is deeper than the central one and possesses similar but opposing potencies. The two sources may be interpreted as part of a deep reservoir from which magma and/or magmatic fluids may transfer.
The histories of the processes occurring in the Campi Flegrei and Vesuvio areas are identified by comparing the temporal functions of the first and second independent components. IC1 is associated with processes that occur in the first few kilometres below Campi Flegrei, while IC2 is linked to deeper processes that involve both Campi Flegrei and Vesuvio. On average, the deep sources seem to have been inactive, or at least poorly active, during the Campi Flegrei subsidence that occurred before the onset of the 2000 mini-uplift (see Figure 2 and Figure 9). The activity began at least two years before the current uplift of Campi Flegrei, during or after the data break that lasted from December 2000 to October 2002. After the data break, the activity of the deep sources steadily increased. As for Campi Flegrei, the source located below the central part of the caldera mostly expanded, while the southwestern source mostly contracted. The existence of the two sources may be related to the likely presence of multiple melt pockets in solid rocks at depths greater than approximately 8 km [36].
The coexistence of a shallower sill-like deformation source at a depth of 3–4 km and a deeper deformation source at approximately 8 km depth appears to be a stable feature at Campi Flegrei. This combination explains the deformation pattern observed along the gulf coast between 1400 and 1536, prior to the Monte Nuovo eruption of 1538. In addition, during the same period, a significant uplift occurred in the southwestern part of the caldera [70]. The same combination explains the post-1538 subsidence phase [3].
During the two mini-uplifts (around 2000 and 2005), the IC1 temporal evolution clearly shows a reversal in its motion. IC2 shows excursions starting/finishing slightly before the beginning/ending of the two mini-uplifts. This is more visible in the temporal function of IC2 when applying vbICA to the Campi Flegrei area only; see Figure 6d. Given the temporal resolution, it is difficult to assess if these excursions are due to two different sources related to the two ICs or to a single moving source. In fact, the adopted ICA algorithm relies on the assumption of non-moving sources, and it could split the deformation signal into two ICs if this assumption was not satisfied (e.g., [72]). In any case, both mini-uplifts can be associated with a single expanding source (S1). This association was previously suggested for the 2005 mini-uplift, based on levelling and cGPS data [73]. If the excursions in IC2 relate to a separate source and they are not an artifact induced by a single moving source, we can interpret them as corresponding to a contraction of the central deep source and an expansion of the southwestern deep source. This seems plausible given that the time lags are in line with the delay between the activation of the deep deformation sources and the start of the ongoing Campi Flegrei uplift. Due to the data break between 2000 and 2002 and the large fluctuations of the IC2 temporal function, making a definitive statement is unfortunately not possible.
The correlation between the transition from subsidence to uplift at Campi Flegrei and the onset of contraction below Vesuvio, i.e., the negative correlation between Campi Flegrei and Vesuvio deformation, is consistent with historical data, suggesting that periods of subsidence at Campi Flegrei are associated with eruptions at Vesuvio and vice versa [74]. This correlation suggests a potential transfer of magma and/or magmatic fluids between the two plumbing systems at a depth of approximately 8–9 km. This possible link between Campi Flegrei and Vesuvio is supported by various geophysical, geochemical, and petrological studies. Seismic data suggest the existence of a partially molten layer beneath Vesuvio and Campi Flegrei at approximately 8 km depth [23,31]. Gravity data suggest the presence of a low-density volume beneath the entire Campanian active volcanic area at depths exceeding 8–9 km [36]. Petrological data suggest that a magma layer extends beneath the entire Neapolitan volcanic area at a depth greater than 8 km (e.g., [35]).
Figure 16a presents a schematic three-dimensional view of the principal processes that may have occurred in the Campi Flegrei and Vesuvio area during 2002–2010, as derived from our analysis of ground displacement data. It is based on a simplified approximation of reality, yet it is able to explain the main spatial and temporal characteristics of the surface deformation observed in the entire Campi Flegrei and Vesuvio region. In the near future, we will attempt to improve our analysis by taking account of known three-dimensional rock heterogeneity, topography, and possible hydrologically induced deformation.
The highlighted interaction between Campi Flegrei and Vesuvio is consistent with repeated evidence that adjacent volcanoes can interact through stress changes caused by variations in pressure, permeability, and magmatic connections. These connections may be due to dykes, a common magmatic system at depth, or direct fluid transfer [75]. Studies have investigated the short-term interactions between Mauna Loa and Kilauea (Hawaii) in 2002 [76,77], Aira and Kirishima (Japan) between 2009 and 2013 [78], and the existing connections between six volcanoes at the Galápagos Islands [79].

5. Conclusions

To better constrain the hazard in the Neapolitan region, it is important to understand the dynamics of the ongoing unrest at Campi Flegrei. In particular, determining whether the unrest is caused by a shallow intrusion of magma or solely due to the influx of magmatic fluids into the hydrothermal system can give rise to different alert scenarios. Ground deformation data can provide valuable insights. However, the effects of a 3–4 km deep sill-shaped deformation source off Rione Terra (S1 in this study) usually dominate, making it challenging to evaluate the contribution of other potential sources of deformation, particularly those at greater depths.
Between 2000 and 2005, Campi Flegrei underwent a shift from subsidence to uplift. The study of the initial stages of subsidence or uplift, or the transition between the two, can better highlight these weak, deeper processes. The transition from subsidence to uplift at Campi Flegrei coincided with the onset of a contraction of a deep volume below Vesuvio. Thus, we analysed the ground displacement time series—obtained from ERS/ENVISAT SAR imagery—for the entire area that includes both Campi Flegrei and Vesuvio from 1993 to 2010. Consistent results were obtained using both simple scatter plots and variational Bayesian independent component analysis (vbICA).
The main findings relate to Campi Flegrei. The mini-uplifts of 2000 and 2005 can be explained as the result of the expansion of the sole S1, whose contraction had previously caused the subsidence of the caldera. Following a period of missing data, including over 2001 and 2002, S1 contracted slowly until 2005, when it began to expand. The contraction and expansion of S1 that occurred after 2002 coincided with an uplift at a nearly constant rate. This primarily affected the eastern area of the caldera and generated vertical displacements that were comparable to those in the central part. Although the latter deformation pattern is complex and difficult to interpret, it is likely related to deformation sources located at an 8–9 km depth.
Considering both Campi Flegrei and Vesuvio, the correlation observed after 2002 between uplift at Campi Flegrei and subsidence at Vesuvio suggests the possibility of the transfer of magma and/or magmatic fluids between the two plumbing systems at a depth of 8–9 km.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/rs16101717/s1, SupplementaryMaterial_uH_uV.pdf including Figure S1: Combination of line-of-sight displacements to obtain horizontal and quasi-vertical displacements, SupplementaryMaterial_PressurizedEllipsoid.pdf including Figure S2: Approximation of a pressurised ellipsoidal cavity through seven point sources, SupplementaryMaterial_ScatterPlots.pdf including Figures S3–S36: Scatter plots of quasi-vertical displacement u V and horizontal displacement u H at various locations in Campi Flegrei compared to u V in the area of maximum pre-2000 subsidence from 1993 to 2010, SupplementaryMaterial_IC3_IC4.pdf including Figure S37: Temporal functions and spatial distributions of IC3 and IC4, SupplementaryMaterial_FEMvalidation.pdf including Figure S38: Comparison between the seven-point-sources approximated approach and the finite element method (FEM) for calculating the radial and vertical superficial displacements caused by a finite pressurised spheroid with a vertical polar axis. Reference [45] is cited in the Supplementary Materials.

Author Contributions

Conceptualisation, A.A. and L.C.; methodology, A.A. and L.C.; software, A.A., A.G. and L.C.; formal analysis, A.A. and L.C.; investigation, A.A.; writing—original draft preparation, A.A.; writing—review and editing, A.A., A.G. and L.C.; visualisation, A.A. and L.C.; supervision, A.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the University of Salerno, grant number 300391FRB20AMORU.

Data Availability Statement

Figures have been produced using open-source software GMT v. 6.3.0 (https://www.generic-mapping-tools.org/, accessed on 21 March 2024), Veusz v. 3.6.2 (https://veusz.github.io/, accessed on 21 March 2024), Paraview v. 5.10.1 (https://paraview.org, accessed on 21 March 2024), LibreOffice v. 24.2.2.2 (https://libreoffice.org, accessed on 21 March 2024), and Inkscape v. 1.1.2 (https://inkscape.org/, accessed on 21 March 2024). ERS-ENVISAT ground displacement time series are taken from [19].

Acknowledgments

We would like to express our gratitude to R. Lanari and M. Manzo from IREA-CNR (Napoli, Italy) for providing us with the DInSAR LOS ground displacement time series. Valerio Acocella is acknowledged for his valuable contributions. He collaborated with A.A. and L.C. at a preliminary stage of this study. Additionally, we extend our appreciation to the authors and maintainers of the open-source software used in this work.

Conflicts of Interest

The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:
SARsynthetic aperture radar
SBASsmall baseline subset
DInSAR       differential synthetic aperture radar interferometry
LOSline of sight
DEMdigital elevation model
GNSSglobal navigation satellite system
NeVoCGPSNeapolitan Volcanoes Continuous GPS network
vbICAvariational Bayesian independent component analysis
ICindependent component
ARDautomatic relevance determination
IAincidence angle
ESAEuropean Space Agency
SNAPSentinel Application Platform

References

  1. Dzurisin, D. Volcano Deformation, Geodetic Monitoring Techniques; Springer: Berlin, Germany, 2007. [Google Scholar]
  2. Guidoboni, E.; Ciuccarelli, C. The Campi Flegrei caldera: Historical revision and new data on seismic crises, bradyseisms, the Monte Nuovo eruption and ensuing earthquakes (twelfth century 1582 AD). Bull. Volcanol. 2011, 73, 655–677. [Google Scholar] [CrossRef]
  3. Trasatti, E.; Magri, C.; Acocella, V.; Del Gaudio, C.; Ricco, C.; Di Vito, M.A. Magma transfer at Campi Flegrei caldera (Italy) after the 1538 AD eruption. Geophys. Res. Lett. 2023, 50, e2022GL102437. [Google Scholar] [CrossRef]
  4. 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. 2010, 185, 48–56. [Google Scholar] [CrossRef]
  5. Rosi, M.; Acocella, V.; Cioni, R.; Bianco, F.; Costa, A.; De Martino, P.; Giordano, G.; Inguaggiato, S. Defining the Pre-Eruptive States of Active Volcanoes for Improving Eruption Forecasting. Front. Earth Sci. 2022, 10, 795700. [Google Scholar] [CrossRef]
  6. Bollettino Mensile Campi Flegrei, February 2024. Available online: https://www.ov.ingv.it/index.php/monitoraggio-e-infrastrutture/bollettini-tutti/bollett-mensili-cf/anno-2024-3/1586-bollettino-mensile-campi-flegrei-2024-02/file (accessed on 18 March 2024).
  7. Tarquini, S.; Isola, I.; Favalli, M.; Mazzarini, F.; Bisson, M.; Pareschi, M.T.; Boschi, E. TINITALY/01: A new triangular irregular network of Italy. Ann. Geophys. 2007, 50, 407–425. [Google Scholar] [CrossRef]
  8. Vitale, S.; Isaia, R. Fractures and faults in volcanic rocks (Campi Flegrei, southern Italy): Insight into volcano-tectonic processes. Int. J. Earth Sci. 2014, 103, 801–819. [Google Scholar] [CrossRef]
  9. Giudicepietro, F.; Ricciolino, P.; Bianco, F.; Caliro, S.; Cubellis, E.; D’Auria, L.; De Cesare, W.; De Martino, P.; Esposito, A.M.; Galluzzo, D.; et al. Campi Flegrei, Vesuvius and Ischia Seismicity in the Context of the Neapolitan Volcanic Area. Front. Earth Sci. 2021, 9, 662113. [Google Scholar] [CrossRef]
  10. Tramelli, A.; Giudicepietro, F.; Ricciolino, P.; Chiodini, G. The seismicity of Campi Flegrei in the contest of an evolving long term unrest. Sci. Rep. 2022, 12, 2900. [Google Scholar] [CrossRef]
  11. Danesi, S.; Pino, N.A.; Carlino, S.; Kilburn, C.R.J. Evolution in unrest processes at Campi Flegrei caldera as inferred from local seismicity. Earth Planet. Sci. Lett. 2024, 626, 118530. [Google Scholar] [CrossRef]
  12. Chiodini, G.; Paonita, A.; Aiuppa, A.; Costa, A.; Caliro, S.; De Martino, P.; Acocella, V.; Vandemeulebrouck, J. Magmas near the critical degassing pressure drive volcanic unrest towards a critical state. Nat. Commun. 2016, 7, 13712. [Google Scholar] [CrossRef]
  13. De Siena, L.; Chiodini, G.; Vilardo, G.; Del Pezzo, E.; Castellano, M.; Colombelli, S.; Tisato, N.; Ventura, G. Source and dynamics of a volcanic caldera unrest: Campi Flegrei, 1983–84. Sci. Rep. 2017, 7, 8099. [Google Scholar] [CrossRef] [PubMed]
  14. Moretti, R.; Arienzo, I.; Civetta, L.; Orsi, G.; Papale, P. Multiple magma degassing sources at an explosive volcano. Earth Planet. Sci. Lett. 2013, 367, 95–104. [Google Scholar] [CrossRef]
  15. Chiodini, G.; Caliro, S.; Avino, R.; Bini, G.; Giudicepietro, F.; De Cesare, W.; Ricciolino, P.; Aiuppa, A.; Cardellini, C.; Petrillo, Z.; et al. Hydrothermal pressure-temperature control on CO2 emissions and seismicity at Campi Flegrei (Italy). J. Volcanol. Geotherm. Res. 2021, 414, 107245. [Google Scholar] [CrossRef]
  16. 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. 2015, 414, 58–67. [Google Scholar] [CrossRef]
  17. Buono, G.; Paonita, A.; Pappalardo, L.; Caliro, S.; Tramelli, A.; Chiodini, G. New insights into the recent magma dynamics under Campi Flegrei caldera (Italy) from petrological and geochemical evidence. J. Geophys. Res. Solid Earth 2022, 127, e2021JB023773. [Google Scholar] [CrossRef]
  18. De Martino, P.; Dolce, M.; Brandi, G.; Scarpato, G. Campi Flegrei cGPS Weekly Positions Time Series [Dataset]; Zenodo: Geneve, Switzerland, 2023. [Google Scholar] [CrossRef]
  19. Amoruso, A.; Crescentini, L.; Sabbetta, I. Paired deformation sources of the Campi Flegrei caldera (Italy) required by recent (1980–2010) deformation history. J. Geophys. Res. Solid Earth 2014, 119, 858–879. [Google Scholar] [CrossRef]
  20. D’Auria, L.; Pepe, S.; Castaldo, R.; Giudicepietro, F.; Macedonio, G.; Ricciolino, P.; Tizzani, P.; Casu, F.; Lanari, R.; Manzo, M.; et al. Magma injection beneath the urban area of Naples: A new mechanism for the 2012–2013 volcanic unrest at Campi Flegrei caldera. Sci. Rep. 2015, 5, 13100. [Google Scholar] [CrossRef] [PubMed]
  21. Castaldo, R.; Tizzani, P.; Solaro, G. Inflating Source Imaging and Stress/Strain Field Analysis at at Campi Flegrei Caldera: The 2009–2013 Unrest Episode. Remote Sens. 2021, 13, 2298. [Google Scholar] [CrossRef]
  22. De Siena, L.; Amoruso, A.; Pezzo, E.D.; Wakeford, Z.; Castellano, M.; Crescentini, L. Space-weighted seismic attenuation mapping of the aseismic source of Campi Flegrei 1983–1984 unrest. Geophys. Res. Lett. 2017, 44, 1740–1748. [Google Scholar] [CrossRef]
  23. Zollo, A.; Maercklin, N.; Vassallo, M.; Dello Iacono, D.; Virieux, J.; Gasparini, P. Seismic reflections reveal a massive melt layer feeding Campi Flegrei caldera. Geophys. Res. Lett. 2008, 35, L12306. [Google Scholar] [CrossRef]
  24. Fournier, N.; Chardot, L. Understanding volcano hydrothermal unrest from geodetic observations: Insights from numerical modeling and application to White Island volcano, New Zealand. J. Geophys. Res. 2012, 117, B11208. [Google Scholar] [CrossRef]
  25. Amoruso, A.; Crescentini, L.; Scarpa, R.; Bilham, R.; Linde, A.T.; Sacks, I.S. Abrupt magma chamber contraction and microseismicity at Campi Flegrei, Italy: Cause and effect determined from strainmeters and tiltmeters. J. Geophys. Res. Solid Earth 2015, 120, 5467–5478. [Google Scholar] [CrossRef]
  26. Troise, C.; De Natale, G.; Schiavone, R.; Somma, R.; Moretti, R. The Campi Flegrei caldera unrest: Discriminating magma intrusions from hydrothermal effects and implications for possible evolution. Earth Sci. Rev. 2019, 188, 108–122. [Google Scholar] [CrossRef]
  27. Cioni, R.; Santacroce, R.; Sbrana, A. Pyroclastic deposits as a guide for reconstructing the multi-stage evolution of the Somma-Vesuvius caldera. Bull. Volcanol. 1999, 60, 207–222. [Google Scholar] [CrossRef]
  28. Scaillet, B.; Cioni, R.; Pichavant, M. Upward migration of Vesuvius magma chamber over the past 20,000 years. Nature 2008, 455, 216–219. [Google Scholar] [CrossRef]
  29. Balcone-Boissard, H.; Boudon, G.; Cioni, R.; Webster, J.D.; Zdanowicz, G.; Orsi, G.; Civetta, L. Chlorine as a geobarometer for alkaline magmas: Evidence from a systematic study of the eruptions of Mount Somma-Vesuvius. Sci. Rep. 2016, 6, 21726. [Google Scholar] [CrossRef] [PubMed]
  30. Di Stefano, R.; Chiarabba, C. Active source tomography at Mt. Vesuvius: Constraints for the magmatic system. J. Geophys. Res. 2002, 107, 2278. [Google Scholar] [CrossRef]
  31. Auger, E.; Gasparini, P.; Virieux, J.; Zollo, A. Seismic Evidence of an Extended Magmatic Sill under Mt. Vesuvius. Science 2001, 294, 1510–1512. [Google Scholar] [CrossRef]
  32. Nunziata, C.; Natale, M.; Luongo, G.; Panza, G.F. Magma reservoir at Mt. Vesuvius: Size of the hot, partially molten, crust material detected deeper than 8 km. Earth Planet. Sci. Lett. 2006, 242, 51–57. [Google Scholar] [CrossRef]
  33. Amoruso, A.; Crescentini, L. DInSAR Data Reveal an Intriguing Contemporaneous Onset of Deep Deflation below Vesuvio and the Ongoing Campi Flegrei Uplift. Remote Sens. 2023, 15, 3038. [Google Scholar] [CrossRef]
  34. D’Auria, L.; Esposito, A.M.; Lo Bascio, D.; Ricciolino, P.; Giudicepietro, F.; Martini, M.; Caputo, T.; De Cesare, W.; Orazi, M.; Peluso, R.; et al. The recent seismicity of Mt. Vesuvius: Inference on seismogenic processes. Ann. Geophys. 2013, 56, S0442. [Google Scholar] [CrossRef]
  35. Pappalardo, L.; Buono, G. Insights into Processes and Timescales of Magma Storage and Ascent From Textural and Geochemical Investigations: Case Studies from High-Risk Neapolitan Volcanoes (Italy). In Crustal Magmatic System Evolution: Anatomy, Architecture, and Physico-Chemical Processes; Masotta, M., Beier, C., Mollo, S., Eds.; John Wiley and Sons: New York, NY, USA, 2021; Chapter 10; pp. 213–235. [Google Scholar] [CrossRef]
  36. Fedi, M.; Cella, F.; D’Antonio, M.; Florio, G.; Paoletti, V.; Morra, V. Gravity modeling finds a large magma body in the deep crust below the Gulf of Naples, Italy. Sci. Rep. 2018, 8, 8229. [Google Scholar] [CrossRef] [PubMed]
  37. Crescentini, L.; Amoruso, A. Effects of crustal layering on the inversion of deformation and gravity data in volcanic areas: An application to the Campi Flegrei caldera, Italy. Geophys. Res. Lett. 2007, 34, L09303. [Google Scholar] [CrossRef]
  38. Kilburn, C.R.; De Natale, G.; Carlino, S. Progressive approach to eruption at Campi Flegrei caldera in southern Italy. Nat. Commun. 2017, 8, 15312. [Google Scholar] [CrossRef] [PubMed]
  39. Ers Mission Description. Available online: https://earth.esa.int/eogateway/missions/ers (accessed on 21 March 2024).
  40. Envisat Mission Description. Available online: https://earth.esa.int/eogateway/missions/envisat (accessed on 21 March 2024).
  41. 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. 2002, 40, 2375–2383. [Google Scholar] [CrossRef]
  42. Amoruso, A.; Crescentini, L. Clues of Ongoing Deep Magma Inflation at Campi Flegrei Caldera (Italy) from Empirical Orthogonal Function Analysis of SAR Data. Remote Sens. 2022, 14, 5698. [Google Scholar] [CrossRef]
  43. De Matteo, A.; Massa, B.; Castaldo, R.; D’Auria, L.; James, M.R.; Lane, S.J.; Pepe, S.; Tizzani, P. An integrated modeling approach for analyzing the deformation style of active volcanoes: Somma-Vesuvius case study. J. Geophys. Res. Solid Earth 2022, 127, e2021JB022338. [Google Scholar] [CrossRef]
  44. Casu, F.; Manzo, M.; Lanari, R. A quantitative assessment of the SBAS algorithm performance for surface deformation retrieval from DInSAR data. Remote Sens. Environ. 2006, 102, 195–210. [Google Scholar] [CrossRef]
  45. Fuhrmann, T.; Garthwaite, M.C. Resolving Three-Dimensional Surface Motion with InSAR: Constraints from Multi-Geometry Data Fusion. Remote Sens. 2019, 11, 241. [Google Scholar] [CrossRef]
  46. Comon, P. Independent component analysis, A new concept? Signal Process. 1994, 36, 287–314. [Google Scholar] [CrossRef]
  47. Hyvärinen, A.; Oja, E. A fast fixed-point algorithm for independent component analysis. Neural Comput. 1997, 9, 1483–1492. [Google Scholar] [CrossRef]
  48. Choudrey, R.; Roberts, S. Variational mixture of Bayesian independent component analyzers. Neural Comput. 2003, 15, 213–252. [Google Scholar] [CrossRef] [PubMed]
  49. Gualandi, A.; Serpelloni, E.; Belardinelli, M. Blind source separation problem in GPS time series. J. Geod. 2016, 90, 323–341. [Google Scholar] [CrossRef]
  50. Gualandi, A.; Liu, Z. Variational Bayesian independent component analysis for InSAR displacement time-series with application to central California, USA. J. Geophys. Res. Solid Earth 2021, 126, e2020JB020845. [Google Scholar] [CrossRef]
  51. Jolliffe, I.T. Principal Component Analysis, 2nd ed.; Springer: New York, NY, USA, 2002. [Google Scholar] [CrossRef]
  52. MacKay, D.J.C. Bayesian non-linear modelling for the prediction competition. ASHRAE Trans. 1994, 100, 1053–1062. [Google Scholar]
  53. Amoruso, A.; Crescentini, L.; Fidani, C. Effects of crustal layering on source parameter inversion from coseismic geodetic data. Geophys. J. Int. 2004, 159, 353–364. [Google Scholar] [CrossRef]
  54. Amoruso, A.; Crescentini, L.; Berrino, G. Simultaneous inversion of deformation and gravity changes in a horizontally layered half-space: Evidences for magma intrusion during the 1982–1984 unrest at Campi Flegrei caldera (Italy). Earth Planet. Sci. Lett. 2008, 272, 181–188. [Google Scholar] [CrossRef]
  55. Amoruso, A.; Crescentini, L. Modelling deformation due to a pressurized ellipsoidal cavity, with reference to the Campi Flegrei caldera, Italy. Geophys. Res. Lett. 2011, 38, L01303. [Google Scholar] [CrossRef]
  56. Davis, P.M. Surface deformation due to inflation of an arbitrarily oriented triaxial ellipsoidal cavity in an elastic half-space, with reference to Kilauea Volcano, Hawaii. J. Geophys. Res. 1986, 91, 7429–7438. [Google Scholar] [CrossRef]
  57. Eshelby, J.D. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc. R. Soc. Lond. Ser. A 1957, 241, 376–396. [Google Scholar] [CrossRef]
  58. Mindlin, R.D. Force at a point in the interior of a semi-infinite solid. Physics 1936, 7, 195–202. [Google Scholar] [CrossRef]
  59. Wang, R.; Lorenzo Martín, F.; Roth, F. PSGRN/PSCMP–A new code for calculating co- and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory. Comput. Geosci. 2006, 32, 527–541. [Google Scholar] [CrossRef]
  60. Judenherc, S.; Zollo, A. The Bay of Naples (Southern Italy): Constraints on the volcanic structures inferred from a dense seismic survey. J. Geophys. Res. 2004, 109, B10312. [Google Scholar] [CrossRef]
  61. Ludwig, W.J.; Nafe, J.E.; Drake, C.L. Seismic refraction. In The Sea; Maxwell, A.E., Ed.; Wiley-Interscience: New York, NY, USA, 1970; Volume 4, pp. 53–84. [Google Scholar]
  62. Ingber, L. Simulated annealing: Practice versus theory. Math. Comput. Model. 1993, 18, 29–57. [Google Scholar] [CrossRef]
  63. Sambridge, M. Geophysical inversion with a neighbourhood algorithm—I. Searching a parameter space. Geophys. J. Int. 1999, 138, 479–494. [Google Scholar] [CrossRef]
  64. Sambridge, M. Geophysical inversion with a neighbourhood algorithm—II. Appraising the ensemble. Geophys. J. Int. 1999, 138, 727–746. [Google Scholar] [CrossRef]
  65. Aki, K.; Richards, P.G. Quantitative Seismology, 2nd ed.; University Science Books: Sausalito, CA, USA, 2002. [Google Scholar]
  66. Petrelli, M.; López, M.A.; Pisello, A.; Perugini, D. Pre-eruptive dynamics at the Campi Flegrei Caldera: From evidence of magma mixing to timescales estimates. Earth Planets Space 2023, 75, 19. [Google Scholar] [CrossRef]
  67. Kilburn, C.R.J.; Carlino, S.; Danesi, S.; Pino, N.A. Potential for rupture before eruption at Campi Flegrei caldera, Southern Italy. Commun. Earth Environ. 2023, 4, 190. [Google Scholar] [CrossRef]
  68. De Siena, L.; Amoruso, A.; Petrosino, S.; Crescentini, L. Geophysical responses to an environmentally-boosted volcanic unrest. Geophys. Res. Lett. 2024, 51, e2023GL104895. [Google Scholar] [CrossRef]
  69. De Martino, P.; Dolce, M.; Brandi, G.; Scarpato, G.; Tammaro, U. The Ground Deformation History of the Neapolitan Volcanic Area (Campi Flegrei Caldera, Somma–Vesuvius Volcano, and Ischia Island) from 20 Years of Continuous GPS Observations (2000–2019). Remote Sens. 2021, 13, 2725. [Google Scholar] [CrossRef]
  70. Amoruso, A.; Crescentini, L.; D’Antonio, M.; Acocella, V. Thermally-assisted magma emplacement explains restless calderas. Sci. Rep. 2017, 7, 7948. [Google Scholar] [CrossRef]
  71. Amoruso, A.; Crescentini, L. An approximate approach to nonisothermal emplacement of kilometer-sized kilometer-deep sills at calderas. J. Geophys. Res. Solid Earth 2019, 124, 1236–1253. [Google Scholar] [CrossRef]
  72. Larochelle, S.; Gualandi, A.; Chanard, K.; Avouac, J.-P. Identification and Extraction of Seasonal Geodetic Signals Due to Surface Load Variations. J. Geophys. Res. Solid Earth 2018, 123, 11031–11047. [Google Scholar] [CrossRef]
  73. Amoruso, A.; Crescentini, L.; Linde, A.T.; Sacks, I.S.; Scarpa, R.; Romano, P. A horizontal crack in a layered structure satisfies deformation for the 2004–2006 uplift of Campi Flegrei. Geophys. Res. Lett. 2007, 34, L22313. [Google Scholar] [CrossRef]
  74. Walter, T.R.; Shirzaei, M.; Manconi, A.; Solaro, G.; Pepe, A.; Manzo, M.; Sansosti, E. Possible coupling of Campi Flegrei and Vesuvius as revealed by InSAR time series, correlation analysis and time dependent modeling. J. Volcanol. Geotherm. Res. 2014, 280, 104–110. [Google Scholar] [CrossRef]
  75. Biggs, J.; Robertson, E.; Cashman, K. The lateral extent of volcanic interactions during unrest and eruption. Nat. Geosci. 2016, 9, 308–311. [Google Scholar] [CrossRef]
  76. Gonnermann, H.M.; Foster, J.H.; Pol, M.; Wolfe, C.J.; Brooks, B.A.; Miklius, A. Coupling at Mauna Loa and Kilauea by stress transfer in an asthenospheric melt layer. Nat. Geosci. 2012, 5, 826–829. [Google Scholar] [CrossRef]
  77. Miklius, A.; Cervelli, P. Interaction between Kilauea and Mauna Loa. Nature 2003, 421, 229. [Google Scholar] [CrossRef]
  78. Brothelande, E.; Amelung, F.; Yunjun, Z.; Wdowinski, S. Geodetic evidence for interconnectivity between Aira and Kirishima magmatic systems, Japan. Sci. Rep. 2018, 8, 9811. [Google Scholar] [CrossRef]
  79. Reddin, E.; Ebmeier, S.K.; Rivalta, E.; Bagnardi, M.; Baker, S.; Bell, A.F.; Mothes, P.; Aguaiza, S. Magmatic connectivity among six Galápagos volcanoes revealed by satellite geodesy. Nat. Commun. 2023, 14, 6614. [Google Scholar] [CrossRef]
Figure 3. Displacement data dates (reprinted from [33]). The acquisition dates of the ascending- and descending-orbit images are represented by red and blue circles, respectively. The green circles indicate the dates of the computed horizontal ( u H ) and quasi-vertical ( u V ) displacements.
Figure 3. Displacement data dates (reprinted from [33]). The acquisition dates of the ascending- and descending-orbit images are represented by red and blue circles, respectively. The green circles indicate the dates of the computed horizontal ( u H ) and quasi-vertical ( u V ) displacements.
Remotesensing 16 01717 g003
Figure 4. Scatter plots of the quasi-vertical displacement u V at selected locations in Campi Flegrei (P09, P06, P15, P17) and in the area of maximum pre-2000 subsidence (P18) from 1993 to 2010. Each location is circular, with a diameter of 600 m. The map in (a) displays all the sampled locations. The map uses a colour scheme to represent locations characterised by similar scatter plots. Cyan indicates locations with noisy almost-zero displacements during Campi Flegrei subsidence and uplift after then; an example is shown in (c). Blue indicates locations with a positive correlation with P18 during Campi Flegrei subsidence and uplift after then; an example is shown in (d). Red indicates locations with a positive correlation with P18 during the whole 1993–2010 period; an example is shown in (e). Magenta indicates locations with a positive correlation with P18 during the whole 1993–2010 period, with some peculiarities; an example is shown in (f). The epoch is identified by the colour of the symbols in each scatter plot, which can be seen in (b). The scatter plots for all the sampled locations are shown in the Supplementary Information Figures S3–S19.
Figure 4. Scatter plots of the quasi-vertical displacement u V at selected locations in Campi Flegrei (P09, P06, P15, P17) and in the area of maximum pre-2000 subsidence (P18) from 1993 to 2010. Each location is circular, with a diameter of 600 m. The map in (a) displays all the sampled locations. The map uses a colour scheme to represent locations characterised by similar scatter plots. Cyan indicates locations with noisy almost-zero displacements during Campi Flegrei subsidence and uplift after then; an example is shown in (c). Blue indicates locations with a positive correlation with P18 during Campi Flegrei subsidence and uplift after then; an example is shown in (d). Red indicates locations with a positive correlation with P18 during the whole 1993–2010 period; an example is shown in (e). Magenta indicates locations with a positive correlation with P18 during the whole 1993–2010 period, with some peculiarities; an example is shown in (f). The epoch is identified by the colour of the symbols in each scatter plot, which can be seen in (b). The scatter plots for all the sampled locations are shown in the Supplementary Information Figures S3–S19.
Remotesensing 16 01717 g004
Figure 5. Scatter plots of the horizontal displacement u H at selected locations in Campi Flegrei (P03, P06, P07, P11) and the quasi-vertical displacement u V in the area of maximum pre-2000 subsidence (P18) from 1993 to 2010. Each location is circular, with a diameter of 600 m. The map in (a) displays all the sampled locations. The map uses a colour scheme to represent locations characterised by similar scatter plots. Cyan indicates locations with noisy almost-zero horizontal displacements during Campi Flegrei subsidence and westward horizontal displacements after then; an example is shown in (c). Blue indicates locations with a positive correlation with the quasi-vertical displacement at P18 during Campi Flegrei subsidence and almost-zero horizontal displacements after then; an example is shown in (d). Red indicates locations with a positive correlation with the quasi-vertical displacement at P18 during Campi Flegrei subsidence and westward horizontal displacements after then; the sole example is shown in (e). Magenta indicates locations with a negative correlation with the quasi-vertical displacement at P18 during Campi Flegrei subsidence and westward horizontal displacements after then; the sole example is shown in (f). The epoch is identified by the colour of the symbols in each scatter plot, which can be seen in (b). The scatter plots for all the sampled locations are shown in the Supplementary Information Figures S20–S36.
Figure 5. Scatter plots of the horizontal displacement u H at selected locations in Campi Flegrei (P03, P06, P07, P11) and the quasi-vertical displacement u V in the area of maximum pre-2000 subsidence (P18) from 1993 to 2010. Each location is circular, with a diameter of 600 m. The map in (a) displays all the sampled locations. The map uses a colour scheme to represent locations characterised by similar scatter plots. Cyan indicates locations with noisy almost-zero horizontal displacements during Campi Flegrei subsidence and westward horizontal displacements after then; an example is shown in (c). Blue indicates locations with a positive correlation with the quasi-vertical displacement at P18 during Campi Flegrei subsidence and almost-zero horizontal displacements after then; an example is shown in (d). Red indicates locations with a positive correlation with the quasi-vertical displacement at P18 during Campi Flegrei subsidence and westward horizontal displacements after then; the sole example is shown in (e). Magenta indicates locations with a negative correlation with the quasi-vertical displacement at P18 during Campi Flegrei subsidence and westward horizontal displacements after then; the sole example is shown in (f). The epoch is identified by the colour of the symbols in each scatter plot, which can be seen in (b). The scatter plots for all the sampled locations are shown in the Supplementary Information Figures S20–S36.
Remotesensing 16 01717 g005
Figure 6. Campi Flegrei area: temporal functions and spatial distributions of the first two ICs. Each temporal function starts at zero, and the value corresponding to the maximum of its modulus is non-negative. The full range of variation, i.e., the difference between the maximum and minimum of the temporal function, is 1. The spatial distributions are presented in centimetres of displacement during the period when the related temporal function underwent its maximum variation. IC1 experienced this from the beginning of 1993 to the beginning of 2005, while IC2 experienced it from the end of 1999 to 2010. (a,d) Temporal functions of IC1 and IC2, respectively. (b,c) Spatial distributions of IC1, showing u V and u H , respectively. (e,f) Spatial distributions of IC2, showing u V and u H , respectively. Grey tones in (b,c,e,f) indicate topography (see Figure 1).
Figure 6. Campi Flegrei area: temporal functions and spatial distributions of the first two ICs. Each temporal function starts at zero, and the value corresponding to the maximum of its modulus is non-negative. The full range of variation, i.e., the difference between the maximum and minimum of the temporal function, is 1. The spatial distributions are presented in centimetres of displacement during the period when the related temporal function underwent its maximum variation. IC1 experienced this from the beginning of 1993 to the beginning of 2005, while IC2 experienced it from the end of 1999 to 2010. (a,d) Temporal functions of IC1 and IC2, respectively. (b,c) Spatial distributions of IC1, showing u V and u H , respectively. (e,f) Spatial distributions of IC2, showing u V and u H , respectively. Grey tones in (b,c,e,f) indicate topography (see Figure 1).
Remotesensing 16 01717 g006
Figure 7. Vesuvio area: temporal functions and spatial distributions of the first two ICs. Each temporal function starts at zero, and the value corresponding to the maximum of its modulus is non-negative. The full range of variation, i.e., the difference between the maximum and minimum of the temporal function, is 1. The spatial distributions are presented in centimetres of displacement during the period when the related temporal function underwent its maximum variation. IC1 experienced this from the beginning of 1993 to 2010, while IC2 experienced it from 2001 to 2010. (a,d) Temporal functions of IC1 and IC2, respectively. (b,c) Spatial distributions of IC1, showing u V and u H , respectively. (e,f) Spatial distributions of IC2, showing u V and u H , respectively. Grey tones in (b,c,e,f) indicate topography (see Figure 1).
Figure 7. Vesuvio area: temporal functions and spatial distributions of the first two ICs. Each temporal function starts at zero, and the value corresponding to the maximum of its modulus is non-negative. The full range of variation, i.e., the difference between the maximum and minimum of the temporal function, is 1. The spatial distributions are presented in centimetres of displacement during the period when the related temporal function underwent its maximum variation. IC1 experienced this from the beginning of 1993 to 2010, while IC2 experienced it from 2001 to 2010. (a,d) Temporal functions of IC1 and IC2, respectively. (b,c) Spatial distributions of IC1, showing u V and u H , respectively. (e,f) Spatial distributions of IC2, showing u V and u H , respectively. Grey tones in (b,c,e,f) indicate topography (see Figure 1).
Remotesensing 16 01717 g007
Figure 8. Whole area including Campi Flegrei and Vesuvio: selection of the number of components through automatic relevance determination (ARD). (a) Ratio of minimum-to-maximum posterior mixing matrix variances in relation to the number of components. The ARD criterion suggests retaining either two or four components. The red dashed line indicates the value of 0.04 used by [50] as a threshold to determine the number of components to retain (two in this case). This is an arbitrary threshold, and a different choice (e.g., 0.01 , green dashed line) would lead to the selection of a different number of components (four in this case). The addition of IC3 and IC4 does not affect the interpretation of the first two ICs. (b) Explained variance of the dataset in relation to the number of components.
Figure 8. Whole area including Campi Flegrei and Vesuvio: selection of the number of components through automatic relevance determination (ARD). (a) Ratio of minimum-to-maximum posterior mixing matrix variances in relation to the number of components. The ARD criterion suggests retaining either two or four components. The red dashed line indicates the value of 0.04 used by [50] as a threshold to determine the number of components to retain (two in this case). This is an arbitrary threshold, and a different choice (e.g., 0.01 , green dashed line) would lead to the selection of a different number of components (four in this case). The addition of IC3 and IC4 does not affect the interpretation of the first two ICs. (b) Explained variance of the dataset in relation to the number of components.
Remotesensing 16 01717 g008
Figure 9. Whole area including Campi Flegrei and Vesuvio: temporal functions of the first two ICs. Each temporal function starts at zero, and the value corresponding to the maximum of its modulus is non-negative. The full range of variation, i.e., the difference between the maximum and minimum of the temporal function, is 1. Selecting either 2 or 4 components has a negligible impact on the temporal functions of IC1 and IC2. IC1 accurately tracks the vertical displacement at Rione Terra from levelling data but may appear reversed due to scaling (see also Figure 2). IC2 remains almost constant until around 2001, after which it increases at a steady rate.
Figure 9. Whole area including Campi Flegrei and Vesuvio: temporal functions of the first two ICs. Each temporal function starts at zero, and the value corresponding to the maximum of its modulus is non-negative. The full range of variation, i.e., the difference between the maximum and minimum of the temporal function, is 1. Selecting either 2 or 4 components has a negligible impact on the temporal functions of IC1 and IC2. IC1 accurately tracks the vertical displacement at Rione Terra from levelling data but may appear reversed due to scaling (see also Figure 2). IC2 remains almost constant until around 2001, after which it increases at a steady rate.
Remotesensing 16 01717 g009
Figure 10. Whole area including Campi Flegrei and Vesuvio: spatial distributions of the first independent component. The spatial distributions are given in centimetres of displacement from the beginning of 1993 to the beginning of 2005. This period corresponds to when the related temporal function underwent its maximum variation. The ground displacement history at each pixel can be obtained by multiplying the temporal function by the displacement shown here. Each panel is divided into two parts by a dashed black line. The colour bars to the left and right of the dashed line differ by a factor of 10 to enhance the visibility of the small displacements in the Vesuvio area. (a) Horizontal displacement u H . (b) Quasi-vertical displacement u V .
Figure 10. Whole area including Campi Flegrei and Vesuvio: spatial distributions of the first independent component. The spatial distributions are given in centimetres of displacement from the beginning of 1993 to the beginning of 2005. This period corresponds to when the related temporal function underwent its maximum variation. The ground displacement history at each pixel can be obtained by multiplying the temporal function by the displacement shown here. Each panel is divided into two parts by a dashed black line. The colour bars to the left and right of the dashed line differ by a factor of 10 to enhance the visibility of the small displacements in the Vesuvio area. (a) Horizontal displacement u H . (b) Quasi-vertical displacement u V .
Remotesensing 16 01717 g010
Figure 11. Whole area including Campi Flegrei and Vesuvio: spatial distributions of the second independent component. The spatial distributions are given in centimetres of displacement from the end of 1999 to mid-2010. This period corresponds to when the related temporal function underwent its maximum variation. The ground displacement history at each pixel can be obtained by multiplying the temporal function by the displacement shown here. (a) Horizontal displacement u H . (b) Quasi-vertical displacement u V .
Figure 11. Whole area including Campi Flegrei and Vesuvio: spatial distributions of the second independent component. The spatial distributions are given in centimetres of displacement from the end of 1999 to mid-2010. This period corresponds to when the related temporal function underwent its maximum variation. The ground displacement history at each pixel can be obtained by multiplying the temporal function by the displacement shown here. (a) Horizontal displacement u H . (b) Quasi-vertical displacement u V .
Remotesensing 16 01717 g011
Figure 12. Marginal probability density functions (MPDFs) of the parameters of the two pressurised spheroidal sources that satisfy IC1 at Campi Flegrei. MPDFs were estimated using NAB and scaled to a maximum of 1 to facilitate plot interpretation. The red curves pertain to the source responsible for large-scale deformation (S1), while the blue curves relate to the source responsible for Solfatara local deformation (S2). The dashed lines relate to the parameter values listed in Table 2. (a) Easting position of the source centres. (b) Northing position of the source centres. (c) Depth to the source centres. (d) Volume change of the sources. (e) Polar-to-equatorial radius (i.e., aspect ratio) of the sources. (f) Equatorial radius of S1.
Figure 12. Marginal probability density functions (MPDFs) of the parameters of the two pressurised spheroidal sources that satisfy IC1 at Campi Flegrei. MPDFs were estimated using NAB and scaled to a maximum of 1 to facilitate plot interpretation. The red curves pertain to the source responsible for large-scale deformation (S1), while the blue curves relate to the source responsible for Solfatara local deformation (S2). The dashed lines relate to the parameter values listed in Table 2. (a) Easting position of the source centres. (b) Northing position of the source centres. (c) Depth to the source centres. (d) Volume change of the sources. (e) Polar-to-equatorial radius (i.e., aspect ratio) of the sources. (f) Equatorial radius of S1.
Remotesensing 16 01717 g012
Figure 13. Modelling of the spatial distribution of the first independent component within the red rectangle in Figure 1. The spatial distributions are given in centimetres of displacement from the beginning of 1993 to the beginning of 2005, which corresponds to when the related temporal function underwent its maximum variation. (a,b) Horizontal and quasi-vertical displacements, redrawn from Figure 10a,b for clarity. (c,d) Computed cumulative contribution of the sill-shaped expansive source near Rione Terra and the expansive source below Solfatara. The centre of the former is marked by a red circle, while the latter is marked by a yellow square. (e,f) Residuals of the model, i.e., the differences between the spatial distribution of the first independent component and the computed cumulative contribution of the two sources.
Figure 13. Modelling of the spatial distribution of the first independent component within the red rectangle in Figure 1. The spatial distributions are given in centimetres of displacement from the beginning of 1993 to the beginning of 2005, which corresponds to when the related temporal function underwent its maximum variation. (a,b) Horizontal and quasi-vertical displacements, redrawn from Figure 10a,b for clarity. (c,d) Computed cumulative contribution of the sill-shaped expansive source near Rione Terra and the expansive source below Solfatara. The centre of the former is marked by a red circle, while the latter is marked by a yellow square. (e,f) Residuals of the model, i.e., the differences between the spatial distribution of the first independent component and the computed cumulative contribution of the two sources.
Remotesensing 16 01717 g013
Figure 14. Marginal probability density functions (MPDFs) of the parameters of the three pressurised spheroidal sources that satisfy IC2. MPDFs were estimated using NAB and scaled to a maximum of 1 to facilitate plot interpretation. The red curves pertain to the deep contracting source below the southwestern part of Campi Flegrei. The blue curves relate to the deep expansive source below Rione Terra. The green curves pertain to the deep contracting source below Vesuvio. The dashed lines relate to the parameter values listed in Table 3. (a) Easting position of the centres of the two sources below Campi Flegrei. (b) Easting position of the centre of the source below Vesuvio. (c) Northing position of the source centres. (d) Depth to the source centres. (e) Volume change of the sources. (f) Polar-to-equatorial radius (i.e., aspect ratio) of the sources.
Figure 14. Marginal probability density functions (MPDFs) of the parameters of the three pressurised spheroidal sources that satisfy IC2. MPDFs were estimated using NAB and scaled to a maximum of 1 to facilitate plot interpretation. The red curves pertain to the deep contracting source below the southwestern part of Campi Flegrei. The blue curves relate to the deep expansive source below Rione Terra. The green curves pertain to the deep contracting source below Vesuvio. The dashed lines relate to the parameter values listed in Table 3. (a) Easting position of the centres of the two sources below Campi Flegrei. (b) Easting position of the centre of the source below Vesuvio. (c) Northing position of the source centres. (d) Depth to the source centres. (e) Volume change of the sources. (f) Polar-to-equatorial radius (i.e., aspect ratio) of the sources.
Remotesensing 16 01717 g014
Figure 15. Modelling of the spatial distribution of the second independent component. The spatial distributions are given in centimetres of displacement from the end of 1999 to mid-2010, which corresponds to when the related temporal function underwent its maximum variation. (a,b) Horizontal and quasi-vertical displacements, redrawn from Figure 11a,b for clarity. (c,d) Computed cumulative contribution of the deep contracting source below the southwestern part of Campi Flegrei, the deep expansive source below Rione Terra, and the deep contracting source below Vesuvio. Based on several tests, the red circle bounds the area where the centre of the first source could be located. The centre of the second source is marked by a blue circle, and the centre of the third source is marked by a yellow square. (e,f) Residuals of the model, i.e., the differences between the spatial distribution of the second independent component and the computed cumulative contribution of the three sources.
Figure 15. Modelling of the spatial distribution of the second independent component. The spatial distributions are given in centimetres of displacement from the end of 1999 to mid-2010, which corresponds to when the related temporal function underwent its maximum variation. (a,b) Horizontal and quasi-vertical displacements, redrawn from Figure 11a,b for clarity. (c,d) Computed cumulative contribution of the deep contracting source below the southwestern part of Campi Flegrei, the deep expansive source below Rione Terra, and the deep contracting source below Vesuvio. Based on several tests, the red circle bounds the area where the centre of the first source could be located. The centre of the second source is marked by a blue circle, and the centre of the third source is marked by a yellow square. (e,f) Residuals of the model, i.e., the differences between the spatial distribution of the second independent component and the computed cumulative contribution of the three sources.
Remotesensing 16 01717 g015
Figure 16. Schematic three-dimensional view of the Campi Flegrei and Vesuvio magmatic systems and their sources during the 2002–2010 period. The bluish regions contracted, while the reddish regions expanded. The sill-shaped deformation source at 3–4 km below Campi Flegrei is labelled as S1. The mini-uplifts of 2000 and 2005 resulted from the expansion of the sole S1, whose contraction had previously caused the subsidence of the caldera. S1 continued to contract slowly until 2005, when it began to expand. The partially molten layer beneath Campi Flegrei and Vesuvio is represented by a light orange parallelepiped and is fed by a deeper reservoir. During 2002–2010, magma was transferred from beneath Vesuvio to beneath Campi Flegrei at a depth of 8–9 km, as indicated by the blue-to-red horizontal arrow in (a,b). The two panels depict two different scenarios. (a) If the spatial pattern of IC2 is not affected by IC1, two volumes located below Campi Flegrei at depths of approximately 9 and 8 km experienced opposite activity: contraction (beneath the southwestern part of the caldera) and expansion (beneath the central part of the caldera). (b) If the spatial pattern of IC2 is affected by contamination from IC1, it is more likely that a region beneath the eastern part of Campi Flegrei has expanded.
Figure 16. Schematic three-dimensional view of the Campi Flegrei and Vesuvio magmatic systems and their sources during the 2002–2010 period. The bluish regions contracted, while the reddish regions expanded. The sill-shaped deformation source at 3–4 km below Campi Flegrei is labelled as S1. The mini-uplifts of 2000 and 2005 resulted from the expansion of the sole S1, whose contraction had previously caused the subsidence of the caldera. S1 continued to contract slowly until 2005, when it began to expand. The partially molten layer beneath Campi Flegrei and Vesuvio is represented by a light orange parallelepiped and is fed by a deeper reservoir. During 2002–2010, magma was transferred from beneath Vesuvio to beneath Campi Flegrei at a depth of 8–9 km, as indicated by the blue-to-red horizontal arrow in (a,b). The two panels depict two different scenarios. (a) If the spatial pattern of IC2 is not affected by IC1, two volumes located below Campi Flegrei at depths of approximately 9 and 8 km experienced opposite activity: contraction (beneath the southwestern part of the caldera) and expansion (beneath the central part of the caldera). (b) If the spatial pattern of IC2 is affected by contamination from IC1, it is more likely that a region beneath the eastern part of Campi Flegrei has expanded.
Remotesensing 16 01717 g016
Table 1. Campi Flegrei multilayered elastic model: a linear interpolation is used to obtain values for intermediate depths between adjacent listed depths.
Table 1. Campi Flegrei multilayered elastic model: a linear interpolation is used to obtain values for intermediate depths between adjacent listed depths.
Depth
(km)
Vp
(km/s)
Vs
(km/s)
Density
(kg/m3)
≥0.001.600.921800
≥0.622.501.442100
≥1.403.201.852270
≥1.553.902.252380
≥2.733.952.282400
≥3.925.203.002580
≥4.035.923.422700
Table 2. Parameters of the two pressurised spheroidal sources that satisfy IC1 at Campi Flegrei. The parameters here are rounded, and their estimated marginal distribution can be seen in Figure 12.
Table 2. Parameters of the two pressurised spheroidal sources that satisfy IC1 at Campi Flegrei. The parameters here are rounded, and their estimated marginal distribution can be seen in Figure 12.
Source 1Source 2
Easting (m)425,900428,000
Northing (m)4,519,1004,519,700
Depth (m)31001300
Volume change ( 10 6  m3)4.91
Polar-to-equatorial radius0.018
Equatorial radius (m)2400
f0.5
Table 3. Parameters of the three pressurised spheroidal sources that satisfy IC2. The parameters here are rounded, and their estimated marginal distribution can be seen in Figure 14.
Table 3. Parameters of the three pressurised spheroidal sources that satisfy IC2. The parameters here are rounded, and their estimated marginal distribution can be seen in Figure 14.
Source 1Source 2Source 3
Easting (m)421,000426,000453,800
Northing (m)4,517,0004,519,4004,518,200
Depth (m)880077009200
Volume change (107 m3)−43.7−0.5
Polar-to-equatorial radius340.3
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Amoruso, A.; Gualandi, A.; Crescentini, L. On the Initial Phase of the Ongoing Unrest at Campi Flegrei and Its Relation with Subsidence at Vesuvio (Italy). Remote Sens. 2024, 16, 1717. https://doi.org/10.3390/rs16101717

AMA Style

Amoruso A, Gualandi A, Crescentini L. On the Initial Phase of the Ongoing Unrest at Campi Flegrei and Its Relation with Subsidence at Vesuvio (Italy). Remote Sensing. 2024; 16(10):1717. https://doi.org/10.3390/rs16101717

Chicago/Turabian Style

Amoruso, Antonella, Adriano Gualandi, and Luca Crescentini. 2024. "On the Initial Phase of the Ongoing Unrest at Campi Flegrei and Its Relation with Subsidence at Vesuvio (Italy)" Remote Sensing 16, no. 10: 1717. https://doi.org/10.3390/rs16101717

APA Style

Amoruso, A., Gualandi, A., & Crescentini, L. (2024). On the Initial Phase of the Ongoing Unrest at Campi Flegrei and Its Relation with Subsidence at Vesuvio (Italy). Remote Sensing, 16(10), 1717. https://doi.org/10.3390/rs16101717

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop