A scalar function ϕ is harmonic if it satisfies Laplace’s equation (Blakely, 1996), ∇ 2 ϕ   = ∂ 2 ϕ   ∂ x 2 + ∂ 2 ϕ   ∂ y 2 + ∂ 2 ϕ   ∂ z 2 = 0 , (1) and the deformation field u ( x , y , z ) = ( u , v , w ) is also a potential field if any biharmonic component satisfies the following relation (Castaldo et al., 2018): ∇ 2 u ( x , y , z ) = ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 ∂ 2 v ∂ x 2 + ∂ 2 v ∂ y 2 + ∂ 2 v ∂ z 2 ∂ 2 w ∂ x 2 + ∂ 2 w ∂ y 2 + ∂ 2 w ∂ z 2 ) = 0 . (2)

To understand for which conditions Eq. 2 is satisfied, we consider Navier’s relation for equilibrium under surface traction (Love, 1906), ( μ ∇ 2 u + ( λ + μ ) ∂ ∂ x ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) μ ∇ 2 v + ( λ + μ ) ∂ ∂ y ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) μ ∇ 2 w + ( λ + μ ) ∂ ∂ z ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) ) = 0 , (3) where, λ and μ are the Lamè’s constants. We arrange the terms between the left- and right-hand members, μ ∇ 2 u = ( − ( λ + μ ) ∂ ∂ x ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) − ( λ + μ ) ∂ ∂ y ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) − ( λ + μ ) ∂ ∂ z ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) ) , (4) and we divide both the members by μ to get: ∇ 2 u = ( − ( 1 + λ μ ) ∂ ∂ x ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) − ( 1 + λ μ ) ∂ ∂ y ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) − ( 1 + λ μ ) ∂ ∂ z ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) ) . (5)

Therefore, u satisfies Laplace’s equation if: ( ( 1 + λ μ ) ∂ ∂ x ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) ( 1 + λ μ ) ∂ ∂ y ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) ( 1 + λ μ ) ∂ ∂ z ( ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ) ) = 0 . (6)

At this stage, we consider the strain-stress relations expressed by Hooke’s Law (Sadd, 2005), ∂ u ∂ x = 1 E ( σ x x − ν σ y y − ν σ z z ) ∂ v ∂ y = 1 E ( σ y y − ν σ z z − ν σ x x ) ∂ w ∂ z = 1 E ( σ z z − ν σ x x − ν σ y y ) , (7) and, by substituting them into Eq. 6, ( ( 1 + λ μ ) ∂ ∂ x 1 E ( σ x x − ν σ y y − ν σ z z + σ y y − ν σ z z − ν σ x x + σ z z − ν σ x x − ν σ y y ) ( 1 + λ μ ) ∂ ∂ y 1 E ( σ x x − ν σ y y − ν σ z z + σ y y − ν σ z z − ν σ x x + σ z z − ν σ x x − ν σ y y ) ( 1 + λ μ ) ∂ ∂ z 1 E ( σ x x − ν σ y y − ν σ z z + σ y y − ν σ z z − ν σ x x + σ z z − ν σ x x − ν σ y y ) ) = 0 , (8) we get: ( 1 + λ μ ) ∇ ( σ x x E + σ y y E + σ z z E ) − 2 ∇ ( ν σ x x E + ν σ y y E + ν σ z z E ) = 0 , (9) that is 1 1 − 2 ν ∇ ( σ k k E ) − 2 ∇ ( ν σ k k E ) = 0 , (10) where, E = μ ( 3 λ + 2 μ ) λ + μ is Young’s modulus, ν = λ 2 ( λ + μ ) is Poisson’s coefficient, and σ k k = σ x x + σ y y + σ z z . Eq. 10 is verified if ν is constant and ∇ ( σ k k E ) = 0 . We note that σ k k already depends on E and ν (Sadd, 2005), therefore, the deformation field satisfies Laplace’s Eq. 2 when: ∇ σ k k = 0 . (11)

In the framework of deforming volcanic bodies, Eq. 11 is always satisfied for a hydrostatic pressure-change within the source embedded in a homogeneous elastic half-space.

The harmonic functions enjoy two crucial field transformations in the wavenumber domain: the upward continuation and the vertical differentiation (Blakely, 1996). The upward continuation is used to form the multi-scale deformation 3D spatial dataset starting from the 2D spatial measured data (i.e., DInSAR measurements). The upward continued field f c is just the deformation field that would have been generated by the same source of the measured data f and upwardly extending the medium by the amount of continuation Δ z i (Castaldo et al., 2018). Specifically, we retrieve f c by antitransforming into the space-domain the following relation (Blakely, 1996): ℱ [ f c   ] = ℱ [ f ] e − Δ z i | k | , (12) where, ℱ [ f c   ] and ℱ [ f ] are the Fourier transforms of f c and f , respectively, and k is the wavenumber.

The vertical differentiation of order p of the ground deformation field ∂ p f / ∂ z p enhances the high-wavenumber contributions of the signal, so as to better characterize the shallowest sources and/or their shallowest parts (Barone et al., 2019). Also, we provide ∂ p f / ∂ z p by antitransforming the following relation (Blakely, 1996): ℱ [ ∂ p f ∂ z p ] = | k | p  ℱ [ f ]   . (13)

A homogeneous function f satisfies the following scaling law in the region R: f ( t x , t y , t z ) = t n f ( x , y , z ) , (20) where, t > 0 and n ∈ ℝ is the homogeneity degree of f . When f is continuously differentiable, Euler’s theorem describes the homogeneity, as follows (e.g., Olmsted, 1961; Vitale and Fedi, 2020): ∇ f ( r ) · ( r − r 0 ) = − n f ( r ) , (21) where, r 0 represents the position of the field source.

n often expresses the falloff rate of the anomaly with the scale. It is an integer and constant for the fields generated by ideal sources (e.g., Fedi et al., 2015). Castaldo et al. (2018) have shown that a point-spherical source with a hydrostatic Δ P (i.e., the Mogi model) generates a homogeneous ground deformation field with n = − 2 . Accordingly, − 2 ≤ n ≤ 1 characterize the fields generated from the 3D concentrated to nonconcentrated sources, respectively.

The homogeneity properties reflect the type of field source. In the case of the harmonic deformation field, which is the first-order derivative of a Newtonian potential, we get the homogeneity degree of the source n s , n s = n − 1 , (22) and, in turn, the structural index of the source N (e.g., Fedi, 2007): N = − n s . (23) N takes on an important meaning since it represents a source parameter characterizing its geometry. Following Fedi et al. (2015), we define N as an integer number varying from 0 to 3 going from a 3D nonconcentrated to concentrated sources. For example, N ∼ 3 characterizes a point-spherical source or mathematically equivalent body concentrated in its center; N ∼ 2 indicates source distribution supported on its top, as the pipe-like one, and finally, N ∼ 1 is related to sill- and dike-like sources, for which their edges are representative of the source distribution.

We show the validity of the multi-scale approach by applying Multiridge and ScalFun methods to synthetic deformation fields. We simulate both harmonic and nonharmonic fields using a Finite Element (FE) approach (COMSOL Multiphysics software), which allows us to model different scenarios or to accurately reproduce fields satisfying the Eq. 11 and Laplace’s equation (Eq. 2).

We consider the simple and geometrically finite sources, characterized by hydrostatic Δ P and embedded in a homogeneous elastic half-space, to model the following cases: (I) the sphere, with radius r = 0.5  km, which is the simplest approximation of a magma chamber; (II) the ascent-pipe, here represented by an ellipsoid with semi-axes r x = r y = 0.5  km and r z = 2.5 km (aspect ratio of the ellipsoid = 0.2 ); (III) the sill emplacement, here approximated by using a rectangular prism with sides L x = L y = 5 km and L z = 1  km. Furthermore, we take into account the additional synthetic cases with different model settings: (IV) the layered half-space test, where a sill-like source is embedded in a heterogeneous medium, characterized by three layers with an increasing E distribution vs. depth (2 GPa for the layer 0–1 km, 7 GPa for 1–3 km, and 10 GPa for the last one); (V) the noisy case, where 30% (with respect to the maximum value) of the random noise is added to the latter simulated field (Supplementary Figure S1); and (VI) the multi-source scenario, where the over-pressurized spherical and pipe-like sources contribute to the same deformation anomaly. The other FE modeled parameters are reported in Supplementary Tables S1, S2.

We start from the analysis of the vertical deformation components ( p = 0 ) of the first three cases (Figures 1A–C). We calculate the 3D spatial multi-scale datasets by the upward continuation of the deformation fields from 0 to 15 km, with a 0.4 km step, and extract the profile (black dashed lines in Figures 1A–C) passing through the maximum value of the anomalies. We use the multiridge subsets (blue and black dots in Figures 1D–F) to represent the ridges (red continuous lines in Figures 1D–F). For the sphere (I), we note that the ridges intersection indicates a source located at x = 50  km and -3 km depth (Figure 1D). Instead, for the pipe (II), the ridges above 5 km scale line up along the same straight line of the previous one and the intersection occurs at x = 50  km and -2.8 km depth (Figure 1E). We observe a similar behavior for the sill (III), where we consider the upward continued field above 8 km scale to retrieve the ridges intersection at the point x = 50  km and -3.1 km depth (Figure 1F). For these cases, the results of the ScalFun method indicate (Figures 1G–I): τ ∼ − 2 and n ∼ − 2 according to the Eq. 25, and N ∼ 3 for Eqs 22, 23.

We also analyze these cases using the vertical differentiation of the ground deformation field. We calculate the first-order vertical derivatives ( p = 1 ) of the deformation field for the sphere (Figure 2A) and of the pipe (Figure 2B), and the second-order vertical derivatives ( p = 2 ) for that of the sill (Figure 2C). Then, we create a 3D spatial multi-scale dataset up to 5 km, with a 0.4 km continuation step, and extract the profile (black dashed line in Figures 2A–C) passing through the maxima of the vertical derivatives. For the sphere (I), the ridges intersection again occurs at x = 50  km and -3 km depth (Figure 2D). The ScalFun method indicates (Figure 2E): τ 1 ∼ − 3 and n ∼ − 2 according to Eq. 25, and N ∼ 3 for Eqs 22, 23. For the pipe (II), we achieve two sets of ridges with different intersections, one of which indicates the point with x = 50  km and -0.51 km depth, that is near the source top (Figure 2F). The ScalFun method characterizes this solution with   N ∼ 2 , since τ 1 ∼ − 2 and n ∼ − 1 (Figure 2G). Finally, for the sill (III), the ridges converge at two different points with coordinates ( x   , z ) = ( 47.6   ,   − 2.4 )  km and ( x   ,   z ) = ( 52.4   ,   − 2.4 )  km, near the body edges (Figure 2H). The ScalFun method characterizes these solutions (Figure 2I) with τ 2 ∼ − 2 , n ∼ 0 , and N ∼ 1 .

Similar considerations occur by analyzing the E-W deformation components and their vertical derivatives, whose results are reported in Supplementary Figures S2, S3.

For what concern the layered scenario (IV), we analyze the vertical component (Figure 3A) of the deformation ( p = 0 ) and observe the ridges intersection at x = 50 km and -1.9 km depth (Figure 3D). The ScalFun method indicates (Figure 3E): τ ∼ − 2 and n ∼ − 2 according to Eq. 25, and N ∼ 3 for Eqs 22, 23. We still analyze the vertical deformation (Figure 3B) component ( p = 0 ) modeled for the noisy test (V) and retrieve a single intersection at x = 50 ± 0.4 km with -1.9 ± 0.15 km depth (Figure 3F). Also, the ScalFun method characterizes this solution with   N ∼ 3 , since τ ∼ − 2 and n ∼ − 2 (Figure 3G). Finally, for the multisource case (VI), we analyze the second order vertical derivative ( p = 2 ) of the vertical deformation (Figure 3C) and retrieve two sets of ridges converging at two different points with coordinates ( x , z ) = ( 46.4   ,   − 4.9 )  km and ( x ,   z ) = ( 53.6   ,   − 1.8 )  km (Figure 3H). The ScalFun method characterizes both the solutions (Figure 3I) with τ 2 ∼ − 4 , n ∼ − 2 , and N ∼ 3 .

A summary of the performed tests is reported in Table 1.

We focus on the conditions for which the pressurized volcanic bodies generate a harmonic and homogeneous ground deformation field. In this framework, we show that Laplace’s equation and the homogeneity law are surely satisfied in the case of hydrostatically pressurized sources embedded in a homogeneous elastic half-space. It is well-known that these conditions are just an approximation of the real volcanic scenario (e.g., Gudmundsson, 2006; Hickey et al., 2016), as well as the case of many AMs. One of these is the Mogi model (Mogi, 1958). Its analytical expression is represented by the gradient of a Newtonian potential in the form ϕ = 1 / r (Castaldo et al., 2018), which is a harmonic and homogeneous function. However, the multi-scale approach is also appropriate to model some classes of nonharmonic fields. Indeed, Castaldo et al. (2018) have applied this methodology to the deformation fields generated by spherical sources embedded in a heterogeneous and layered elastic half-space by retrieving negligible errors on the estimated parameters. Also, in this work, we propose the layered half-space model setting by considering a different source, represented by a sill-like one.

Here, we show that the proposed multi-scale strategy is based on two methods: The Multiridge method provides information on the center of deformation sources by evaluating the ridges intersection, which can also occur, for the pipes, sills, and dykes, at their top and/or boundaries. The ScalFun method characterizes the source shape through the estimates of the structural index N . For example, with N ∼ 3 , we are analyzing the fields generated by a mathematically equivalent body concentrated in its center, while N ∼ 2 characterizes the pipe-like sources and ridges intersections related to their top. N ∼ 1 is finally related to the edges or top vertexes of the planar bodies, such as sills- and dykes-like sources (e.g., Reid et al., 1990; Florio et al., 2014). We remark that these methods are based on the properties of continuous harmonic and homogeneous functions. The applicability depends evidently on the 2D distribution of deformation measurements, which should be spatially dense, and this occurs with the DInSAR dataset. On the contrary, the Leveling and Global Navigation Satellite System (GNSS) data generally lead to unreliable interpretation of the multi-scale results.

We perform several tests on the simulated harmonic and homogeneous fields. In particular, we analyze the vertical and E-W deformation components. Although the simulated N-S one should be available, we do not consider it since this study is conceived for the analysis of DInSAR measurements, which do not provide this information, due to the near-polar orbit of the SAR sensors. However, the harmonic and homogeneous properties of the deformation field also allow the involvement of the LOS components in the framework of the multi-scale procedure (Castaldo et al., 2018). In addition to the field component, we employ their vertical differentiations of order p since the harmonic properties allow us to compute them in the wavenumber domain. This operator increases the analysis resolution by emphasizing the signals high-wavenumber contributions that are related to the shallowest features of the sources, as the pipe’s top and the sill’s edges. It is important to note that the vertical differentiation also enhances the high-wavenumber noise of the dataset, which affect the lower scales. Therefore, when the signal-to-noise ratio is low, it is recommendable to perform the multi-scale analysis at rather large scales (Fedi et al., 2009). The scale-invariance vs. p is a property of the harmonic and homogeneous fields and, by using the vertical differentiation operator, we can verify this condition. The sphere test (I) represents the simplest example of this property, where both for the field components (vertical and E-W) and for their vertical derivative analysis, we observe ridges converging to the center of the body, with the structural index N ∼ 3 . Differently from the previous source, the pipe case (II) shows that the ridges at larger scales converge to the source center with N ∼ 3 , while at lower scales, they tend to intersect at its top (shallowest singular point). Accordingly, the vertical derivative analysis shows a ridges intersection nearby its top with N ∼ 2 . The N estimates and the related depths of each ridges intersection (center and top) allow us to interpret the responsible source as a vertically elongated one (e.g., pipe-like). Similarly, in the sill test (III), the ridges at larger scales intersect at the source center finding N ∼ 3 , while at lower scales, they are strongly influenced by the body edges. Indeed, by analyzing the vertical derivative of the field components, we retrieve the ridges intersections nearby the prism edges, with N ∼ 1 . We can correctly interpret the deformations generated by a planar body (e.g., sill-like).

It is worth noting that our simulations have shown that the Multiridge and ScalFun methods are valid tools to study the harmonic deformation fields of homogeneous sources. They have provided, in the worst cases, estimates on the depth to the source center affected by ∼ 7% error with respect to the real value and on its top/edge position with 4% of error with respect to the source extent. The retrieved errors in controlled conditions may mainly depend on the replication of the harmonic field properties and the data sampling.

We obtain comparable errors on source parameters by analyzing the nonharmonic fields, generated by a sill embedded in a layered half-space (IV) and with the same random noise (V). In both the simulations, we retrieve the source depth to the center with only 5% error with respect to the expected values. This finding demonstrates that the Multiridge and ScalFun methods in the same conditions can be extended to the analysis of nonharmonic fields, as in the case of layered heterogeneous medium with weak physical parameters contrasts. Furthermore, we observe that the noisy dataset does not affect the multi-scale results, at least for the source center estimate.

Finally, with the last test (VI), we show the flexibility of the proposed methodology, which can be also used to analyze a multi-source scenario. This occurs by exploiting the vertical derivatives of the field that allow us to detect the parameters of both the modeled sources, with a maximum estimated error of 10% at the depth to the center.

We use the Multiridge and ScalFun methods for analyzing three different American volcanoes. We show the applicability of the multi-scale approach in different volcanic environments, since the same geodynamics settings and geochemical features do not characterize the considered systems.

Unlike the previous synthetic cases, we evaluate R 2 parameter for each ridge and the related intersection uncertainties, since different signal contributions (i.e., instrumental/experimental errors, high-wavenumber noise in the dataset, effect due to the draped-to-flat upward continuation at the areas with high slopes, and effects due to other nonvolcanic sources), can affect the multi-scale results. We can avoid some of these contributions that are predominant at the lower scales by focusing the analysis on the rather larger scales of the multi-scale dataset.

We perform the first application at Uturuncu volcano (Bolivia), which is located on the largest updoming zones on the Earth. In this framework, it is well known that the deep APMB source has a crucial role for the observed multi-decade decimeter-scale deformation signal (Allmendinger et al., 1997; Schmitz et al., 1997; Schilling et al., 1997; Chmielowski et al., 1999; Pritchard and Simons, 2002; Zandt et al., 2003; Pritchard and Simons, 2004; Sparks et al., 2008; Fialko and Pearse, 2012; del Potro et al., 2013; Henderson and Pritchard, 2013; Hickey et al., 2013; Walter and Motagh, 2014; Ward et al., 2014; Comeau et al., 2015; Comeau et al., 2016; Perkins et al., 2016; Gottsmann et al., 2017; Kukarina et al., 2017; Ward et al., 2017; McFarlin et al., 2017; Henderson and Pritchard, 2017; Pritchard et al., 2018; Barone et al., 2019), despite the evidences of geophysical, seismological, and petrological data analyses also highlighting the hypothesis of shallow sources (Sparks et al., 2008; Jay et al., 2012; Muir et al., 2014; Alvizuri and Tape, 2016; Comeau et al., 2016; Lau et al., 2018; Barone et al., 2019). We study the cumulative vertical deformation related to the transient 2007 unrest, for which the AMPB is supposed to be not responsible of the observed ground deformation field (Barone et al., 2019). Our results infer the depth to the source center of 4.5 km b.s.l., with a structural index N ∼ 3 , pointing out the accordance with the supposed scenario of the existence of a shallow hydrothermal active system beneath Uturuncu volcano (e.g., Gottsmann et al., 2017; Lau et al., 2018). This also matches with the low resistivity zones, related to the saline fluids of the magmatic and/or meteoric origin (Comeau et al., 2015; Comeau et al., 2016), and with other geophysical and petrological data analysis (Sparks et al., 2008; Jay et al., 2012; Muir et al., 2014; Alvizuri and Tape, 2016). However, the multi-scale approach does not provide any information on the nature of the shallow source below Uturuncu volcano, for which the hypothesis of a hydrothermal system is still the most reliable option explaining the geophysical, seismological, and geodetic anomalies (Hudson et al., 2022).

We perform the second analysis at Okmok volcano (Alaska, United States), where the multi-scale approach has been already applied to the single LOS component. In this work, we consider the same dataset of Castaldo et al. (2018), but we apply the Multiridge and ScalFun methods to its second-order vertical derivative. Here, the higher resolution analysis shows the scale-invariance vs. p . of the multi-scale results. Indeed, they do not significantly change by increasing p , we estimate the depth to the source center, with N ∼ 3 , at 3.1 km b.s.l. in the case of p = 2 , which agrees with those obtained with p = 0 (3.4 km b.s.l. and N ∼ 3 ). It follows that the analyzed dataset is in good approximation with the harmonic and homogeneous field, and the analyzed source consists of a 3D concentrated body. The retrieved source (Table 3) agrees with those proposed by several authors modeling the 1997–2008 intereruptive stage (e.g., Lu et al., 2005; Biggs et al., 2010; Masterlark et al., 2010; Masterlark et al., 2012; Xue et al., 2020), despite the fact that our analysis is related to a smaller time interval (2003–2004).

Finally, we apply the multi-scale methods to evaluate the depth, horizontal location, and morphological features of the deformation source related to the 2013 unrest at Fernandina volcano. Different from the previous real cases, we analyze the E-W deformation component. We specify that the coherent measurement points (Figure 8A) are not well distributed at the summit caldera (uncorrelated signals), where the maximum vertical deformation is expected. A better data distribution is observed at the volcano flanks, where the maxima and minima of the E-W deformation anomaly are located. Our findings point out two clusters of solutions with different N values and we associate both to a unique volcanic body. The deformation pattern is interpreted as generated by a vertically elongated source with 1.5–1.7 km b.s.l. depth to the center ( N ∼ 3 ) and 0.5–0.6 km b.s.l. depth to the top ( N ∼ 2 ). We specify that our findings are robust and give recurring solutions by analyzing different profiles and different p -order cases. We compare them with the geodetic studies based on the optimization/inverse methods (Table 4), finding that the suggested volcanic source is compatible with the scenario of Bagnardi and Amelung (2012), where the magmatic feeding system consists of connected deep and shallow sill-like reservoirs. Our vertical source just images a conduit body for magma or fluid ascent from the deeper to shallower portions of the volcanic system. The retrieved results also agree with the pipe-like source proposed by Pepe et al. (2017). However, we note that our and latter works refer to the modeling of uplift that has involved the central part of the summit caldera. Instead, Bagnardi and Amelung (2012) have mainly investigated the unrest events that have led to the intrusion of lateral dykes (Jonsson et al., 1999; Chadwick et al., 2011).

We conclude the work by listing the main features of the proposed multi-scale approach.

The methodology represents an alternative procedure with respect to the optimization/inverse methods, which, differently from our approach, generally need to fix information, such as the medium elastic parameters, the source physical parameters, and its geometry.

The methods provide an unconstrained source geometrical model, while they do not provide information on the source physical parameters. Their results can also be used as constraints for the inverse/optimization modeling procedures, aiming at estimating a more complete set of source parameters.

The methods are suitable to:

• Analyze any deformation components;

and their results:

• Do not depend on lower variations of the physical parameters of the elastic half-space;

Although the presented multi-scale approach is restricted to the analysis of sources as point-spherical, pipe-, and sill-like ones, future developments will concern the use of Multiridge and ScalFun methods in the case of sources with complex geometry.