Multiridge Method for Studying Ground-Deformation Sources: Application to Volcanic Environments

Volcanic phenomena are currently monitored by the detection of physical and chemical observations. Generally, the ground deformation field is the most relevant shallow expression of the geometric and physical parameters variations in the magmatic reservoir. In this study, we propose a novel method for the direct estimation of the geometric parameters of sources responsible for volcanic ground deformation detected via the DInSAR technique. Starting with the biharmonic properties of the deformation field, we define an approach based on the Multiridge and ScalFun methods to achieve relevant information about both the positions and shapes of active sources, such as the Mogi source. Our methodology is definitely different from the methods currently used for modeling ground-deformation sources, mainly based on forward or inverse techniques. In fact, (i) it does not require any assumptions about the source type, and (ii) it is not influenced by the distribution of medium elastic parameters or (iii) the presence of high-frequency noise in the dataset. For synthetic cases, we accurately estimate the depth to the source within a 3% error. Finally, we study the real case of the Okmok volcano ground-deformation field and achieve results compatible with those in previous works.

In recent decades, the development of remote sensing technologies has provided relevant efforts in the context of earth science. Specifically, the Differential SAR Interferometry (DInSAR) technique has significantly contributed to the characterization and imaging of the deformation of Earth's surface, allowing for the analysis of a large variety of geological phenomena, such as landslides 1,2 , earthquakes 3,4 and volcanic events 5,6 .
In a volcanic environment, the accumulation processes of shallow magma in the crust often cause slight movements at the surface, which can be measured using standard land-surveying techniques (leveling) or satellite geodesy methodologies (e.g., GPS and DInSAR). In this context, these aforementioned processes often result from variations in several geometric and physical parameters of magmatic sources (e.g., shape, size, depth, dip, and magnitude of pressure change ΔP) 7 . The estimation of these parameters is performed using two different approaches: (i) forward modeling 8,9 , which is based on a trial-error procedure and strongly depends on the experience of the scientist, and (ii) inverse modeling [10][11][12] , which is faster than the previous approach. Current methods of inversion are essentially based on optimization techniques, where the final solution is that obtaining the best-fit between the measured and modeled data; this optimization/inversion procedure is performed by using algorithms such as Monte-Carlo, Gibbs, Simulated Annealing, and Levenberg-Marquardt 10,13,14 . In the case of the optimization/inversion of geodetic data, several authors have evaluated model parameters via posterior probability density distributions 15,16 . However, we remark that these methods are often unable to provide unique solution 17 .
Starting with these considerations, we propose a novel method based on the theory of Laplace's equation: where ∇ 2 represents the Laplace operator and V represents the potential (or any of its any-order derivatives). The solutions to Laplace's equation, which is the simplest example of elliptic partial differential equations, are the so-called harmonic functions, and the theory is often referred to as the Potential Field Theory 18 (PFT). We study the conditions under which the theory of the deformation field reduces to the PFT and then show how to obtain the main geometric parameters of the active source.
The PFT can be involved to analyze the deformation fields by considering Love's argumentation 19 . We further assume that modeling the magmatic reservoir is performed through the approximation of the elasticity theory 7 and under the assumption that an active volcanic source is not affected by shape changes 19,20 . Among the sources of deformation, we first consider those that also satisfy the homogeneity equation 18 , such as a Mogi source 21 . In this framework, we use methods that have been applied to analyze potential fields 22 , namely, the Multiridge 23 and ScalFun 24 methods. By using these, the geometric features of the active source, such as the position and source type, can be easily detected; the main advantages of the proposed methods are that they are determined independently from (i) other physical features of the source (such as the pressure variation; ∆P), (ii) physical-elastic parameters of the medium (such as the shear modulus, µ and the Poisson ratio ν), and (iii) a low signal-to-noise ratio.
To show the applicability of the proposed methodology, we first carry out a preliminary analysis on the harmonicity and homogeneity properties of the deformation field, based on the simplest volcanic source. Then, we study several test configurations by analyzing the field components and the satellite Line Of Sight (LOS) projection of the deformation field. Finally, we consider the deformation field of the Okmok volcano (Aleutian Islands), which has already been interpreted with standard methods in other works 12,25-28 .

Methods: Potential Theory for the Deformation Field
The method is applied to interpret the large number of deformation points achieved by the SBAS-DInSAR technique 29 . In particular, we consider the Multiridge 23 and ScalFun 24 methods, which are two well-assessed methods, to analyze homogeneous potential fields. For the first approximation, we assume that the ground deformation field results only from volume changes in subvolcanic bodies (i.e., magma chambers) in a linear elastic medium.
Harmonic functions. A harmonic function (V ) can be defined as a function that satisfies Laplace's equation 18 . In particular, a potential field (e.g., x y z F( , , )) is a vector field where the scalar function V (i.e., the potential function) exists 18 : where the potential function V is a harmonic function that satisfies Laplace's equation with its partial derivatives of any order at any point not occupied by the sources of field F: The interpretation of the ground deformation field via the PFT-based methods derives from Love's work on the elasticity theory 19 . According to his arguments, when the modeled deformation field is not caused by a source shape-change, the deformation field u can be represented by the coefficients of the partial differential equation (PDE), with respect to the coordinates of a single scalar function (φ) 19 : By combining equation (4) with Navier's equation, to describe the features of u in a homogeneous and isotropic medium under equilibrium conditions: 2 we find that the displacement potential φ is a biharmonic function since φ satisfies 20 : 2 2 4 This result may be also obtained under the Helmholtz formulation 20 . We show that the ground deformation field modeled with the Mogi source also encompasses harmonic properties. The deformation field modeled by a hydrostatic pressure change within a spherical cavity embedded in an elastic half-space, with a radius much smaller than its depth 21 , is given by: It is simple to argue that u is the gradient of the Newtonian potential in the form 1/r, which is a harmonic function 18 . Therefore, all of the components of the deformation field u are also harmonic, and we find that: In the case of the LOS satellite data analysis, the projection of the modeled deformation (u LOS ) can be simply calculated by combining the ground deformation field components (u v w , , ) with the LOS unit vector as follows: where c x , c y , and c z are the components of the LOS vector c. The operation of equation (9) also returns a harmonic function 18 , since we consider the mean values of LOS vector components.
Multiridge method. The Multiridge 23 method is a multiscale method based on the analysis of so-called ridges, which are defined as lines passing through the maxima of a field and its derivatives at different scales. We emphasize that this method can only be applied in cases when the field can be expressed by harmonic functions. We consider a coordinate system where the xand y-directions are represented by the North-South and East-West directions, respectively, and the z-axis is definite and positive downward. If P x y z ( , , ) is a generic field generated by a simple source Q x y z ( , , ) 0 0 0 and by taking into account the cross section = y y 0 , it is possible to obtain three equations for the ridges relative to the minima of the horizontal derivative of the field 23 : where γ is the angular coefficient of the straight line γ β = tan , where β represents the angle of the ridge line from the vertical z-axis. Since the three ridges intersect at the source center x z ( , ) 0 0 , its position can be simply individuated with this geometric method.
Specifically, the Multiridge method mainly consists of three phases: (1) The creation of a multiscale dataset by performing upward continuation 18,30 from the original measurement level to different levels, (2) Individuation and representation of the edges, and (3) Representation and continuation of the ridges down to the source-region, to individuate the correct depth of the source at the intersection of more ridges.
Regarding the application of the upward continuation operator to the deformation field (1), we specify that we are clearly not considering that the deformation field could propagate into the air. The upward-continued deformation field instead corresponds to that which would have been produced by the same source in a region that was upper extended by the same amount of upward continuation. For example, if the analyzed field is observed by a source with a depth of = z 2 0 km, the upward-continued field towards a 1-km altitude is generated by the same source but with a center depth at = z 3 0 km. The scientific soundness of this procedure is shown in the Supplementary Information ( Figure S1).
For the concerns of this study (2), we clarify that more than one subset can be defined: a multiridge subset, R 1 , that is individuated by the zeros of the horizontal derivatives; a multiridge subset, R 2 , that is individuated by the zeros of the first-order vertical derivative; and a multiridge subset, R 3 , that is individuated by the maxima of the field 22 . In our analysis, we have selected the ridges of subsets R 1 and R 2 . Moreover, since the method involves the field at different levels (i.e., multiscale), high-frequency noise can be easily recognized in the multiridge subsets. For this reason, in the real-case data analysis, we have to exclude the edges related to noise.
For the representation and continuation of the ridges (3), we specify that each ridge is determined by a best-fit linear regression within a 95% confidence interval; in particular, we calculate the determination coefficient R 2 , where R is the correlation coefficient, which represents a statistical measure of how the data (edges) are close to the fitted regression line (ridges). Moreover, we evaluate the solution uncertainties (intersection at the ridges) by considering the error on the best-fit linear regression coefficients (intercept and slope constants). We remark that in the case of irregular (not flat) data measurement levels, we have to perform another procedure before applying the Multiridge method. Since the method is based on a level-to-level algorithm, and the ground deformation field is measured at the topographic surface (which is generally detailed), the dataset must be processed to numerically generate a ground deformation field that could have been measured in the case of constant distance between the source and measurement surface. This transformation relocates the analyzed field onto the constant reduction level, which is performed in this study by using the CWT-domain algorithm 31 . The accuracies and details of this procedure are shown in the Supplementary Information (Figures S2 and S3).

Homogeneous functions.
Consider again the ground deformation field modeled through the Mogi source, and let us evaluate its homogeneity properties.
The function U is said to be homogeneous, with degree n, if satisfies the Euler's equation 18 : The importance of the homogeneous functions derives from the fact that simple source fields are homogeneous functions of degree n which, for many sources, reflects the falloff rate of the potential field anomaly with distance 24 . By inserting equation (11) into equation (7), it is obvious that each component of u for Mogi's source model is homogeneous of degree = − n 2. The homogeneity degree of the field (n) may be used to estimate the homogeneity degree of the source (n s ) 32 which, for a Mogi model, is given by: Since n s is a source parameter determined by a field parameter (n), it may be convenient to refer to its opposite, the so-called Structural Index (N ): s which is a parameter utilized within the framework of homogeneous functions to study gravity and magnetic field anomalies 33,34 . We conclude that the Mogi model source is characterized by =− n 3 s and = .
N 3 ScalFun method. The ScalFun method is based on the properties of the scaling function, which was introduced into the framework of the DEXP theory 24 to estimate the homogeneity degree of the observed field (n); this is, in turn, related to the homogeneous property of the source (n s ). For any p th vertical derivative for the Newtonian potential of a pole source f z ( ) p at = x x 0 and = y y 0 , we define the scaling function τ p of order p as: p p 0 where = − + n p ( 1) represents the degree of homogeneity of f p . Equation (14) can be rewritten by = z q 1 ; therefore, τ p becomes a function of q: p 0 which means that τ → = . q n ( 0) p Hence, in a plot diagram of τ p as function of q, the intercept gives an estimate of the homogeneity degree n. Starting from the z 0 source depth retrieved by using the Multiridge method, we can use equation (15) to estimate n, n s and N; these values give us information about the geometry of the source.

Result I: Application of the Potential Theory to the Synthetic Deformation Field
To show the validity of the proposed method, we perform several synthetic tests based on the analysis of the deformation field generated by a spherical source. We chose Cartesian reference system, with the origin of the system located at the point O(0, 0); the xand y-axes are oriented in the E-W and N-S directions, respectively, and the z-axis is negative downward. For all of the following tests, the deformation field is simulated with a grid of 120 km x 120 km sampled at 0.1 km intervals.

Analysis of the vertical component of the deformation field.
We generate the vertical component of the ground deformation field on a flat surface ( = z 0) by using the Mogi model, as produced by the overpressured (Δ = P 10 MPa) spherical magma chamber (Fig. 1a). The active source depth is −2 km and is located 60 km along both the xand y-directions; its radius is 0.3 km, while the medium shear modulus and Poisson's coefficient are 1 GPa and 0.25 [-], respectively. We apply the Multiridge method to the vertical component, analyze the AA' E-W profile passing for the maximum value of the field (Fig. 1b). The results allow us to easily identify a source at an approximate −2.05 ± 0.01 km depth and located at 60 ± 0.02 km along the x-direction (Fig. 1c). The estimation reliability is supported by the high values of R 2 . Then, we apply the ScalFun method to the central ridge (cyan) and achieve a homogeneity degree of ≈ − n 2 (Fig. 1d). This value corresponds to the Structural Index value of = N 3, suggesting that the representative geometry of the source is spherical 32 , which is in agreement with the features of the true source.
Subsequently, we consider the previous test adding to the modeled vertical component a 10% of high-frequency noise (with respect to the maximum value of the anomaly) (Fig. 2a,b). The results show that the estimates are not conditioned by the lower signal-to-noise ratio since (also in this case) we identify a source depth of approximately −2.04 ± 0.02 km, and a location of 60 ± 0.03 km along the x-direction (Fig. 2c). Then, we apply the ScalFun method to the central ridge (cyan) in order to evaluate the homogeneity degree of the field, ≈ − n 2 ( = N 3) which, in turn, gives us an evaluation of the geometric shape of the source (Fig. 2d).
For the third case study, we consider the same source model parameters as those of the previous tests, but it is embedded in a heterogeneous elastic space. In particular, using a 3D finite element approach, we run the model by considering the heterogeneities of the elastic parameters (increasing shear modulus with depth) of the medium ( Figure S4). The modeled vertical component is reported in Fig. 3a. By analyzing the AA' profile (Fig. 3b), we retrieve a source depth of approximately −2.03 ± 0.04 km, and a location of 60 ± 0.05 km along the x-direction (Fig. 3c). The application of the ScalFun method confirms the expected result for a Structural Index equal to ≈ N 3 (Fig. 3d). Analysis of the LOS deformation field. We further apply our methodology to the modeled LOS deformation field, which is calculated first on a flat surface at = z 0 (Fig. 4); then, the Okmok topographic surface is considered (Fig. 5) and, finally, the high-frequency noise is added (Fig. 6). All of the synthetic fields are projected along the ascending and descending orbits of the LOS satellite sensor, and we consider a look-angle of 23°.
We consider the same source model parameters as those of the first test, changing only the horizontal position beneath the caldera floor: 58 km in the x-direction and 53 km in the y-direction. By combining the simulated components of the ground deformation field with the LOS unit vectors ([−0.346, −0.081, 0.935]: mean values for ascending orbit; [0.346, −0.081, 0.935]: mean values for descending orbit), we obtain the projected field along the ascending (Fig. 4a-c) and descending (Fig. 4b-d) orbits. In both cases, the retrieved source depth is approximately −2.05 ± 0.02 km, and the location equal to 58.02 ± 0.02 km along the x-direction (Fig. 4e,f). The application of the ScalFun method to both the ascending and descending cases confirms the expected value of the Structural Index ( = N 3) (Fig. 4g,h). Then, we consider the previous test configuration projected along the descending LOS by changing only the surface measurements from = z 0 to z equal to the Okmok volcano topography. We process the simulated dataset to relocate the analyzed field at the new constant reduction level, which is chosen to be 1.5-km a.s.l. After this transformation, the expected depth is equal to −3.5 km (i.e., the distance between the constant reduction level and source depth). We apply the Multiridge and ScalFun methods to the processed field shown in Fig. 5a by considering the AA' profile (Fig. 5b): the first method allows us to identify the source at an approximate −3.53 ± 0.02 km depth (Fig. 5c), and the horizontal position at 58.01 ± 0.02 km along the x-direction; the second method reveals the homogeneity degree of the field, = − n 2 (Fig. 5d). This value corresponds to a source with a Structural Index of = N 3, suggesting that the source geometry is well represented by a spherical geometry. Finally, we consider the aforementioned test configuration, with 10% high-frequency noise (with respect to the maximum value), on the descending LOS deformation field (Fig. 6a,b). The achieved results confirm that the estimated geometric parameters of the source (−3.55 ± 0.04 km for depth and 58.01 ± 0.02 km for horizontal position) are not influenced by the presence of high-frequency noise in the dataset (Fig. 6c,d).
We point out that for all of the performed tests summarized in Table 1, the same parameter uncertainties have also been achieved by analyzing the N-S profile.

Result II: Application of the Potential Theory to the Real Deformation Field
Geological setting of the Okmok volcano. The Okmok volcano is an active caldera field located on the oceanic crust of the central Aleutian arc (Alaska -USA) 35 that represents the surface expression for the subduction of the Pacific Plate as it moves northwards beneath the North American Plate 36 .
In particular, it is a dominantly a basaltic shield volcano covering most of the northeastern end of Umnak Island in Alaska 26 . The Okmok physiography is dominated by a central caldera, with a diameter of 10 km; the rim and caldera floor have elevations of approximately 900 and 400 m a.s.l., respectively 28 . This physiography is the result of two different and large (≈15 km 3 ) caldera-forming events 35 caused by catastrophic pyroclastic eruptions that occurred approximately 12.0 and 2.05 kyr ago, respectively 35,37 . These eruptions began with a small Plinian rhyodacite event, followed by the emplacement of a dominantly andesitic ash-flow unit, whereas effusive interand post-caldera lavas have been more basaltic 35 .
Subsequent eruptions produced a field of small cones, which was once entirely filled with water, and lava flows, including several historically active vents. Observed pillow lavas and other textures are consistent with the subaerial eruptions within the caldera 38 .
For the past 200 years, the Okmok volcano has had an effusive and basaltic eruption every 10-20 years, generally from the intracaldera cones 38 . The most recent eruption, in 2008, originated from several new vents and occurred near the eastern rim of the caldera 28 , while the three previous eruptions (1945, 1958 and 1997) originated near the southwest rim of the caldera 39 . In particular, the 1997 eruption began with steam and ash plumes and progressed into moderate Strombolian activity, producing explosive ash plumes and lava that flowed toward the center of the caldera. Geochemical analyses of the erupted products are consistent with the primitive magma from the depth and brief storage of the shallow reservoirs 35 .

Analysis of the DInSAR measurements.
We analyze the Okmok volcano ground deformation pattern retrieved by processing the ENVISAT SAR images. The interferogram is related to the period 15 July 2003-29 June 2004, with images acquired by the ENVISAT satellite (ESA) along the descending orbit (Track 115). The ENVISAT dataset was processed by using the online P-SBAS web tool available within ESA's Grid Processing On Demand (G-POD) environment, which is within the framework of ESA's Geohazards Exploitation Platform (GEP) 40,41 . The P-SBAS results from the ENVISAT data were spatially averaged (i.e., multilooked) to obtain a pixel size of approximately 80 m by 80 m on the ground. The ENVISAT satellite look angle is 23°, and the mean values of the LOS unit vector are [0.346, −0.081, 0.935] relative to the satellite descending orbit.
After the processing steps, the descending interferogram is unwrapped to retrieve the LOS deformation measurements (Fig. 7a). In the considered period, the DInSAR measurements shown an uplift phenomenon, with a maximum deformation value of approximately 12 cm. Since the measurement surface is detailed (i.e., the Okmok volcano topography), the LOS deformation dataset is processed to obtain the ground deformation field evaluated at the constant reduction level, specifically at 1.5-km a.s.l. (Fig. 7b).
At this stage, we applied our methodology to the DInSAR measurements to evaluate the geometric parameters of the active source; we applied both the Multiridge and ScalFun methods along the E-W (Fig. 8a) and N-S (Fig. 8b) profiles to retrieve the following results: (1) Multiridge method: A source was located at a −4.9 ± 0.06 km depth (Fig. 8c,d) from the constant reduction level (1.5-km a.s.l.) along the z-axis, which corresponded to a depth of −3.4 ± 0.06 km b.s.l. with horizontal UTM coordinates equal to 690.9 ± 0.08 km East (Fig. 8c) and 5924 ± 0.07 km North (Fig. 8d); (2) ScalFun: A homogeneity degree of ≈ − n 2 (Fig. 8e,f) was computed on the central cyan ridge for both the E-W and N-S profiles; this value corresponds to a Structural Index of ≈ N 3, suggesting that the source geometry related to the measured ground deformation field can be well-approximated by a spherical reservoir.

Discussion and Conclusion
In this work, we have proposed and validated a novel method based on the PFT to study the volcanic environment and estimate the geometric parameters of the source model responsible for the ground deformation field.
Multiridge and ScalFun methods provide information on depth, horizontal position and source shape without any a priori information; it is also simple, fast and not computationally expensive. This feature makes it very different from current methods, which are usually based on optimization/inversion procedures and that, most of times, assuming a source-model.
First, we study the harmonicity and homogeneity properties of the theoretical deformation field generated by the Mogi model. We show relevant results since the model equations are simply harmonic and homogeneous functions, with a homogeneous degree of the field of = − n 2, which corresponds to a Structural Index value of = N 3. These findings are valid for all deformation components (both horizontal and vertical) and LOS projections of the ground deformation field; therefore, this methodology may provide information increasing the reliability of the results. Although this study is limited to the analysis of a spherical source, it shows that when retrieving for = N 3, the only source responsible for the ground deformation field is a spherical source. In the future, other more complex geometries of deformation sources will be analyzed.
Then, to show the applicability of the Multiridge and ScalFun methods, we have performed several synthetic tests (Table 1) by considering different model configurations. The geometric parameter results show a very good performance of the methodology, with an uncertainty in the estimated source position of approximately 3% respect to the expected value. In all cases, the ScalFun method retrieves the expected source type by estimating a Structural Index of = N 3. The reliability of the solutions is supported by the high values of R 2 calculated by the linear regression of each ridge. The achieved uncertainties (Table 1) increase with the complexity of the model setting. In particular, the best solutions are retrieved when the vertical component of the deformation field is simulated on a flat surface, while the worst solution when the deformation field is LOS projected and evaluated on real topography with high-frequency noise.
All of the achieved results emphasize that the proposed methods provide solutions that do not depend on the following: different model pressure changes, ΔP; the physical-elastic features of the medium in which the source is embedded (e.g., the shear modulus, µ, and Poisson ratio, ν); and the presence of high-frequency noise in the analyzed dataset. On the other hand, the goodness of fit of the results is influenced by the sampling of deformation anomalies (such as the density of the measurement points) and the entirety of their measurements.
Finally, we apply the proposed methodology to the real case of the Okmok volcano ground deformation pattern retrieved by the DInSAR analysis. In particular, we use the interferogram related to the period 15 July 2003 -29 June 2004, with images acquired by the ENVISAT satellite along the descending orbit. The Okmok volcano deformation field has been studied by many authors 12,[26][27][28] , which all interpreted that deformation was due to the inflation or deflation of a spherical magma chamber. In most of these studies, the source model type (i.e., Mogi model) was a priori assumed, and the elastic parameters (i.e., shear modulus and Poisson ratio) were fixed before inverting the data. Only in one case 28 did the authors use seismic tomography as a priori information to set the heterogeneous distribution of the elastic parameters. For these authors, the source depth ranged from 3.1 to 3.5 km, while the source position ranged from 690.3 km to 690.72 km for the East UTM coordinate and 5923.6 km to 5923.98 km for the North UTM coordinate ( Table 2).
Our results are in good agreement with those of the aforementioned works, indicating a depth source equal to 3.4 ± 0.06 km and a horizontal position at 690.9 ± 0.08 km E and 5924 ± 0.07 km N (Table 2). Moreover, we estimated a Structural Index of = N 3 and, therefore, we can state that the geodetic source geometry is well represented by the spherical model. Additionally, for the Okmok real case, we calculated the R 2 of the best-fit linear regression for each ridge to support the robustness of the proposed solutions. These values are not as high as those retrieved for the synthetic tests, which is probably due to the resolution and sampling step of the measurements. This limitation may be overcome with the measurements retrieved by the new generation of satellite images (e.g., Sentinel-1A and -1B).
In conclusion, since our method does not depend on physical model parameters, we need to apply a subsequent optimization/inversion procedure to fully interpret the ground deformation measurements and determine the source pressure variation, ∆P, and potentially the elastic parameters (such as the shear modulus, µ, and Poisson ratio, ν) of the medium. In the meantime, we stress that the source parameters estimated with the PFT procedure, depth to the center, horizontal position and source-type are decisive to constrain the entire interpretation procedure.  Table 2. Source locations for the Okmok magma chamber retrieved from the DInSAR measurements. *Depths are below sea level.