Dipole Source to Model Ground Deformation Due To Magma Migration

Although magma migration is a necessary condition for the occurrence of a volcanic eruption, static analytical models for ground deformations consider stationary sources that undergo a change only in forces acting in the Earth's interior, such as changes in pressure or volume in a geometrically defined body. The variations in pressure or volume are mostly linked to the fluid dynamics or the superficializing of magma batches. Detection and tracking of magma motion in volcanic areas represent fundamental capabilities in predicting possible eruptions and therefore in reducing their associated risks. In the 3‐dimensional space, a magma migration defines a starting volume in depression and an arrival volume in overpressure. These end members can be modeled with two bonded point sources. Such a double‐point source can be considered the equivalent of a physical dipole. Here, the analytical solutions to predict ground deformations caused by a dipole source are presented and analyzed. The dipole analytical source allows modeling the effects on ground deformation of magma migration and vice versa. Estimating the magma migration from the measured ground deformation through data inversion. The dipole solution is compared with the basic volumetric source to discern its distinctive peculiarities and range of applicability. In addition, a discrepancy analysis with the equivalent numerical model was conducted to evaluate the approximation errors of the analytical model. To evaluate the practical value of the proposed model, two case studies on Sinabung volcano and Mt. Etna are presented. The analysis of continuous deformation data from GNSS (Global Navigation Satellite System) stations on Mt. Etna shows that a dipole source activates days before the December 2018 eruption, thus allowing the identification of the final magma rise.


Introduction
Understanding the structures and the dynamics of magma plumbing systems is critical for predicting volcanic eruptions and assessing volcanic hazards (Scarpa & Tilling, 2012).Ground deformations around volcanoes come to the aid of this objective by reflecting buried volcanic processes occurring in the plumbing system that are transmitted to the surface through the crust (Sigmundsson et al., 2018).The magma motions triggering an eruption are included among these processes.It follows that measurements of ground deformation, either stationbased (e.g., GNSS (Global Navigation Satellite System) stations, tiltmeters, strainmeters) or remote sensing (e.g., InSAR satellite radar) (Dzurisin, 2007), contain precious information to reconstruct the processes occurring beneath a volcano.
Ground deformation in active volcanic areas usually follows cycles characterized by the replenishment/depletion of magma reservoirs.The common assumption for such functioning model is the presence of a deeper and a constant magma source (Reverso et al., 2014).
Analytical source modeling for volcanic ground deformation has shown its value in investigating magma plumbing systems due to its ability to provide insights into the underlying physical processes that govern magma transport and eruption.
Classic volcanic analytical source models predict ground surface effects from forces acting in the Earth's interior, such as changes in pressure or volume in geometrically defined bodies (Dzurisin, 2007;Segall, 2013).These sources include bodies of various shapes, and the typical associated dynamics refer to inflation or deflation for the volumetric sources, and opening or closing for dyke-like sources.The diversity of volcano deformation processes requires a flexible representation of the wide variation in source geometries (Nikkhoo et al., 2016).The need to still develop analytical models in elastic media lies in their ability to give quick answers and concise information from partial knowledge with reduced complexity.In this context, the multitude of volcano deformation processes requires further analytical models capable of covering this variety (Nikkhoo et al., 2016).
Most of the classic ground deformation analytical models are defined "static" models though all of them intrinsically take into account the time in one or more parameters that represent a variation (e.g., volume change ΔV, pressure change ΔP, strike slip, etc.).Any parameter describing a variation happens in time, thus the "static" adjective typically refers only to the fact that the time variable is not explicit in the equations.In addition, it is interesting to note that the static nature of the most widely used analytical models also refers to their position in space.Bearing in mind that any data collected in different times describes a dynamic process, in static modeling dynamic data (i.e., changes in certain parameters over time) are used to fit a source model that excludes any possible dynamic change in force location within the period considered by the data.In other words, the variation in source position is not expressed in the formulations of traditional analytical models.Conversely, magma (hereinafter, the term magma will be used in the generic sense of fluid, gas or multi-phase mass able to create a pressure gradient) moves either to accumulate in some chamber or to reach the surface in the case of eruptions.The static analytical models contrast with the dynamics of magma motions.Even a simple magma accumulation needs a rising of some magma batch.To reach the surface, magma has to move, and this movement, associated with pressurization/depressurization dynamics due to possible multiple physico-chemical processes, can deform the ground surface.In an analytical source describing a magma confined chamber, the variation parameters (e.g., volume/pressure changes) can be associated to magma transfer only if the other end-side (source or hole) is much deeper and its effects at the surface are negligible.In all the other cases, the variation in some parameters can be due only to internal changes (e.g., magma differentiation, gas exsolution) or medium changes (e.g., stress change, pore pressure change).In this context, none of the classic analytical deformation sources considers the dynamic change in magma position.Therefore, a possible variation in the position of a single magmatic source is usually not deduced with a single measurement of the surface deformation field.Instead, the source location variation is typically inferred by considering and statically inverting multiple measurements in time (Bonaccorso & Aloisi, 2021;Cannavo' et al., 2019;Cannavò et al., 2015;González et al., 2013).Magma accumulation implies magma motion, and for a finite magma quantity, a motion implies a simultaneous increasing in pressure in a volume and a decreasing in pressure in another volume.This characteristic of magma transfer/motion is what the dipole source aims to model.In other words, the model describes processes that simultaneously and differentially involve two points in space, mimicking a magma transfer between the two points.Some models with different complexities already exist (Kozono, 2021;Melnik & Costa, 2014;Reverso et al., 2014;Sigmundsson et al., 2018) that consider articulated magma plumbing systems with magma transfers in multi-chamber plumbing systems through conduits.In the case of the Soufrière Hills volcano on Montserrat, for example, magma transfer between two chambers connected by a vertical conduit was inferred over a long period using magma efflux data and surface deformation data obtained through GPS during both eruptive and inter-eruptive periods (Elsworth et al., 2008).For the same volcano and with the same double chamber plumbing system, Melnik and Costa (2014) created a numerical model to analyze the possible volcano dynamics.
Either analytical or numerical, the proposed models in literature consist of dynamic models with explicit time variable, for which by setting constraints, parameters and initial variables, they return the time evolution of certain quantities.What is lacking is a simple, albeit approximate, analytical tool that can identify and coarsely model the magma transfer process causing ground deformation.
Here, an analytical solution for surface deformation due to migration of a volumetric source is presented.The formulation is valid for small movements of the volumetric source in an elastic half-space.The proposed model is described by an oriented axis of motion that connects the initial and final positions of the source, which can be likened to a volumetric dipole.The two vertices that define the physical dipole are two spherical point sources of radial symmetry.A dipole source, though approximate, supplies possible footprints in ground deformation data for magma transfer/motion simultaneously linked to an overpressure (or volume change).This allows modeling the effects on ground deformation of magma migration and vice versa estimating the magma migration from the measured ground deformation through data inversion.In the following, after presenting the model development and analysis, two case studies are illustrated to identify a dipole source at Sinabung volcano and the final phase of an unrest at Mt. Etna.

Methods
An electric dipole consists of two equal charges separated by a small distance.The associated electric dipole moment is a vector quantity that points from the negative charge to the positive charge and has a magnitude given by the product of the charge magnitude and the separation vector between the charges.The electric field due to an electric dipole is directly proportional to its moment and inversely proportional to the cubic of the distance of the measuring point.
Similarly, a dipole magma source can be described as a differential source, defined as consisting of two point sources of equal intensity and opposite sign, separated by a small distance (see the sketch in Figure 1).
Assuming an isotropic elastic half space, a pressure change ΔP within an embedded spherical cavity with radius α much smaller than its depth d (i.e.α ≪ d) produces a surface displacement given by Equation 1(also known as Mogi's model) (Delaney & McTigue, 1994;Dzurisin, 2007;Mogi, 1958): where G is the shear modulus, ν is the Poisson ratio, and ‖r‖ = x 2 + y 2 + d 2 ) 1 2 is the norm of the row vector r = [x, y, d] that is the Euclidean distance between the point source located at d > 0 depth on the origin coordinates, and the probing point at the ground surface with coordinates (x, y).The equivalent cavity volume change ΔV can be approximated with though due to magma compressibility, ΔV may not reflect the volume of actual moving magma.Thus, Equation 1 can be rewritten as: It has been shown that analogous solutions of the deformation field given by a point source of pressure can be similarly derived from other different conceptual models such as three orthogonal force dipoles or three orthogonal tensile dislocations provided that the parameters are properly chosen (Bonafede & Ferrari, 2009).
Given Equation 3, in analogy with an electric dipole, the dipole moment m for a pair of opposite point sources can be defined as a vector representing the source strength and directed along the dipole axis which is a line directed from the negative pressure variation point to positive pressure variation point.Let us define δ in Equation 4 as the dipole axis (i.e., the vector distance joining the two sources) whose length is ‖δ‖ = δ (i.e., the distance between the two sources), θ is the azimuth angle (clockwise from North) of the dipole axis directed from the deflation source to the inflation one and φ is the elevation angle measured from the vertical (0°at the Zenith, 90°for horizontal dipole).The dipole moment can be expressed as: where ΔV is the volume change of each of the two point sources.The negative sign in Equation 5 accounts for the reference system of the source with vertical coordinate positive downward.The dipole strength is given by the module of the moment ‖m‖ = ⃒ ⃒ ΔV 1 ν π δ ⃒ ⃒ .
In the hypothesis of non interacting sources, the surface displacements at a point c produced by a dipole source can be approximated by the superposition of the effects of the two point sources separately, and adding the hypothesis of sufficiently small δ it is possible to approximate the formulation as shown in Equation 6.
where u + p (c) and u p (c) represent respectively the displacements produced by the positive and negative poles of the dipole source located at a point c, which consist of two point sources located at the extremes of the dipole axis (i.e., u p (c ± δ/2) in the second equivalence).
Defining ɛ = δ/d as the fraction of the dipole axes length over its depth, the above feasibility assumptions are equivalent to ɛ ≪ 1.In the text, this assumption will be referred to as the small movements hypothesis.
The limitations of the hypothesis on non-interacting sources will be analyzed in Section 3.2 devoted to comparison with results obtained by numerical finite element methods.
To make ∇u p explicit it is possible to use Equation 3 and calculate: where δ ij is the Kronecker delta.In analogy to the electric field, the deformation field induced by a dipole source, whose barycentre is located at depth d, can be expressed in terms of its dipole moment m using Equations 3 and 5-7 to derive Equation 8: where r = r/‖r‖ are the relative normalized coordinates.Equation 8 can therefore be considered the dipole equation for calculating the related deformation field.
On the other hand, we define the tilt T as the angle variation along the coordinate axes which is reckoned positive when the vertical displacement decreases in the x or y direction.The surface tilt components in Equation 9 are obtained from Equations 8 and 7 by taking the negative of the partial derivative of the vertical displacement u dp v with respect to the x and the y coordinate.

Generalized Dipole
In the case where the two sources at the poles display different strengths as they move away from the dipole center, thus unbalancing the dipole, it is necessary to consider an additional term to account for this asymmetry in Equation 8.In particular, in ideal conditions, assuming the same thermodynamic pressure-volume Earth and Space Science 10.1029/2023EA003003 CANNAVÓ work (i.e., W = PΔV) (Callen, 1960;Gudmundsson, 2014) in the two poles of a dipole source (W + = W ), an equal volume change deviation dV/2 with opposite sign in the two poles and a homogeneous space (thus pressure proportional to the depth, P ∝ d), the volume variation ratio can be written as: where dV is the difference between the volume variations at the two poles and δ z is the length of the vertical component of the dipole axis.The variation in ΔV is thus converted in a depth variation from Equation 10: However, during magma migration, volumes in general are not conserved depending on stiffness of the host rock, shape of the sources involved, magma compressibility, depth of the sources and external stress field (Rivalta & Segall, 2008).These and other many chemical and physical processes can make a dipole unbalanced.To account for this imbalance, we must abandon the previous assumption of equal volume change deviation at the two poles; so, maintaining the assumptions of conservation of thermodynamic pressure-volume work and pressure proportional to depth, it is possible to write W + = W as: where an increase in volume at lower pressures was taken as a reference and ξ ≠ 1 represents a dimensionless coefficient modulating the imbalance of the volume change deviation at the two poles.Equation 12 can be rearranged as following: For ξ away from the unfeasible 1, considering the assumption of d ≫ δ, from Equation 13 it can be derived that: that is a generalization of Equation 11 from which it differs by a proportionality factor.
Thus to generalize the dipole formulation it is possible to introduce a single coefficient κ = 1/(1 + ξ) that controls the volumetric change as the depth changes.For κ > 0 the shallower part of the dipole varies by a larger volume than the deeper part (κ = 1 in the ideal hypothesis of Equation 11).Conversely, for κ < 0 it is the deeper part that undergoes a greater volume change than the shallower part.
By differentiating Equation 3 for ΔV, considering Equation 14, the defined κ, and considering that the change in volume is positive at negative change in depth, the asymmetric general dipole model due to differential depth becomes: where êz is the row unit vector of the vertical axes.
For κ = 0 the formulation collapses into the symmetric dipole in Equation 8. κ < 0 indicates a reducing in strength (i.e., ΔV) for the shallower monopole.This makes it possible to contemplate cases where the deeper chamber deflation is predominant during the eruption in the dual-magma-chamber system as found in several eruptions of different volcanoes (Kozono, 2021).

Analysis of Solutions
The ideal dipole source comprises two point monopoles of equal strength and opposite polarity separated by a distance δ that is very small compared to the distance of considered effects.It is worth noting that Equations 8 and 15 show that the deformation strength at the surface is proportional to the product of volume change and dipole axis length, which becomes a source strength factor (δΔV) called hypervolume from here on.Because of this ambiguity, in an inversion process, it is not possible to retrieve ΔV and δ separately.
Figure 2 shows the deformation patterns generated by a dipole for different values of elevation angle (φ).In particular, Figure 2b refers to a symmetric dipole source (κ = 0 in Equation 15), while 2a refers to an asymmetric one (κ = 1 in Equation 15).In the latter, it is evident that the shallower dipole end source becomes predominant as the dipole becomes vertically oriented.
Figure 2c shows the deformation pattern in the case of asymmetric dipole with a negative κ, when the deeper dipole end becomes increasingly prominent.
From Figure 2 it is clear that a signature of a dipole source is the possible simultaneous presence of uplift and subsidence zones, thus an analysis on the extreme values that these deformations induce has been carried out.To simplify the calculation, but without losing generality, an x-oriented dipole source centered on the coordinate axis has been considered.The minimum and maximum vertical deformations can be located only along the x axis.
Figure 3 shows the normalized positions of the minimum (in blue) and maximum (in red) vertical displacements induced by a symmetric (κ = 0) dipole source rotating along the y axis with an angle given by the elevation angle (φ).The positions are normalized for the depth of the dipole source.
The line in black shows the absolute distance between the points with extreme vertical deformations.It can be seen that the maximum distance is reached for a pure vertical dipole (i.e., elevation 0°± n180°) and it equals twice times the depth location, while the minimum distance is reached for a perfectly horizontal dipole (i.e., elevation 90°± n180°) and it equals the depth location.As expected, in a completely horizontal dipole (i.e.side by side dipole extremes at same depth d) the zero vertical displacement is reached at the source epicenter (see green line in Figure 3).To better visualize the characteristics of surface deformations induced by a vertical symmetric dipole, all deformation components scaled with the maximum vertical deformation are shown in Figure 4.

Analysis of Discrepancy and Numerical Comparison
As is the practice in the literature when presenting a new model (e.g., Nikkhoo & Rivalta, 2023;Pascal et al., 2013), various comparisons are made in this analysis to check the deviation of the dipole model from known composed models of which this is a particular case.
In particular, a comparison with a two point-source superposition model is performed to evaluate the discrepancies of solutions (i.e., the differences in the predicted deformations) and the range of validity of the dipole model.The ground deformation predicted by a dipole source is compared with the predicted deformation due to two point sources (Mogi, 1958), identical but opposite in strength (ΔV), which are either superposed or juxtaposed with the same distance of the dipole axis (Pascal et al., 2013).
As expected, the behaviors of the two compared models are almost identical (see Figure 5), with a discrepancy at the peak values.
From Equations 3 and 15, it is possible to calculate analytically the difference between the solutions of the two models compared and evaluate this discrepancy in the results as the considered parameters change.By varying the average depth (d) of the sources and their distance (δ) it can be shown that the maximum relative discrepancy between the two models depends essentially on the ratio between the two-point source distance and their average depth (i.e., the ɛ ratio for the dipole).
Figure 6 shows the maximum relative discrepancies for the horizontal and vertical deformation components in the two extreme cases of vertically oriented (superposed) sources and horizontally oriented (juxtaposed) sources.It can be deduced that the horizontal component has the highest discrepancy that  This analysis therefore confirms that, if the feasibility assumption is met, the dipole model solution is valid compared to the point source analytical model (Mogi, 1958).

Finite Element Analysis Comparison
To further investigate the validity of the proposed dipole model, a comparison was made with the solutions obtained from a finite element model (FEM) by using the Partial Differential Equation Toolbox of MATLAB software (Inc, 2022) environment.
In particular, the solutions between the dipole model and the equivalent FEM model were compared when varying the dipole depth d and axis length δ parameters as done previously with the dual source analytical model.
In the case of the relatively more realistic FEM model, a homogeneous and isotropic elastic domain was considered large enough to avoid boundary effects; two spherical cavities of radius a located at the ends of the dipole and oppositely pressurized were built for each model.
For comparison with the realistic case of the FEM model, it was therefore necessary to also introduce the radius parameter a as a variable into the discrepancy analysis and to use the analytical formulation of the dipole moment described with the pressure variation (i.e., Equations 2 and 5).
To analyze the discrepancy, the maximum relative error was considered (i.e., the maximum error normalized for the maximum predicted deformation).

Earth and Space Science
10.1029/2023EA003003 CANNAVÓ As previously found in the comparison with the 2-source analytical model, by varying the average depth (d) of the sources and their distance (δ) and the radius a, it can be shown that the maximum relative discrepancy between the two models depends essentially on the two ratios: the ɛ and the ratio between radius and axis length (i.e., a/δ).In the formulation with only the volume change, the assumption is reduced to only that on ɛ.
Figure 7a shows the results of this comparison for a vertical dipole (2 superimposed cavities in FEM).
It is possible to see that for a/δ < 0.3 and ɛ < 0.3 the discrepancy between the dipole vertical solution and the realistic FEM vertical solution the discrepancy remains well below 10%.
Figure 7b shows the results of the same comparison but for a horizontal dipole (2 juxtaposed cavities in FEM).
In this case, it can be seen that the vertical component of the dipole deviates more from that simulated numerically.
Figure 8 shows the compared results for a dipole with ratios ɛ = 0.23 and a/δ = 0.22 and the deformation measured by simulating an equivalent FE model.The two extreme cases of a vertical dipole and a horizontal dipole are considered for a broader comparison.It is possible to note that the misfit between the two solutions is low when the assumptions on ɛ and a/δ are met.

Earth and Space Science
10.1029/2023EA003003 CANNAVÓ Two evidences from the discrepancy analysis support the hypothesis of non-interactivity of the pair of sources in a dipole: (a) that assuming the proposed model is correct, where the solutions are quite similar to those of two independent non-interacting analytical sources (Figure 6); (b) that assuming the propose model is correct, where the solutions are similar to those simulated numerically with two independent but interacting sources (Figure 7).

Identifying a Dipole
The deformation pattern at the surface generated by a dipole source is axial symmetric along the dipole axis and degenerates into a central symmetry in a pure vertical dipole (φ = 0).To characterize the dipole source and highlight its peculiarities with respect to the single point source, the predicted deformation fields have been compared.
As mentioned in the previous section, a possible signature of a dipole source can be identified in the changing of sign in the vertical component of the deformation field.In this section, other peculiarities are analyzed to recognize and distinguish a dipole from a single pressure source.Since the deformation pattern at the surface for a tilted dipole (sin θ ≠ 0) is not centrally symmetric (see Figure 2) and thus easily distinguishable from a single point source, only the pure vertical dipole is taken into account in this analysis.
In a purely vertical generic dipole, the horizontal and vertical components of Equation 15 can be written in a simplified form as the following:

Analysis of Vertical Deformation
From Equation 16b, it can be calculated that the above-mentioned signature of the presence of a dipole (and thus a magma transfer) is the sign change of vertical component at a radius of x = 4).The sign change vanishes for κ ∉ ] 2,1[ when one of the two sources becomes predominant.
The maximum vertical deformation is reached at x = 0, while the minimum opposite vertical deformation is and its absolute value is 2 5/ 2 of the maximum value (see Figure 4) which is 0% for κ = 1 (i.e., there is no sign change in the deformation field), it is ≈ 1.8% for κ = 0 and reaches ≈20% with κ = 1.It can be noted as, for small κ, the opposite (in sign) vertical deformation is very small compared to the maximum vertical deformation, making it more difficult to recognize the presence of a vertical dipole.As κ decreases approaching 2 the opposite vertical deformation component becomes the predominant one inducing very atypical deformation patterns.
To compare the predicted vertical surface deformations of a point source and a vertical dipole source for the same volume variation ΔV, two constraints have been imposed: that they have the same maximum value (reached at projection of the source center on the surface), and that they assume the same value at an horizontal distance of d from the source center.The former condition is satisfied for ] otherwise, and in this case the latter condition is satisfied for whose approximate explicit solutions are d d ≈ 1.6d m for κ = 0, d d ≈ 1.4d m for κ = 1, and d d ≈ 2.16d m for κ = 1 where d d and d m are the depths of the dipole source and the point source respectively.The ratio ɛ becomes 2.56, 0.98 and 4.68 for κ 0, 1 and 1 respectively.The high values of ɛ make the mathematical correctness of the found solutions compromised, violating the assumption of small movements in the model formulation.Thus, it is not possible to compare the dipole solutions with the point source ones under the two restrictive constraints.
To overcome this problem, the constraints were relaxed to have equal maximum vertical deformation amplitude with the point source and make the small movements hypothesis hold by setting a small ratio between the dipole axis and depth (i.e.ɛ).
In a vertical dipole the constraint leads to: Equation 18 shows that for a feasible dipole (i.e.ɛ ≪ 1) the source has to be shallower than a point source with the equal maximum vertical deformation for the same volume change ΔV.
Analyzing the limit case to consider the above hypothesis still valid for ɛ = 1/3, the asymmetric (κ = 1) vertical dipole source produces the same maximum vertical deformation as the point source at the same depth and with the same reference ΔV.This is despite the fact that the upper part of the dipole is shallower than the compared point source.However, Figure 9 shows that the vertical deformation decays much quicker with the distance from the source compared with the one produced by a single equivalent point source.As expected, the effect of the deeper part of the source becomes ever more predominant as the asymmetry factor k decreases (e.g., green line in Figure 9).
For a vertical dipole source, the decay of the vertical component of deformation with distance from the source is greater than that of a point source of equal ΔV and equal depth d up to a radius of: Considering the case of equal maximum vertical deformation met for Equation 18, the radius of faster decay simplifies into x < ̅̅̅̅̅̅̅ ̅ 2/3 √ ⋅ d from the source center projection on the surface.
For a more in-depth comparison analysis, the ratio of volume changes between a vertical dipole source and a point source is calculated when both sources produce the same maximum vertical deformation: This means that there is a quadratic relationship between the volume change ratio and the depth ratio for the two kinds of sources (providing they show the same maximum vertical deformation) indicating that, compared to a point source, a dipole source can have a quite different volume change on moving away from the point source depth.

Analysis of Horizontal Deformation
The maximum horizontal displacement is reached at a horizontal distance d/2 from the center for a symmetric dipole (κ = 0), while it is reached at x = d Following the analysis carried out for the vertical component of deformation, to compare the horizontal deformations of the point and a vertical dipole source with the same ΔV, two similar constraints have been imposed: that their maxima are reached at the same distances from the source, and that they have the same maximum value.
The former condition is satisfied for between the maxima occurs for that with the first constraint becomes In this case, the ratio ɛ ≈ 0.53 for κ = 1 thus unfeasible for the hypothesis of small movements.
For κ = 0 the tight conditions are satisfied for ε = 1/ 24 , resulting in ɛ ≈ 0.90 and confirming the incorrectness of meeting both conditions simultaneously.For κ = 1 the ɛ even reaches the value 2.0.Also in this case, the more feasible constraints have been imposed, that is, to show equal maximum horizontal deformation amplitude with the point source and preserve the small movements hypothesis.
For the horizontal displacement, given Equation 21, the relaxed constraints lead to Figure 10 shows how the maximum horizontal deformation (for a ratio ɛ = 1/3) occurs closer to the source projected center compared with the single point source.
It is possible to observe, looking at the horizontal deformations, that the two types of sources can be indistinguishable in the near field, while the difference becomes more evident in the far field where the dipole source shows a smaller deformation.The ratio of the horizontal deformations due to a vertical dipole source and a single point source is given by Equation 23that simplifies into Equation 24 for κ = 0.
In the hypothesis of ɛ ≪ 1, the ratio in Equation 24 is smaller than 1 meaning that a symmetric vertical dipole source generally produces a smaller horizontal deformation compared to an equivalent single point source for the same ΔV.
As for the volume change ratio comparison, given a vertical dipole source and a point source with the same maximum horizontal deformation, the ratio between the two expected volume changes is: Even if in this case the coefficient is more complicated than in Equation 20, but providing they show the same maximum horizontal deformation, the relationship between the volume change ratio and the depth ratio for the two kind of sources is still quadratic.
The dipole source can show a characteristic signature also on tilt data.Following Equation 9, a vertical symmetric dipole changes polarity at distance 2 times the depth apart from the source location projected on the free surface (see Figure 11).
In the same case, the maximum tilt deformation is reached at distance ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ 27

Cases of Study
The existence and usefulness of a dipole source in a volcanic environment were researched for some case studies already known in the literature, where ground deformation data were available.The cases cover the volcanoes Sinabung and Etna.

Sinabung 2013 Eruption
Sinabung volcano is a 2.5km-high andesitic strato-volcano located in northern Sumatra of Indonesia.In mid-September 2013, Sinabung showed phreatic eruptions that led to the main vulcanian eruption occurring on November 23, continuing with lava emission and pyroclastic flows throughout December and January 2014 (Hotta et al., 2019;Pallister et al., 2019).Hotta et al. (2019) used continuous GPS observations at Sinabung to model magma migration and effusion before and during the onset of eruptive activity in September 2013.They divided the whole period into four sub-periods of coherent deformation: two extensional stages followed by two contractional stages.The Sinabung study case considers the deformation field measured by 5 GPS stations.They provide the calculated displacements and their uncertainties at the four GPS stations located in the area.The Period 2 between October and December 2013 represents the inflation stage that brought the magma to the surface and the ensuing lava lake formation in December.This is the period of magma migration toward shallower crust, thus the most likely to fit a dipole source.Indications on a possible dipole source are also given by the vertical components of the displacements: three stations registered an uplift and one a slight lowering that can indeed be the signature of a dipole.To check the capacity of a dipole to infer the measured deformation field, the GPS data has been inverted to find the source parameters able to minimize the Weighted Root Mean Square Error of the predicted deformation field.The weights used in the inversion were the reciprocals of the propagated uncertainties associated to the measurements.In particular, following the approach used by Hotta et al. (2019), the relative displacements between all the station pairs were considered to calculate the error function during the inversion, as shown in Equation 26 where u i,j , v i,j and w i,j represent the relative east-west, north-south and up-down displacements at station i with respect to station j, respectively; ˆrefers to the calculated values; and σ indicates the associated uncertainty.This error function avoids the need to define a reference system for the ground displacement measurements.
The optima source parameters were found by using the Simplicial Homology Global Optimization (SHGO) algorithm (Endres et al., 2018a) with Sobol' sequences sampling (Joe & Kuo, 2008) and are shown in Table 1.The predicted deformations are compared with those measured in Table 2.
The estimated values of dipole parameters indicate a source beneath the volcano compatible in location with the source obtained in Hotta et al. (2019) for the same period, when they found an inflation source at a depth of 0.9 (0.4-2.1) km bsl with a volume change of 0.39 (0.18-0.60)Mm 3 .
To assess how well the parameters are constrained by the data, the confidence intervals have been estimated.Assuming data errors are normally distributed, the approximated confidence intervals have been calculated by using the Ftest approach (Arnadóttir & Segall, 1994) with a confidence level of 90%.For the considered period, the results show a magma transfer dipole oriented almost WE (azimuth = 265°).The estimated zenith angle of the dipole source (41°) approaches 45°, meaning that the horizontal component of the magma motion is significant.This inflation source can be considered to correspond to magma migration toward shallower depths and the summit.It is worth noting that the three-dimensional orientation of the dipole axis is compatible with the distribution of the 3 sources derived by Hotta et al. (2019) to track the magma migration.
The mean GPS station altitude is ∼1,300 m a.s.l., thus the dipole source depth in a half-space is ∼2,100 m.Considering as rule of thumb a maximum CANNAVÓ movement equal to one-third of the depth in order not to violate the hypothesis of δ ≪ d, it is possible to derive the volume change from the estimated hypervolume, ΔV = 303 ⋅ 10 6 /700 = 0.43 ⋅ 10 6 m 3 which matches the estimates of Hotta et al. ( 2019) (Table 2).

Etna 2018 Eruption
Mt. Etna volcano (Italy) is one of the most active volcanoes on Earth with its recurrent explosive and effusive activities.The last activity phase characterized by frequent summit paroxysms started in January 2011, with more than 100 paroxysmal eruptive episodes to date, and 1 significant intrusive episode in December 2018 (Aloisi et al., 2020).The 2018 intrusion suggests a build-up of energy in the preceding period that could imply a movement of magma from the deep Paonita et al. (2021).For this reason, in this case study, the last part of the inflation period leading to the intrusion has been taken into account.Daily solutions of the positions for the GNSS stations belonging to the INGV Etnean permanent network were computed for the period July-December 2018, following the same indications used in Cannavò (2019).A total of 18 stations were considered for their data continuity and steadiness (see Figure 12), and their insensitivity to instability of the eastern flank likely not directly related with volcanic activity (Aloisi et al., 2011;Cannavó et al., 2014).
A sliding time window of 120 days has been used to estimate the ground velocity field.This window enables a coherent inflationary field by smoothing out noise components from the GNSS positioning time series.With this time window, the first useful velocity field for analysis is in October.
The velocity fields and their associated uncertainties have been calculated on a daily basis using the Theil-Sen estimator for a robust linear regression (Sen, 1968;Theil, 1950).
In preliminary inversions for the whole considered period, the inferred dipole source showed an almost constant sub-vertical direction (small φ) at high depths (4 7 km b.s.l.).This may be considered plausible given the greater isotropy at depth.Hence, to reduce the degrees of freedom of the model, in the subsequent inversions, the zenith angle was set to 0. Such a pure vertical dipole implies that also the azimuth parameter becomes redundant and can be considered fixed (e.g. to 0), reducing the degrees of freedom even further.The varying-depth first order topographic correction Williams and Wadge (2000) has been applied to take into account the different altitudes of the GNSS station locations.
With the aim of evaluating the possible presence of a dipole source able to better explain the measured data, this type of source has been compared with the most common source of inflation, the point source (Mogi, 1958;Taylor et al., 2021).For each source type, independent inversions were carried out for all the calculated daily velocity fields to estimate the values of source parameters along the considered time window.A Monte Carlo simulation of 300 inversions was performed for each day, with each having a sample of the velocity field distribution (considering the velocity field distributed as a multivariate Gaussian with the covariance matrix associated to the estimated velocity uncertainty) and random but equal initial conditions for the two types of sources (for the common parameters).Earth and Space Science 10.1029/2023EA003003 An algorithm for optimizing (minimizing) the RMS (root mean square) of the error between measured and predicted deformations was used to invert the data and identify model parameters.The algorithm used is the SHGO (Endres et al., 2018b) which is appropriate for solving general purpose nonlinear and blackbox optimization problems for global optimality.The search domain is constrained within the ranges reported in Table 3.
The simultaneous solutions for the two types of sources have been compared to select the best one.The wellestablished Akaike information criterion (AIC) (Bozdogan, 1987) has been considered to identify the best model.In problems with small data sets the AIC tends to overfit unbalancing on the model with larger parameters (Claeskens & Hjort, 2008).To address such potential problem, the corrected version for small sample sizes of AIC (the AICc) has been used as the criterion to select the most likely source model between dipole source and point source.
For one Monte Carlo simulation (i.e., for each day), the number of times each of the two models was selected as the most plausible source was counted.Figure 13a shows the percentage of times each of the two models resulted as the winner, and therefore more plausible, according to the AICc test.The goodness of fit for the selected models was calculated through the root mean squared error whose development is shown in Figure 14 indicating a temporal improvement in the fitting.The quality of fit is also verified by means of velocity time series analysis as shown in Figure 15, where the measured and predicted deformation values for the vertical component at ESLN station (Figure 12) are compared.Comparison of the value distributions shows that after an initial period of unclear overlap, the models predict the measured values with good statistical accuracy.More figures comparing solutions at different stations can be found in the, Figures S1-S18 in Supporting Information S1.
Figure 13b shows the distributions of depths for the most likely sources.It is worth noting that source depth remains stable within the associated uncertainty throughout the whole period, giving a robust indication on the  If for most of the period a spherical source proves the most plausible, thus indicating a stasis phase, it can be noted that days before the eruption the dipole source becomes statistically predominant, thus indicating that a magmatic ascent phase toward the surface can better explain the measured data.Earth and Space Science 10.1029/2023EA003003 Due to the delay introduced by the window considered for the deformation field estimation, the time evolution of solutions should not be considered synchronized with the real physical processes.It shows only when the processes become detectable with the available measurements and their associated uncertainties.Thus the real magma ascent may have started even before the 5 days to eruption identified by this analysis.
As for the other model parameters, it is necessary to recall that the hypervolume of a dipole source is the product of the volume variation ΔV and the axis length δ. Figure 16 shows the curve of admissible values and the associated uncertainty band for these parameters given the mean hypervolume estimated in the 5 days preceding the eruption.For cases in which the point source has found to be the most likely, the estimated mean volume variation is 2.5 ± 0.5 ⋅ 10 6 m 3 .
The hypothesis of small movements in the dipole formulation indicates that the ratio (ɛ) between dipole axis and depth has to be small.For the near-limit ratio ɛ of 1/3, the estimate volume variation is ≈17.5 ± 5 ⋅ 10 6 m 3 , which is much higher than the volume variation estimated by the spherical source in the period preceding the dipole activation.Smaller ratios implicate even larger volumes (see Figure 16), while a hypothetical longer axis (in this case vertical by hypothesis) reduces the required volume change.Against this variation in volume involved, the total erupted volume was ≈3 ⋅ 10 6 m 3 (Calvari et al., 2020).This discrepancy in volumes has already been highlighted in previous works (Bonforte et al., 2019;Cannavo' et al., 2019).It is also significant to note that the estimated κ parameter for the winning dipole models in the periods close to the eruption takes on slightly negative values at 0.35 ± 0.73, indicating a lower magma volume rise than the volume change estimated in the deepest deflating part of the source.This effect can be likely due to the dissipation of a considerable amount of energy in the process of creating the ascending pathway.Nevertheless, the activity following this eruption involved paroxysmal episodes and silent effusive activity (De Beni et al., 2021), possible signs indicative of resident magmatic presence in the shallow zone stored during the 2018 eruption.

Discussion and Conclusions
The evolution of surface displacement in volcanic zones results from the intricate interplay between the dynamics of pressure sources and involved rheologies.While previous investigations primarily focused on the temporal evolution of displacement linked to magma dynamics centered around one reservoir supplied by a deeper source (irrelevant to deformation) (Edmonds et al., 2019;Reverso et al., 2014;Segall, 2019), recent studies have shown that in some cases the displacement field can be adequately explained only by taking simultaneously into account more than one magma reservoir at distinct depths, all influencing the deformation field (Bagnardi & Amelung, 2012;Chadwick et al., 2010;Elsworth et al., 2008;Foroozan et al., 2010;Hautmann et al., 2013;Melnik & Costa, 2014;Reverso et al., 2014).The prevailing trend in the literature involves proposing complex plumbing systems featuring two-connected-magma chamber models.In these models, the positions of the two chambers are either estimated independently or predefined (e.g., Reverso et al., 2014).This is for example, the case of the recent de Sagazan et al. (2024), in which the authors interpret the latest submarine volcanic activity in Mayotte with two compliant deep magma pressurized magma reservoirs even to explain the local seismicity.
Besides, there are even more complex models that add different physical phenomena in the plumbing system and different rheologies to predict the deformation measured at the surface (e.g., Alshembari et al., 2022;Anderson & Segall, 2011;Kozono, 2021;Marsden et al., 2019).The downside is that these models require knowing many parameters that are often quite difficult to estimate accurately in real applications (Hickey & Gottsmann, 2014); this implies a difficulty in calibrating the models and affects precision and accuracy of their results.In addition, their predictions, taking into account the uncertainties associated on many parameters, often do not deviate much from simpler models (Taylor et al., 2021).In the category of the currently most reliable models are the numerical models, which, however, also suffer from the need for intensive computation, making them challenging to apply for precise inverse modeling.

10.1029/2023EA003003
The prevalence of analytical models in volcanic geodetic literature can be attributed to their perceived simplicity and the fact that these models have been developed, shared, and refined over many years (e.g., Battaglia et al., 2013).Moreover, these simpler models that are quicker to use are often useful for providing first-order estimates of source parameters for urgent answers in near-real time before or during an eruption, or when a volcano that has not been studied before starts to exhibit signs of unrest (e.g., Pritchard & Simons, 2002).
The dipole model presented here originates from the necessity to describe in a single analytical formulation the deformation effects of differential pressure/volume variation in space that can be attributed to magma (or fluids in general) in a volcanic region.A peculiarity, which is also a benefit, of the model is that with a single measurement of the deformation field a possible magma transfer can be estimated.Through data inversion, indeed, the proposed model attempts to quickly answer the question of whether there has been movement of a pressure source (e.g., a magma batch) and if so in which direction and to what magnitude, from the measured deformation data.
Due to differential pressure sources that can mimic magma movement, the dipole model has the advantage of being able to generate a considerable variety of deformation patterns.The spectrum of possible deformation patterns generated by Equation 15 is basically controlled by the parameter κ and the two angles θ and φ.Complex patterns may show asymmetries and changes of polarity in the deformations.Some distinguishing characteristics of a deformation pattern generated by a dipole are: possible change of sign in values on moving away from the source; narrow field of action, that is, high deformation gradient decay on leaving the source.However, in the case of vertical dipole, the deformation patterns generated by a dipole may not be easily distinguishable from those of a single pressure source.In this regard, the best deformation data to better constrain such kind of source should be the space-continuum data provided by InSAR.
A disadvantage of dipole model is the difficulty in separating the spatial dimensions of the source (i.e., the variation in volume and the amplitude of the dipole's axis of motion); in fact, only their product can be known.On the other hand, even the simple point model has shortcomings due, for example, to ambiguity between depth and volume (pressure).Another weakness in the dipole's approximations is to not explicitly account for any deformation induced by the conduit dynamics connecting the two ends of the dipole.Nevertheless, the formulation explicitly requires that the dipole axis (i.e., conduit) extent be negligible with respect to the size of the domain.This restrictive assumption mitigates possible conduit effects neglected in the model.
Comparison analysis with finite element models showed the range of applicability of the analytical approximation of the dipole model, confirming the good agreement of the predicted deformations under the small movements hypothesis (i.e., small dipole axis with respect to depth).Despite the possible uncertainties in identifying a dipole source using sparse GNSS measurements, a procedure was developed to statistically verify the effectiveness of this type of source.The procedure was applied to a known case study that has already been extensively studied in the literature such as the December 2018 eruption at Mt. Etna.The model selection strategy plays a central role in the quantitative evaluation of results, so it requires to be defined in a way that is as robust as possible.In this view, the procedure is designed to mitigate statistical fluctuations and consists of comparing by means of Monte Carlo simulation the dipole source against another type of source and choosing the most plausible one according to an objective criterion such as AIC (Burnham et al., 1998).
The dipole model has enabled the interpretation of the displacement measured by GNSS at Mt. Etna in the period preceding the 2018 eruption.This case study has revealed that the dipole source can be useful in explaining some deformation patterns better than a simple single inflationary or deflationary source.In particular, the analyzed case shows that the dipole is activated in advance just as the eruption is approaching, suggesting a shift from a static to a dynamic phase of unrest.The analysis indicates not only the ability of a dipole model to predict a real deformation field, but also its unique quality to describe magma motion (e.g., preceding the intrusion of the surface dyke that has led to the eruption).This ability to identify an activation phase of magma movements can also be decisive from a monitoring point of view, and in terms of risk mitigation (Rosi et al., 2022;Sparks & Cashman, 2017).

Earth and Space Science
10.1029/2023EA003003 CANNAVÓ An interesting development of this model could concern the introduction of the mass component and the calculation of the predicted gravity field variation.This could be useful in quantifying the masses involved and better characterizing the cause of the unrest as for example, in Battaglia et al. (2003) and Nikkhoo and Rivalta (2022).In future works, it may also be useful to compare the solutions of complex numerical dynamical models provided by newly published papers that have further explored cases of binary sources with the dipole analytical solutions presented here.For example, recently, Furst et al. (2023) developed a viscous model based on Weertman (1971) crack concept, which describes a magmatic intrusion that has many similarities with a dipole (tensile stress/compression) approach.

Figure 2 .
Figure 2. Deformation field due to a dipole source for different values of elevation angles.θ is the azimuth angle (clockwise from North) of the dipole axis directed from the deflation source to the inflation one and φ is the elevation angle measured from the vertical.Color is associated to the vertical displacement.Distances are normalized by the source depth d.The asymmetry of color scale is necessary to make the vertical negative component visible.The black arrows show the horizontal pattern of the deformation field.

Figure 3 .
Figure3.The line in black shows the absolute distance between the points with extreme vertical deformations (in red the position of maximum, in blue the position of the minimum).It can be seen that the maximum distance is reached for a pure vertical dipole and it equals twice the depth location, while the minimum distance is reached for a perfectly horizontal dipole (i.e., elevation 90°) and it equals the depth location.In green the location of the zero vertical deformation.

Figure 4 .
Figure 4. Profiles of axisymmetric displacements (vertical (red), horizontal (blue), magnitude of the total displacement (black)) for a vertical symmetric dipole source (κ = 0) as a function of horizontal distance in source depths.The displacements are scaled by 2ΔVδ 1 ν πd 3 that is the maximum vertical deformation at x = 0 (thus the y-axis is dimensionless).The dashed vertical lines mark the maximum horizontal displacement x = d/2, the distance ( ̅̅ ̅ 2 √ d) where vertical displacement changes sign, and the distance (|x| = 2d) where the vertical deformation assumes the maximum opposite value.

Figure 5 .
Figure 5.Comparison of the deformation components predicted by a dipole and a two-point source.The source parameters are: d = 2,000 m, δ = 500 m and ΔV = 10 * 10 6 m 3 .The deformations are shown as a function of the horizontal distance normalized by the source depth.In the insets, the schematic position of the sources.

Figure 6 .
Figure 6.Maximum relative discrepancies of predicted deformation between a dipole model and two point-sources superposition model for source positions as shown in the insets of Figure 5.

Figure 7 .
Figure 7. Relative discrepancies of expected deformations between a dipole model and a finite element model model with two spherical-point sources.The color of dots indicates the ratio between the radius of the spherical pole and the dipole axis length.In the insets, the schematic position of the dipole and the direction of the deformation component considered.

Figure 8 .
Figure 8.Comparison between solutions of dipole models and the equivalent finite element model models for parameters ɛ = 0.23 and a/δ = 0.22.

Figure 9 .
Figure 9.Comparison of vertical predicted deformation by a point source (in black), a feasible (i.e., ɛ = 1/3) vertical symmetric (κ = 0) dipole source (in red), a feasible asymmetric (κ = 1) vertical dipole source (in blue) and a feasible asymmetric (κ = 0.5) vertical dipole source (in green).The feasible dipoles are constrained to have the same maximum vertical deformation and the same volume variation as the point source.The deformations are scaled by the maximum value and are shown as a function of the horizontal distance normalized by the source depth for point source d m .
vertical dipole.For κ ∈ [ 1,1] the maximum horizontal displacement is reached between ≈0.43 and ≈0.54 times the depth d, thus κ has little influence on the spatial scope of the vertical dipole.

Figure 10 .
Figure 10.Comparison of horizontal deformation predicted by a point source (in black) and feasible (i.e., ɛ = 1/3) vertical dipole sources with κ = 0 (in red), κ = 1 (in blue) and κ = 0.5 (in green).The feasible dipoles are constrained to have the same maximum horizontal deformation and the same volume variation as the point source.The deformations are scaled by the maximum value and are shown as function of the horizontal distance scaled by the source depth for point source d m .

Figure 11 .
Figure 11.Profile of axisymmetric tilt deformation for a vertical symmetric dipole source as a function of horizontal distance in source depth.The tilt deformation is scaled by the maximum tilt deformation.The vertical dashed lines represent the maximum/minimum, the polarity changes and the maximum/minimum with opposite sign.

Figure 13 .
Figure 13.Results from Monte Carlo simulations to assess the likely acting source during the considered period for Mt.Etna.(a) Fraction of success according to the AICc criterion for the two models.(b) Violin plot of distributions (along the vertical axis) of depth locations for the most likely sources.

Figure 14 .
Figure 14.Violin plot of distributions (along the vertical axis) of Root Mean Squared Error of the fits for the selected best models in the Monte Carlo simulations.The vertical dashed red line indicates the time when the dipole source becomes statistically predominant.

Figure 15 .
Figure 15.Violin plot of distributions (in vertical) of model fits for the vertical component of ESLN station.In gray, the measured data with the associated uncertainty distributions; in red the predicted values for the selected best models.The vertical dashed red line indicates the time when the dipole source becomes statistically predominant.

Figure 16 .
Figure16.Curve of relationship between volume variation (y-axis) and dipole axis length (expressed as fraction of source depth, x-axis) for the estimated dipole hyper-volume in the last days before the eruption.The light-blue area represents the associated uncertainty band.

Table 1
Estimated Source Parameters for a Dipole in a Half-Space (Equation15)

Table 2
Deformation Field Measured and Predicted Figure 12.GNSS permanent stations with their names at Mt Etna used in the analysis.

Table 3
Constraints of Model Parameters Used for Inversions