We report the solution for the displacement field in a half-space domain due to a stress-free strain ε 0 within a spherical geometry source (Rinaldi et al., 2010). We assume that the observation points Q x Q , y Q , z Q are placed on the surface of the half-space domain. Since the model has an axial symmetry, a 2D axi-symmetric domain is used in cylindrical coordinates Q ′ ρ , ϑ , z Q , where ρ = x C − x Q 2   + y C − y Q 2 and ϑ = arctan ⁡ ⁡ y C − y Q / x C − x Q . The dependence on ϑ disappears because of the axial-symmetry condition. The axis of symmetry is vertical and passes through the center 0 , z C of the source. The displacement, in the radial U r and vertical U z components, respectively, generated by a spherical source at the observation point Q ′ ρ , z Q is given as follows (Mindlin, 1936; Mindlin and Cheng, 1950; Davies, 2003; Rinaldi et al., 2011): U r Q ′ = 1 + ν Δ V s 3 π R ¯ ′ 3 ρ , (8) U z Q ′ = 1 + ν Δ V s 3 π R ¯ ′ 3 z Q − z C . (9) Here, Δ V S = ε 0 V 0 S is the stress-free volume change, where V 0 S represents the initial volume of the sphere, ν is the Poisson’s ratio, and R ¯ ′ = ρ 2 + z Q − z C 2 is the Euclidean distance between the observation point Q ′ ρ , z Q and the center 0 , z C of the deformation source.

The horizontal deformation components, U x and U y , respectively, at the observation point Q in Cartesian coordinates can be easily derived as follows: U x Q = 1 + ν Δ V s 3 π R ¯ 3 x Q − x C , (10) U y Q = 1 + ν Δ V s 3 π R ¯ 3 y Q − y C . (11)

It is worth noting that formulations (8) and (9) are similar to the Mogi solution, which provides ground deformation generated by a spherical cavity under the action of an overpressure Δ P M on its wall, that reads as (Mogi, 1958) U r Q ′ = a 3 1 − ν Δ P M μ R ¯ ′ 3 ρ = 1 − ν Δ V ¯ M π R ¯ ′ 3 ρ , (12) U z Q ′ = a 3 1 − ν Δ P M μ R ¯ ′ 3 z Q − z C = 1 − ν Δ V ¯ M π R ¯ ′ 3 z Q − z C , (13) where Δ V ¯ M = π a 3 Δ P M μ is the actual volume change and a is the radius of the Mogi source. Because of this similarity, it is not possible to distinguish between the magmatic and thermo-poro-elastic processes using deformation data alone (Lu et al., 2002), since they engender the same deformation pattern. In the Mogi model, the actual volume change Δ V ¯ M represents the actual radial expansion of the source wall, whereas the volume change Δ V S = 3 1 − ν 1 + ν Δ V ¯ S , which appears in the thermo-poro-elastic equation (Eqs 8, 9), represents the stress-free volume change if the source could expand freely under the action of the stress-free strain ε 0 . The Mogi solution is identical to the solution for a thermo-poro-elastic spherical source embedded within an elastic medium if an overpressure Δ P S = 1 + ν 1 − ν μ Δ V S 3 π a 3 acts from within (Belardinelli et al., 2019).

Ground deformation induced by thermo-poro-elastic strain changes within a cylindrical source has been described by Wang (2000). The mathematical formulation exploits the analogy between thermo-poro-elastic deformation and gravity change problems. For the gravity changes, closed analytical solutions were derived but only at the observation point along the axis of symmetry of the cylindrical source (Telford et al., 1981). Benefiting from these solutions, recently Todesco (2021) estimated the expected poro-elastic vertical deformation at Campi Flegrei induced by a cylindrical-shaped hydrothermal system. However, this closed-form analytical solution only allows the computation of the vertical deformation at a point along the axis of symmetry. Semi-analytical solutions can be derived using spherical harmonic series of Legendre polynomials (Na et al., 2015). Examinations of both gravity (Na et al., 2015) and, more recently, thermo-poro-elastic deformation (Mantiloni et al., 2020) have shown that the solutions are accurate under the assumption that the height h of the cylinder is much smaller than its radius R ( h / R ≪ 1 ). The analytical solutions for thin disk-shaped sources (Mantiloni et al., 2020) can be easily generalized to model thick cylindrical sources by stacking two or more disks (Nespoli et al., 2021; Belardinelli et al., 2022).

We provide an alternative way to compute the displacement induced by pore-pressure and temperature changes for a cylindrical source by extending the formulation proposed by Hemmings et al. (2016) for the gravity changes. Starting from this formulation, we have derived the semi-analytical solution to compute the analogous radial and vertical displacements. In the cylindrical coordinate system, each elementary source ( ρ i , z i ), with which the source is composed of, represents a portion of an elementary source ring centered in ( 0 , z i ) and radius ρ i (Supplementary Figure S1). The distance between the observation point Q ′ and the elementary source at ( ρ i , z i ) is, in cylindrical coordinates, a function of ϑ , that is, s ϑ = ρ i 2 − 2 ρ i ρ cos ϑ + ρ 2 + z Q − z i 2 . Under this configuration, from Eq. 7, we derived the displacement potential ϕ = 1 + ν 3 π ε 0 i ρ i d ρ i d z i ∫ 0 2 π 1 s ϑ d ϑ for a single ring and obtained the displacement field at Q ′ by computing the gradient of the potential (Eq. 6): U i r Q ′ = 1 + ν 3 π ε 0 i ρ i d ρ i d z i ∫ 0 2 π ρ − ρ i cos ϑ s ϑ 3 d ϑ , (14) U i z Q ′ = 1 + ν 3 π ε 0 i ρ i d ρ i d z i z Q − z i ∫ 0 2 π 1 s ϑ 3 d ϑ , (15) where V i = ρ i d ρ i d z i d ϑ is the volume of each elementary source ( ρ i , z i ) on the ring and G r ρ i , z i = 1 + ν 3 π ρ i d ρ i d z i ρ − ρ i cos ϑ s ϑ 3 and G z ρ i , z i = 1 + ν 3 π ρ i d ρ i d z i z Q − z i s ϑ 3 are the Green’s functions. Equations 14, 15 represent, respectively, the radial and vertical components of the displacement field due to a single ring. The integrals in the equations are solved numerically using a recursive adaptive Simpson quadrature scheme.

Finally, the displacement, recorded at the observation point Q ′ , is obtained by summing up the contributions of each elementary ring: U r Q ′ = 1 + ν 3 π ε 0 ∑ ρ i , z i ρ i d ρ i d z i ∫ 0 2 π ρ − ρ i cos ϑ s ϑ 3 d ϑ , (16) U z Q ′ = 1 + ν 3 π ε 0 ∑ ρ i , z i ρ i d ρ i d z i z Q − z i ∫ 0 2 π 1 s ϑ 3 d ϑ . (17)

Equations 16, 17 represent the semi-analytical formulations to compute the ground displacements generated by a cylindrical source consisting of many elementary rings as a function of the position of the observation point Q ′ ρ , z Q .

This approach only allows us to compute the surface displacements (i.e., z Q = 0 ). In order to compute displacement, strain, and stresses in the interior points of the half-space domain, it is no longer possible to use the scalar potential because the displacement field is not irrotational (Mantiloni et al., 2020).

In order to verify the results and check the accuracy of the numerical integration, solutions from Formulas 16, 17 were compared with the numerical thermo-poro-elastic results (Supplementary Material) calculated using COMSOL Multiphysics software (COMSOL, 2012), which solves Eqs 1–3 with a finite-element (FE) discretization (Currenti and Napoli, 2017; Stissi et al., 2021).

The accuracy of the semi-analytical solution depends on the discretization of the source in elementary rings with finite thickness. The smaller the element, the better the solution. For shallow and/or large cylindrical sources, a more accurate solution is obtained by considering the cylindrical source as the superposition of smaller sources (Supplementary Figure S3). The surface displacement is calculated, according to Eqs 16, 17, as the sum of the displacements induced by the single smaller sources. The comparison shows a good agreement between the proposed semi-analytical and numerical solutions (Supplementary Figures S2, S3).

In the following section, we applied the derived solutions to invert the ground deformations observed at the onset of the 2021 unrest episode in Vulcano Island.

Figure 2 shows the daily horizontal ( U x and U y and vertical U z deformation data measured by the GPS monitoring network at Vulcano Island from January 2019 to December 2021. Raw daily GPS data were processed using GIPSY software (Bertiger et al., 2020) in a precise point positioning mode applied to the ionospheric-free carrier phase and using JPL’s final orbit and clock products. Solutions were first aligned to IGS14 by applying a daily seven-parameter Helmert transformation. Finally, to remove the, albeit minor, regional long-term tectonic movements, each time series of local positions was corrected from residual trends estimated up to August 2021 as a linear trend. To reduce the noise, a moving centered median filter was applied with a fixed window length of 35 days. No significant variations were observed before September 2021 after which the crisis started with sudden increases in gas emission (Aiuppa et al., 2022) and seismicity (INGV Report, 2022; Federico et al., 2023). Concurrent ground deformations of a few centimeters were recorded at almost all the stations from the beginning of September until mid-October, when the deformation ceased. The average deformation rate at the closest station to the summit (IVCR) was approximately 4.06 x 10⁻² cm/day. Data disclosed a clear radial pattern centered in the La Fossa Cone area (Figure 3A). The ground deformation observed between 2 September and 13 October indicated a continuous expansion of the volcano edifice. During the entire period, a continuous amplification of the deformation was recorded at all the stations with a persistent pattern (Figure 3A). Indeed, the ratio between ground deformation components at different stations was fairly constant (Figure 3B), suggesting that the deformation source, which was observed from the onset of the Vulcano Island crisis, is almost spatially stationary.

The horizontal and vertical displacements observed at Vulcano Island have been inverted to constrain the source and gain insights into the deformation process. We explore both the spherical and cylindrical sources in order to find the best-fitting solution.

The inverse method combines the derived forward models with GA (Tiampo et al., 2000; Currenti et al., 2007; Carbone et al., 2008), in order to find the model parameters that minimize the misfit between the observed and computed deformations. Following an evolutionary scheme, GA iteratively explores the model parameter space and tries to find the global optimal solution. The misfit is quantified using an objective function defined as the root mean square error (RMSE) between the observed data U obs = ( U obs x , U obs y , U obs z ) and the deformation computed by the forward model U calc = ( U calc x , U calc y , U calc z ): RMSE = 1 n ∑ i = 1 n 1 U o b s x i − U c a l c x i 2 + ∑ i = 1 n 2 U o b s y i − U c a l c y i 2 + ∑ i = 1 n 3 U o b s z i − U c a l c z i 2 . (18) Here, n = n 1 + n 2 + n 3 , where n 1 , n 2 , and n 3 represent the number of observation points for each deformation component.

For the spherical source, the model parameter vector is represented by p = x C , y C , z C , Δ V S , where x C , y C ,   z C ) are the coordinates of the center of the source and Δ V S is the stress-free volume change, while for the cylindrical source, it is represented by p = { x C , y C , d , R , h , ε 0 , where x C ,   y C are the coordinates of the center of the source, d is the depth of the top of the cylinder, R is the radius, h is the height, and ε 0 is the stress-free strain (Tables 1, 2).

For each set p of parameters, the forward model computes an initial population of solutions U calc p whose fit with the observed data is evaluated by means of the objective function (Currenti et al., 2007) to evaluate the “best” sets of parameters in the entire population, i.e., the model that minimizes the objective function. The best sets of parameters are modified by applying evolutionary rules for the generation of a new set of parameters that on average are better than the previous parameters. The procedure is iterated, and the initial population evolves over several generations (Currenti et al., 2005), until the algorithm converges to the global optimal solution.

It is worth noting that the deformation solutions (Eqs 8, 9 and Eqs 16, 17) are linearly proportional to Δ V S for the spherical source and ε 0 for the cylindrical source. For this reason, the parameters Δ V S and ε 0 are omitted in the set of the inversion parameters and are computed at each step by the least squares method: Δ V S  or  ε 0 = U calc × U calc T − 1 × U calc × U obs T . (19)

In such a case, the dimension of the search domain for optimization is reduced to p = x C , y C , z C and p = { x C , y C , d , R , h for the spherical and cylindrical geometries, respectively.

An extensive search was performed on the model parameter space, whose ranges are reported in Tables 1, 2. Since the deformation pattern clearly points to a deformation source located below the La Fossa Cone, the search range for the source position is limited within a 4-km-wide box centered in the summit area. GA is initialized with a random population consisting of 100 individuals, and it iterates till it converges. A single GA inversion takes on average 3.5 s for the sphere and 12 s for the cylinder. The GA inversion is performed 5,000 times to obtain an estimate of the model uncertainty. The models obtained after the convergences are used to appraise the results by computing the 1D and 2D marginal distributions. The 1D marginal distributions, given by the histograms of the model parameters in the solution set, provide the confidence intervals, whereas the 2D marginal distributions, calculated for selected pairs of parameters, offer further information about their trade-offs (Sambridge, 1999).

The computed displacements for the optimal spherical solution (Figure 4), whose parameters are reported in Table 1, generally fit with the observed data with an RMSE of approximately 3 mm. The solution indicates a deformation source centered in the La Fossa Cone, as already suggested from the radial pattern of the data (Figure 4A), at an average depth of 720 m below the ground surface. The 1D marginal distributions (Figure 5) indicate that the model parameters are well constrained. The mean, the median, and the optimal solutions are well within the 95% confidence intervals which are very narrow for all the parameters (Table 1). The 2D marginal distributions of the estimated source parameters from the inversions show that they concurrently converge toward the optimal solution.

For the cylindrical source, the computed displacements of the optimal solution (Figure 6; Table 2) are very similar to those obtained for the spherical source (Figure 4) with a comparable RMSE (3 mm) (Tables 1, 2). The source position is almost the same, and the 95% confidence intervals of the two solutions show considerable overlap. The top of the source lies at 360 m depth providing a cylinder center depth of approximately 860 m, which is almost the same center depth estimated for the spherical source (720 m). Nonetheless, the cylinder center depth is slightly outside the confidence interval of the sphere depth, which ranges between 646 and 849 m. The radius R and height h parameters are not well constrained (Figure 7A) and have larger confidence intervals. The examination of the 2D marginal distributions shows that the radius and height parameters are not fully independent (Figure 7A). The smaller the radius, the taller the cylinder. The marginal distributions show that these parameters have a wide range of variability, and, starting from the optimal solution, there are several pairs of parameters ( R , h ) values and associated stress-free strain ε 0 values, for which the computed deformation fits the observations as well (Supplementary Figures S5, S6). This implies that diverse cylindrical sources with different radii, heights, and strain variations can be equally responsible for the observed ground deformations.

Under the action of thermo-poro-elastic effects, both the optimal spherical and cylindrical sources undergo a comparable volume change. For the spherical source, an optimal value of Δ V S = 1.07 x 10⁵ m³ is obtained. For the cylindrical source, the volume change is also well constrained (Figure 7C). Its distribution assumes values between 0.9 x 10⁵ m³ and 1.4 x 10⁵ m³ centered around a value of Δ V C = ε 0 V 0 C = 1.15 x 10⁵ m³ (Figure 7C), which is very similar to that of the spherical source.

Despite the reasonable and similar fits of the spherical and cylindrical source models, discrepancies between the observed and computed deformations are obtained, for both geometries, at VCSP, IVUG, and IVLT stations. At VCSP and IVLT stations, although the amplitudes of the computed deformations are comparable to those observed, the computed radial deformation vector appears to be rotated by approximately 40° clockwise with respect to the observations. Moreover, at IVUG, the model predicts a strong attenuation. We further investigated whether these discrepancies could be explained by the volcano topography, with an edifice that extends from −1,180 m b.s.l. to 497 m a.s.l. We performed finite-element modeling for the spherical source using COMSOL (Supplementary Figure S7). The COMSOL numerical results are similar to the analytical solutions. The only slight difference is observed at the IVCR station, closest to the La Fossa cone, where the maximum radial deformation is observed. In fact, at this station, the COMSOL deformation solution appears rotated clockwise with respect to both the semi-analytical solution and the observed data. Overall, the RMSE is approximately 3 mm, similar to that achieved for the analytical solution. Therefore, the surface topography alone is not sufficient to justify the rotation at VCSP and IVLT and the larger displacement at IVUG.

A sensitivity analysis is performed to evaluate the impact of the parameters on the solutions of the spherical and cylindrical sources. We applied the Morris method (Morris, 1991), a widely used one-at-a-time (OAT) global sensitivity analysis (GSA) (Feng et al., 2019). In particular, the Morris method represents an excellent tool for identifying the influential parameters of a model with a large number of inputs and quantifying the response of the model to the change in the input parameters (Conca et al., 2016; Liu et al., 2020).

Starting from a sample of the input parameters x 1 * , … , x i − 1 * , x i * , x i + 1 * , … , x n * , a trajectory is built within the input space by modifying one parameter at a time x 1 * , … , x i − 1 * , x i * ± Δ , x i + 1 * , … , x n * and evaluating, for each variation, the elementary effect E E i j ( i th parameter and j th trajectory) given by the ratio between the variation of the output solution y and the variation Δ of the input: E E i j = y x 1 * , … , x i − 1 * , x i * + Δ , x i + 1 * , … , x n * − y x 1 * , … , x i * , … , x n * Δ  or  E E i j = y x 1 * , … , x i * , … , x n * − y x 1 * , … , x i − 1 * , x i * − Δ , x i + 1 * , … , x n * Δ . (20)

Having iterated this procedure for a r number of trajectories for each parameter i , the absolute value of the mean μ i * and the standard deviation σ i of the variations is computed, respectively, as follows: μ i * = 1 r   ∑ j = 1 r E E i j , (21) σ i = 1 r − 1   ∑ j = 1 r E E i j − μ i 2 , (22) where μ i = 1 r   ∑ j = 1 r E E i j .

These quantities provide a measure of the sensitivity of the model. High values of μ i * indicate that the i th parameter has a great influence on the output solution.

The results of the Morris method applied to our models are shown in a ( μ * , σ ) plot (Figure 8). For the cylindrical source model, the results are in agreement with the marginal distributions of the GA solutions set (Figure 7A). The points of the plot in the lower left ( R and h ) have smaller mean μ * values, indicating that a variation in these parameters causes negligible effects on the output (Conca et al., 2016). Indeed, these parameters show a wide range of variability in the marginal distributions (Figure 7A). On the other hand, the points at the top right of the plot ( x C , y C , and d ) show larger values of the μ * , which denotes that a variation in these parameters strongly influences the output of the model (Conca et al., 2016). In fact, these parameters are well constrained as shown in the 1D marginal distributions (Figure 7A).

For the spherical source model, the mean values μ * of the parameters do not differ significantly from each other. All the parameters are well constrained and equally concur with the solution as already seen from the inspection of the marginal distributions.

The thermo-poro-elastic deformation is the result of the volumetric strain variations linked to changes in pore-pressure and temperature (Eq. 4). Thermal effects can be much greater than those due to pressure (e.g., Nespoli et al., 2021; Belardinelli et al., 2022) but are usually much slower than the pore-pressure buildup (Coco et al., 2016; Currenti and Napoli, 2017). However, thermal effects may be faster when advection processes develop. In the 2021 Vulcano Island crisis, although the thermo-elastic effect could have also contributed to the deformation (temperature increase up to 50°C; INGV Report, 2022), we focused on the estimation of pore-pressure increase. Assuming that the observed deformation is driven only by the pore-pressure variation at the source, we can estimate the pressure change from the relationship between Δ P and ε 0 (Eq. 4). This relationship suggests that the pressure variation also depends on the elastic properties of the rock, which are closely related to lithology, rock type, degree of fracturing, water content, depth, confining pressure, temperature, compaction, and hydrothermal alteration (Heap et al., 2014; Juncu et al., 2019; Todesco, 2021). In the literature, no data are available on the elastic properties of Vulcano rocks. For this reason, we refer to the values reported in similar volcanic environments, according to which K s generally does not exceed 30 GPa (Rinaldi et al., 2010; Todesco et al., 2010; Currenti and Napoli, 2017; Juncu et al., 2019), while K d and β vary from ∼10⁻¹ to ∼10 GPa and from 0.5 to 1, respectively. In order to satisfy this last condition, we have chosen K d in the range between 0.1 and 15 GPa.

For the spherical source, the inversion enabled to constrain the source volume change Δ V S = 1.07 x 10⁵ m³, from which we can derive the stress-free strain ε 0 = Δ V S / V 0 S using reasonable values of the source radius. Using the end-member values of the elastic parameters and a radius ranging from 500 to 700 m, which corresponds to a volume V 0 S between 5.2 x 10⁸ and 1.4 x 10⁹ m³, the pore-pressure change varies approximately from 0.01 MPa to 6 MPa.

For the cylindrical source, the stress-free strain ε 0 is not well constrained. Solutions with larger ε 0 and lower source volumes (Figure 7B) fit the observations as well (Supplementary Figure S8). However, larger ε 0 values correspond to unrealistic pore-pressure changes. In addition, the optimal solution (Table 2) could require large overpressure (0.3–80 MPa for lower and higher K d values). A compromise is found by selecting those solutions whose source volume V 0 C ranges between 5 x 10⁸ and 9.5 x 10⁸ m³. This solution leads to an increasing attenuation of the deformation at the IVCR station as the radius increases, causing similar radial displacements at the IVCR and VCSP stations. As an example, Figure 9 shows the displacements induced from a cylindrical source, with a radius of 400 m and a volume of 5 x 10⁸ m³, which undergoes a volume variation Δ V C of 1.25 x 10⁵ m³ and a pressure change from 0.02 to 7 MPa for lower and higher K d values, respectively.