Compound dislocation models (CDMs) for volcano deformation analyses

SUMMARY

Volcanic crises are often preceded and accompanied by volcano deformation caused by magmatic and hydrothermal processes. Fast and efficient model identification and parameter estimation techniques for various sources of deformation are crucial for process understanding, volcano hazard assessment and early warning purposes. As a simple model that can be a basis for rapid inversion techniques, we present a compound dislocation model (CDM) that is composed of three mutually orthogonal rectangular dislocations (RDs). We present new RD solutions, which are free of artefact singularities and that also possess full rotational degrees of freedom. The CDM can represent both planar intrusions in the near field and volumetric sources of inflation and deflation in the far field. Therefore, this source model can be applied to shallow dikes and sills, as well as to deep planar and equidimensional sources of any geometry, including oblate, prolate and other triaxial ellipsoidal shapes. In either case the sources may possess any arbitrary orientation in space. After systematically evaluating the CDM, we apply it to the co-eruptive displacements of the 2015 Calbuco eruption observed by the Sentinel-1A satellite in both ascending and descending orbits. The results show that the deformation source is a deflating vertical lens-shaped source at an approximate depth of 8 km centred beneath Calbuco volcano. The parameters of the optimal source model clearly show that it is significantly different from an isotropic point source or a single dislocation model. The Calbuco example reflects the convenience of using the CDM for a rapid interpretation of deformation data.

Geomechanics, Kinematics of crustal and mantle deformation

1 INTRODUCTION

Since their early development, analytical models of crustal deformation have been used to interpret volcano deformation processes (Segall 2010; Lisowski 2007). Despite their simplicity, the widespread application of these models has proven that they are notably useful tools as first-order approximations for studying physical processes at volcanoes. These models are based entirely on the mathematical theory of elasticity applied to the deformation of the Earth's surface and its interior due to forces or dislocations (Love 1944). The diversity of volcano deformation processes (see Segall 2010) requires a flexible representation of the wide variation in source geometries in elastic media. Despite being in use for decades, some basic analytical models based on either concepts of dislocations (Okada 1985, 1992) or forces (Davis 1986) include numerical artefacts and structural imperfections that are to be addressed in this paper.

The simplest models of planar volcanic sources, namely, dikes and sills, are those based on rectangular dislocations (RDs) with only a prescribed uniform opening (Segall 2010). One limitation of the existing RD solutions (Okada 1985, 1992) is that they lack full rotational degrees of freedom. The reason for this limitation is that in these solutions, two parallel edges of an RD can be dipping at any arbitrary angle, but the other two edges, which are perpendicular to the first two, have to be parallel to the free surface. The existence of the artefact singularities along the edges and below and above the vertices is the other problem with these RD solutions. Okada (1992) reported this problem and addressed it in a small neighbourhood of the artefact singularities. The predefined size of the neighbourhood, however, makes the Okada (1992) analytical solution scale dependent, which means that the solution depends on the dimensions of the RDs and on the distribution of the calculation points in a model. Bradley & Segall (2012) addressed the artefact singularities along the edges of RDs but not below and above their vertices. Nevertheless, neither of these attempts fully addressed the problem of artefact singularities.

The widely known point-source models of pressurized cavities are based on force dipoles. The simplest case in this group of cavity models is the centre of dilatation, referred to as the Mogi model (Mogi 1958). The Mogi model is a point-source approximation that is suitable only for deep spherical sources (Segall 2010). Comparatively, the Davis (1986) point ellipsoidal cavity model is a point-source approximation for deep triaxial ellipsoidal sources, which also includes the spherical cavities as a special case. Nonetheless, compared to the Mogi model, the Davis (1986) solution has not yet been adopted as a routine modelling tool.

These issues also affect other analytical and numerical solutions. It has been shown that a certain combination of a few dislocations can reproduce the same displacement and stress fields as that of an isotropic pressurized point source (Bonafede & Ferrari 2009). Had this model been generalized, it could have replaced most of the volcano deformation analytical models, which are based on either concepts of dislocations or forces.

Given the problems and imperfections in the aforementioned analytical solutions, developing new solutions that first enhance the functionality and improve the performance of these analytical solutions and then make the application of such solutions more convenient is of relevance.

In the following, we first develop new solutions for the surface, internal displacement and stress fields of RDs with full rotational degrees of freedom in both full-space and half-space. We address the problem of artefact singularities along the edges and below and above the vertices of the RDs. Using our RD solutions, we develop a compound dislocation model (CDM) as an alternative generalized source model of pressurized cavities and compare the CDM with analytical and numerical solutions. Finally, we apply the CDM to the 2015 April 22 co-eruptive displacements of Calbuco volcano in Chile. The application of the CDM to the Calbuco case study supports a vertically elongated lens-shaped source, which is significantly different from both spherical and planar volcanic sources. We show that the CDM can be simply integrated in any optimization algorithm to invert the observed displacements for the unknown parameters of various deformation sources.

2 METHODS

In this section, we develop new analytical solutions for RDs in a full-space and half-space based on the earlier works of Comninou & Dundurs (1975) and Nikkhoo & Walter (2015). We then elaborate a CDM, which is a generalization of the Bonafede & Ferrari (2009) dislocation source.

2.1 New analytical, artefact-free solutions for rectangular dislocations

We develop new artefact-free solutions for the displacement and stress fields associated with RDs with full rotational degrees of freedom in an elastic full-space and half-space.

In the first step, we construct two different but mathematically equivalent configurations of angular dislocations in a full-space (Fig. 1a). The double dislocation lines in Fig. 1(a) show the lines of artefact singularities. These lines in the first configuration are extended in directions that completely differ from those in the second configuration (Fig. 1a). The edges of the RD and their perpendicular bisectors, as well as the bisectors of the RD external angles, mark the boundaries between the two complementary partitions (white and dark shadings in Fig. 1b). Clearly, the artefact singularities of each configuration remain inclusively inside one of the two partitions. Additionally, the orthogonal projection of any calculation point onto the RD plane falls into only one of the two partitions. In the new RD solution, we first specify the partition that includes the projection of the calculation point, and then we use the configuration that is located in the other partition. In this way, we always avoid all of the artefact singularities in our RD solution. Moreover, we calculate the Burgers function associated with the RD as the sum of the Burgers functions of two triangular dislocations or four angular dislocations that form exactly the same RD (Nikkhoo & Walter 2015).

Figure 1.

Addressing the problem of artefact singularities and numerical instabilities for rectangular dislocations (RDs). (a) Two different configurations for the RD, named by ABCD, are depicted. Apart from the placements of the semi-infinite dislocation lines, these two configurations are equivalent and can be used interchangeably. The calculation point P in the plane of the RD intersects one of the dislocation lines of the configuration I. For the ‘alternative’ configuration II, the point P is located at the maximum possible distance from the dislocation lines. Through the selection of the appropriate configuration, any artefact singularity or numerical instability can be avoided. (b) The partitions in the plane of the RD formed by the RD angle bisectors and the RD perpendicular bisectors (the solid blue lines). The red short dashed lines and the green dashed lines represent the lines of singularities of the RD configurations I and II in (a), respectively. Each line of singularity lies in only one of the partitions. To avoid artefact singularities and numerical instabilities, we use configuration I for the calculation points in which their orthogonal projection onto the RD plane falls in the dark partitions, and for the other calculation points, we use configuration II. Either configuration can be used for the points inside the RD.

In the second step, we apply the method of images directly to the new RD full-space solution and develop an artefact-free solution for an RD in a half-space (see Nikkhoo & Walter 2015). For this purpose, using the artefact-free solution for RDs in a full-space, we calculate the main and image dislocation contributions. The sum of these two together with the free surface potential function contribution (Comninou & Dundurs 1975) completes the solution for the RD in a half-space. Using this technique, we are able to calculate displacement and stress fields at depth and on the surface of a half-space.

As the final step, using the Comninou & Dundurs (1975) solution for an angular dislocation in a half-space, we develop another solution only for the surface displacement and stress fields associated with an RD with full rotational degrees of freedom in a half-space. Implementing two configurations and switching between them in a method similar to the full-space case, we remove the artefact singularities that appear when one of the RD edges is in the vicinity of or breaches the free surface (see Nikkhoo & Walter 2015).

The geometrical parameters of the Okada (1985, 1992) and the new RD solutions are illustrated in Fig. 2. The strike angle (α) and the dip angle (δ) determine the orientation of the plane upon which the RD is located. We introduce a new parameter that controls the orientation of the RD in this ‘extended RD plane’. This parameter, which we call the ‘plunge angle’ (θ), is the angle between the upper edge of the RD and the intersection of the free surface with the extended RD plane (see Fig. 2). The possibility of having a plunging upper edge in the new RD solutions allows them to possess full rotational degrees of freedom. In other words, the new RD solutions that we developed here extend the well-known Okada (1985, 1992) solutions by addressing the numerical artefacts and the geometrical limitation problem in their structure.

Figure 2.

The geometry of the Okada (1985, 1992) solutions (dashed line) and the new solution (solid line) for a rectangular dislocation (RD) in a half-space. E₀, N₀ and d are the coordinates of the RD centroid, and L and W are the RD length and width, respectively. δ and α, which define the orientation of the RD plane, are the dip and strike angles, respectively. The plunge angle, θ, which only exists in the new solution, can change the orientation of the RD within the RD plane. Note that the upper edge of the Okada (1985, 1992) RD solutions is always parallel to the free surface. Because of having the plunge angle in its geometrical configuration, the new RD solution possesses full rotational degrees of freedom.

2.2 The compound dislocation model

The displacement field of three mutually orthogonal square dislocations with the same prescribed uniform opening is equivalent to that of the Mogi model (Bonafede & Ferrari 2009). We generalize this concept and use three mutually orthogonal RDs to develop a source model that can represent planar and volumetric sources of various aspect ratios. As a result of utilizing the RD solutions of Section 2.1 in its configuration, this generalized model may possess any arbitrary size and orientation in space (Fig. 3). We refer to this new model as the compound dislocation model (CDM).

Figure 3.

The geometrical structure and rotation angles of a compound dislocation model (CDM) composed of three rectangular dislocations (RDs), labelled as ‘A’, ‘B’ and ‘C’ (the yellow, green and blue planes). The origin of the XYZ coordinate system is at the Earth's surface and the positive X, Y and Z axes point to the east, north and up, respectively. The origin of the xyz coordinate system is located on the CDM centroid and the x, y and z axes are normal to ‘A’, ‘B’ and ‘C’, respectively. The xyz coordinate system is fixed to the CDM. The a, b and c are the semi-axes of the CDM along the x, y and z axes, respectively. The ω_X, ω_Y and ω_Z are angles of rotation about the X, Y and Z axes, respectively. Positive values of these angles correspond to clockwise rotations. (a) The initial orientation of the CDM. (b) The orientation of the CDM after applying ω_X. (c) The orientation of the CDM after applying ω_X and ω_Y. (d) The final orientation of the CDM. The numerical values of the rotation angles for this given example are denoted on the bottom right of each panel.

In the following, we adopt two coordinate systems (see Fig. 3). The x, y and z axes, which lie on the CDM axes, establish a Cartesian coordinate system, which its origin is located on the CDM centroid. The lengths of the CDM semi-axes along the x, y and z axes (Fig. 3a) are equal to a, b and c, respectively. The RDs labelled as ‘A’, ‘B’ and ‘C’ are perpendicular to the x, y and z axes and their dimensions are (2b, 2c), (2a, 2c) and (2a, 2b), respectively. The origin of the XYZ coordinate system is at the Earth's surface and the positive X, Y and Z axes point to the east, north and up, respectively. Therefore, the XYZ coordinate system is a local Earth-fixed Cartesian coordinate system and the XY plane represents the free surface. We use the XYZ coordinate system as the reference frame to define the orientation and location of the CDM and the xyz coordinate system. The clockwise rotations of ω_X, ω_Y and ω_Z about the X, Y and Z axes, respectively, specify the orientation of the CDM in space (see Fig. 3). Note that before applying the rotations the xyz and XYZ coordinate systems are coincident (see Fig. 3a) and after applying the rotations these coordinate systems still share the same origin (see Fig. 3d). The translation vector that moves the CDM to its final location in space is (X₀, Y₀, −d), where X₀, Y₀ and d are the east and north coordinates and the depth of the CDM centroid, respectively, in the XYZ coordinate system. Considering the angles of rotation in Figs 2 and 3 it can be shown that ω_X = θ_A, |$\omega _Y = \frac{\pi }{2}-\delta _A$| and ω_Z = α_A, where θ_A, δ_A and α_A are the plunge angle, dip angle and strike angle of the first RD of the CDM, which is labelled as ‘A’ in Fig. 3. The RDs that form the CDM all have the same amount of uniform opening u. Therefore, the CDM has a total number of 10 parameters that include the CDM centroid location (X₀, Y₀, d), the CDM rotation angles (ω_X, ω_Y, ω_Z), the CDM semi-axes (a, b, c) and the opening (u).

If we consider the long-wavelength components of the displacement signal associated with a CDM in the far field only, then the model will effectively be a point source (Aki & Richards 2002). Then, the CDM can be considered as a system of force dipoles located at its centroid, with a moment tensor equal to the integral of the moment density over the CDM surface (Aki & Richards 2002). Due to the linearity of the integral operator, the moment tensor of the effective point source is equal to the sum of the moment tensors of the individual RDs in the model. Without loss of generality, we derive the moment tensor of the CDM in Fig. 3 with respect to the ‘CDM-fixed’, xyz coordinate system. For this CDM, the areas of the RDs ‘A’, ‘B’ and ‘C’ are equal to S_A = 4bc, S_B = 4ac and S_C = 4ab, respectively, and u is the uniform opening of all three RDs. Based on the moment tensor density of a tension crack after Aki & Richards (2002), the moment tensor of the CDM in the xyz coordinate system will be

\begin{eqnarray} M^{\text{CDM}} &= &\left(\begin{array}{c@{\quad}c@{\quad}c}\left(\lambda +2\mu \right) u S_{A} & 0 & 0\\ 0 & \lambda u S_{A} & 0\\ 0 & 0 & \lambda u S_{A}\end{array}\right) \nonumber \\ &&+ \left(\begin{array}{c@{\quad}c@{\quad}c}\lambda u S_{B} & 0 & 0\\ 0 & \left(\lambda +2\mu \right) u S_{B} & 0\\ 0 & 0 & \lambda u S_{B}\end{array}\right) \nonumber\\ &&+ \left(\begin{array}{c@{\quad}c@{\quad}c}\lambda u S_{C} & 0 & 0\\ 0 & \lambda u S_{C} & 0\\ 0 & 0 & \left(\lambda +2\mu \right) u S_{C}\end{array}\right), \end{eqnarray}

(1)

where μ, known as the shear modulus, and λ are the Lamé coefficients.

The contraction of eq. (1) is

\begin{eqnarray} M^{\text{CDM}} = u \left(\begin{array}{ccc}\lambda S + 2\mu S_{A} & 0 & 0\\ 0 & \lambda S+2\mu S_{B} & 0\\ 0 & 0 & \lambda S+2\mu S_{C}\end{array}\right), \end{eqnarray}

(2)

where S = S_A + S_B + S_C = 4(ab + bc + ac).

Moreover, eq. (2) can be written as

\begin{eqnarray} &&{M^{\text{CDM}} =}\nonumber\\ &&\left(\begin{array}{c@{\quad}c@{\quad}c}\lambda \: \Delta V + 2\mu \: \Delta V_{A} & 0 & 0\\ 0 & \lambda \: \Delta V+2\mu \: \Delta V_{B} & 0\\ 0 & 0 & \lambda \: \Delta V+2\mu \: \Delta V_{C}\end{array}\right),\nonumber\\ \end{eqnarray}

(3)

where ΔV_A = uS_A, ΔV_B = uS_B and ΔV_C = uS_C are the potencies of the three RDs that compose the CDM and ΔV = ΔV_A + ΔV_B + ΔV_C is the potency of the CDM. The potency is the product given by area × slip, which describes the strength of a dislocation source and holds the dimensions of volume change (Aki & Richards 2002). The potency of a tensile dislocation is given by area × opening. From eq. (3), we can write

\begin{eqnarray} \Delta V &=& \frac{M_{kk}^{\text{CDM}}}{3K}, \end{eqnarray}

(4)

where |$M_{kk}^{\text{CDM}}$| is the sum of the diagonal components of M^CDM and |$K=\lambda +\frac{2}{3}\mu$| is the bulk modulus. The fact that |$M_{kk}^{\text{CDM}}$| is the first invariant of M^CDM shows that eq. (4) is independent of the coordinate system.

In the case of a CDM that simulates an inflating crack-like or volumetric cavity, the potency ΔV represents the volume available to host new fluids intruding into the cavity from the outside (Bonafede & Ferrari 2009). The potency, however, must not be confused with the actual volume change δV, which represents the expansion that a deformation source applies to the surrounding elastic medium (see Aki & Richards 2002; Kumagai et al. 2014; Ichihara et al. 2016).

2.2.1 The point CDM

We develop another source model by replacing the finite RDs in the CDM configuration with point tensile dislocations (Okada 1992). This model is a point-source version of the CDM and we refer to it as the point CDM. Similar to other point sources such as the Mogi (1958) and Davis (1986) models the point CDM is exactly equivalent to a single moment tensor that represents a system of three mutually orthogonal force dipoles in a half-space (see Aki & Richards 2002; Bonafede & Ferrari 2009). In the far field, the point CDM and CDM are equivalent and therefore share the same moment tensor in eq. (3). The point CDM, from a geometrical point of view and by using the concept of the calculus of infinitesimals, can be visualized as the CDM in Fig. 3(a), but with the semi-axes a, b and c replaced with ‘differential’ semi-axes da, db and dc, respectively. Because the point CDM is inherently a point source, the differential semi-axes tend to zero, that is |$da \longrightarrow 0$|⁠, |$db \longrightarrow 0$| and |$dc \longrightarrow 0$|⁠, respectively. However, the potencies of the point dislocations, that form the point CDM, and the ratios of the differential axes, that determine the geometrical shape of the point CDM, can be different from zero. The parameters that specify the location and spatial orientation of the point CDM are the same as those of the CDM. These parameters and the potencies of the point dislocations are the 9 parameters that uniquely specify the point CDM. Considering the definition of potency, it is straightforward to show that the axes ratios of the point CDM can be calculated as ratios of the point dislocation potencies.

3 EVALUATION OF THE COMPOUND DISLOCATION MODEL

In this section, we evaluate the CDM through comparisons with analytical and numerical solutions. We compare the CDM to point and finite ellipsoidal sources and assess the ability of the CDM to simulate deformation sources of various geometries and depths in the near field and far field.

3.1 The CDM and point ellipsoidal source comparison

Moment tensors reveal the far-field behaviour of their corresponding deformation sources (Aki & Richards 2002) and can therefore be used for a quantitative assessment of these deformation sources. To better understand the performance of the CDM, we compare its moment tensor to that of the Davis (1986) point ellipsoidal cavity. Because the CDM and point CDM share the same moment tensors (see Section 2.2), the results of the comparison in this section, which are presented in Fig. 4, apply to both the CDM and point CDM.

Figure 4.

The moment tensor spectrum of the compound dislocation model (CDM) and point ellipsoidal cavity models in a Poisson solid (μ = λ) as a function of M₂/M₁ and M₃/M₁, where M₁, M₂ and M₃ are the eigenvalues of the moment tensors and M₁ ≥ M₂ ≥ M₃ holds (Trasatti et al. 2009). The lower right-triangle area on the plot represents the universal set of moment tensors. The CDM possesses a wider moment tensor spectrum (light and dark grey area), which encompasses that of the ellipsoidal cavity (dark grey area) as a subset.

The equivalent moment tensor M^E for the point ellipsoidal cavity after Davis (1986) is

\begin{eqnarray} M^{E} = \sigma ^{T} V, \end{eqnarray}

(5)

where σ^T is the transformation stress after Eshelby (1957) and V is the volume of the ellipsoid (see Segall 2010).

To calculate σ^T, we implement the Carlson (1995) algorithms for numerical computation of elliptic integrals and the Eshelby (1957) transformation strain equations in a MATLAB function. Using eqs (2) and (5), we calculate the moment tensors M^CDM and M^E that cover the entire range of various aspect ratios of the CDM and the Davis (1986) point ellipsoidal cavity. By applying the Trasatti et al. (2009) convention to the eigenvalues of the calculated moment tensors, we can represent the domains of moment eigenvalue ratios that correspond to the CDM and point ellipsoidal cavity in a Poisson solid, where μ = λ (Fig. 4). Clearly, the CDM possesses a wider domain, which as its subset also encompasses the domain of the ellipsoidal cavity model. This implies that corresponding to any arbitrary point ellipsoidal cavity, there is a CDM with the same moment tensor; however, the reverse statement does not necessarily hold. To derive the direct relationship, we assume that the diagonal moment tensor

\begin{eqnarray} M &=& \left(\begin{array}{c@{\quad}c@{\quad}c}M_{11} & 0 & 0\\ 0 & M_{22} & 0\\ 0 & 0 & M_{33}\end{array}\right) \end{eqnarray}

(6)

represents a point ellipsoidal cavity with arbitrary dimensions in the same xyz coordinate system introduced in Section 2.2. The potencies of the RDs that compose a CDM with the same moment tensor M are uniquely determined through comparing eqs (3) and (6) as

\begin{eqnarray*} \Delta V_{A} &=& \frac{1}{2\mu }\left(M_{11}-\frac{\lambda }{3K}M_{kk}\right), \end{eqnarray*}

\begin{eqnarray*} \Delta V_{B} &=& \frac{1}{2\mu }\left(M_{22}-\frac{\lambda }{3K}M_{kk}\right), \end{eqnarray*}

\begin{eqnarray} \Delta V_{C} &=& \frac{1}{2\mu }\left(M_{33}-\frac{\lambda }{3K}M_{kk}\right), \end{eqnarray}

(7)

where M_kk = M₁₁ + M₂₂ + M₃₃.

We find that eqs (3) and (7) establish a one-to-one relationship between any point ellipsoidal cavity and a CDM in the far field (see also Fig. 4). Moreover, it is straightforward to calculate the aspect ratios of the CDM from eq. (7). Using eqs (5) and (7) and the MATLAB function that we mentioned above, we also develop MATLAB functions for forward-model calculations associated with point ellipsoidal cavities.

The uniform pressure p on the walls of a pressurized ellipsoidal cavity with a volume of V is related to the actual volume change δV^E of the cavity through

\begin{eqnarray} \frac{\delta V^E}{V}-\epsilon _{kk}^T &=& -\frac{p}{K}, \end{eqnarray}

(8)

where |$\epsilon _{kk}^T$| is the sum of the diagonal components of the transformation strain (Segall 2010). By using the isotropic form of Hooke's law (Segall 2010), and assuming the bulk modulus to be the same inside and outside the cavity, it can be shown that

\begin{eqnarray} \epsilon _{kk}^T &=& \frac{\sigma _{kk}^T}{3K}, \end{eqnarray}

(9)

where |$\sigma _{kk}^{T}$| is the sum of the diagonal components of the transformation stress. Combining eqs (4), (5), (8) and (9) yields

\begin{eqnarray} \Delta V^E &=& \delta V^E + \frac{pV}{K}, \end{eqnarray}

(10)

where ΔV^E is the potency of the ellipsoidal cavity.

The potency and actual volume change of the ellipsoidal cavity, and consequently, its corresponding CDM, are related through eq. (10). It is straightforward to show that for the Mogi model, eq. (10) can be written as |$\Delta V^{\text{Mogi}} = \left(1+\frac{4\mu }{3K}\right) \delta V^{\text{Mogi}}$|⁠, which in a Poisson solid with μ = λ will be reduced to ΔV^Mogi = 1.8δV^Mogi (Bonafede & Ferrari 2009).

3.2 The CDM and finite ellipsoidal source comparison

To evaluate the efficiency of the CDM in a more general sense, we compare the surface displacement fields associated with finite ellipsoidal cavities and their corresponding CDMs, the results of which are presented in Figs 5 –7. For this purpose, we calculate the surface displacements associated with the ellipsoidal cavities using a numerical approach based on the boundary element method (BEM) after Kuriyama & Mizuta (1993). In the following comparisons, the moment tensor and depth to the centre of the corresponding ellipsoids and CDMs are identical. We calculate the moment tensors of the ellipsoidal cavities using eq. (5). Then, using eq. (7), we derive the potency and aspect ratios of the corresponding CDM uniquely. However, similar to other finite source models in the far field, there is a trade-off between the CDM dimensions and opening. Therefore, using the known potency and aspect ratios of the CDM from the previous step, we determine the optimized CDM dimensions and opening that provide the best fit to the BEM displacements. However, because the point CDM is dimensionless, the displacements associated with it can be directly calculated from the potency and aspect ratios. Following this approach for the comparisons minimizes the optimization-related uncertainties and guarantees that any deviation from the accurate BEM calculations is due to the mathematical and geometrical structures of the CDM and point CDM.

Figure 5.

Comparison of the BEM (solid line), CDM (dashed line) and point CDM (dotted line) normalized displacements along the X axis (left) and the Y axis (right) for the first ellipsoidal cavity with 0.75, 0.5 and 1.5 km semi-axes along the X, Y and Z axes, respectively.

Figure 6.

Comparison of the BEM (solid line), CDM (dashed line) and point CDM (dotted line) normalized displacements along the X axis (left) and the Y axis (right) for the second ellipsoidal cavity with 1, 2 and 0.75 km semi-axes along the X, Y and Z axes, respectively.

Figure 7.

Comparison of the BEM (solid line), CDM (dashed line) and point CDM (dotted line) normalized displacements along the X axis (left) and the Y axis (right) for the third ellipsoidal cavity with 2, 3 and 0.25 km semi-axes along the X, Y and Z axes, respectively.

We conduct the comparisons for various prolate, equidimensional and oblate cavities in the near field. Here, we only detail three cases, which represent all of the three geometrical categories of ellipsoidal cavities. However, utilizing triaxial ellipsoidal cavities and considering the model performance along the X and Y axes simultaneously, we cover a wider spectrum of cavity aspect ratios. In the following comparisons, the Lamé coefficients are μ = λ = 33 GPa, the depth to the centre of the cavity is d = 3 km, and the uniform pressure on the cavity walls is P = 10 MPa. We summarize the results of the comparisons in Figs 5–7 in which profiles of the BEM, CDM and point CDM horizontal and vertical surface displacements along the X and Y axes are presented. In each case, the distances are normalized by the depth to the centre of the cavity and the horizontal and vertical displacements are normalized by the maximum vertical displacement of each method.

For the first comparison, weuse an ellipsoidal cavity with semi-axes of 0.75, 0.5 and 1.5 km along the X, Y and Z axes, respectively. In this case, the point CDM and CDM displacements are equivalent, and both are in reasonable agreement with the BEM displacements (Fig. 5). The fit to the BEM displacements along the X and Y axes of the model are very similar; however, the quality of the fit for the horizontal displacements along the longest horizontal axis of the cavity, that is, along the X axis, is slightly better (Fig. 5, left). As the second case, we consider an ellipsoidal cavity with semi-axes of 1, 2 and 0.75 km along the X, Y and Z axes, respectively. In this comparison, the CDM displacements are in very good agreement with the BEM displacements, and in general, the CDM performs better than the point CDM (Fig. 6). In particular, along the major axis of the ellipsoid, which aligns with the Y axis, the CDM clearly shows a better fit to the BEM displacements than the point CDM does (Fig. 6, right). In the third comparison, we consider an ellipsoidal cavity with semi-axes of 2, 3 and 0.25 km along the X, Y and Z axes, respectively. Fig. 7 clearly shows that the CDM and BEM displacements match perfectly, whereas the point CDM overestimates the near-field displacements and underestimates the far-field displacements. This case shows that in the near field, the CDM can be used as an appropriate source model for rather thin magma bodies, particularly those that resemble planar intrusions such as dikes and sills.

The first ellipsoidal cavity in our comparisons resembles a prolate ellipsoidal source. For this cavity, the ratio of the semi-major axis to the depth to the centre, which we refer to as the characteristic ratio, is equal to 0.5. Prolate cavities in the near field can be approximated by a point source if their characteristic ratio is smaller than 0.5 (Yang et al. 1988). This explains the rather good agreement between the BEM displacements and those of the CDM and point CDM (Fig. 5). However, the equivalence of the CDM and point CDM displacements shows that for prolate ellipsoidal sources of pressurization, the point CDM can perfectly replace the CDM.

The second ellipsoidal cavity, however, represents a near equidimensional volumetric source. For a spherical cavity, the point-source approximation (Mogi 1958) is appropriate if the characteristic ratio of the cavity is smaller than 0.5 (McTigue 1987). The characteristic ratio of the second cavity in the XZ-plane is equal to 0.3, and as mentioned earlier, both the CDM and point CDM perform very good in this cross-section (Fig. 6, left). Additionally, the characteristic ratio of the cavity in the YZ plane is equal to 0.7. Therefore, in this plane, the point CDM does not perform as well as in the XZ plane. However, the CDM, due to its finite spatial extent, performs considerably better than the point CDM (Fig. 6, right).

The third ellipsoidal cavity resembles a horizontal sill. The point-source approximation of a horizontal penny-shaped crack (Fialko et al. 2001) is fully adequate for characteristic ratios smaller than 0.2 (Segall 2010). Therefore, in this case, the characteristic ratios of 0.7 and 1 completely explain the poor performance of the point CDM in the XZ and YZ planes, respectively (Fig. 7). However, again due to the finite spatial extent of the CDM as in the second test and also because of the aspect ratios of the third cavity, the performance of the CDM is excellent (Fig. 7).

We performed other tests using dipping cavities with arbitrary orientations in space. Similar to the three cases above, the results of these tests show that the performance of the CDM is a function of the depth and aspect ratios of the cavity and that the cavity orientation does not significantly alter the performance. According to the results of the evaluation tests that we discussed here, in the near field, the CDM and point CDM perform equivalently only for prolate ellipsoidal sources. However, in all the other tests, the performance of the CDM is better than the performance of the point CDM. In the following, we will apply the CDM to explain the deformation data measured at a real volcano.

4 THE 2015 CALBUCO CO-ERUPTIVE DEFORMATION MODELLING

After more than four decades of being dormant, Calbuco volcano (72.614°W, 41.326°S) in the southern volcanic zone in Chile erupted on 2015 April 22. The summit of the Calbuco volcano is approximately two kilometres above the mean sea level. This dome building volcano is an andesitic stratocone that produces pyroclastic and blocky lava flows (López-Escobar et al. 1995).

We estimate the surface displacements due to the 2015 eruption from the Sentinel-1A satellite ascending acquisitions of April 14th and April 26th and descending acquisitions of April 21st and May 3rd. We process the data following the Prats-Iraola et al. (2012) interferometric approach. For this purpose, we use the satellite state vectors and the 90 m resolution topographic data generated from NASA's Shuttle Radar Topography Mission (SRTM) to align the radar images and to remove the topographic component from the interferometric phase. By applying this method to the two Sentinel-1A pairs, we retrieve the two ascending and descending differential InSAR interferograms shown in Fig. 8. Due to the small perpendicular baselines of 43 m and 96 m for the Sentinel-1A ascending and descending track pairs, respectively, topographic errors would be minimal. Therefore, the phases observed in the differential interferograms are primarily attributed to deformation and minor atmospheric artefacts. However, in both cases, the topography correlated atmosphere has minimal effect on the Calbuco deformation because the main deformation fringes, which are coherent, fall in rather flat areas that have minor atmospheric noise (see Fig. 8).

Figure 8.

Sentine1-1A ascending (left) and descending (right) interferograms representing the ground deformation due to the 2015 April 22 eruption of the Calbuco volcano in Chile. The ascending images are acquired on April 14th and April 26th, and the descending images are acquired on April 21st and May 3rd.

The broad range of the observed surface deformation with a radius exceeding a few tens of kilometres about the summit implies that the deformation source is rather deep or horizontally elongated. In either case, considering the tests in Section 3, the CDM can be applied to the 2015 Calbuco co-eruptive displacements as a first source model.

We invert the subsampled Sentinel-1A ascending and descending line-of-sight (LOS) displacements (see Fig. A1) to infer the unknown parameters of the deformation source. For this purpose, we implement the CDM in a non-linear inversion scheme based on the genetic algorithm (GA; Haupt & Haupt 2004). The objective function that is minimized in the inversion is the L1-norm of the model residuals. We use the reciprocal of the observation variances as the observation weights in the inversion. Due to a better data quality and coherence, the standard deviation of the ascending data is on average half the standard deviation of the descending data. Therefore, the influence of the ascending data on the modelling outcomes is four times greater than the influence of the descending data. In addition to the unknown CDM parameters in the inversion, we also account for two more parameters, which represent the phase biases of the ascending and descending interferograms. We introduce these biases to the simulated ascending and descending LOS displacements as independent constant shifts. Consequently, these biases in the interferograms, which may exist due to atmospheric noise, unwrapping errors or local site effects, cannot significantly alter the modelling results. In the inversion procedure that we follow, the solution space is chosen to be very wide initially in order to span the likely parameter solution space and find the global minimum residual solution. We then restrict the solution space boundaries based on the preliminary inversion results. To estimate the uncertainties associated with the parameters of the optimal solution, we repeat the inversion 5000 times. We use a population size of N_pop = 60 chromosomes, a mutation rate of m_r = 0.25 and N_itr = 55 iterations in the GA (see Haupt & Haupt 2004). The Poisson ratio in the model is 0.25. We show the ascending and descending observations and predicted displacements and residuals from the best-fitting CDM in Fig. 9 (see also Fig. A1). The largest absolute value of the residuals for both the ascending and descending cases is about 1 cm. Additionally, the 1-D and 2-D marginal distributions of the estimated CDM parameters from the inversions, and also the geometry of the best-fitting CDM are illustrated in Fig. 10. The marginal distributions show that, in the restricted solution space, any trade-off between the estimated CDM parameters in the inversions are fairly well controlled. The best-fitting CDM parameters and their 95 per cent confidence bounds are summarized in Table 1. Using the estimated CDM dimensions and opening from Table 1 in ΔV = 4u(ab + ac + bc) (see Section 2.2), we estimate a potency of |$\widehat{\Delta V}=-0.099$| km³ for the best-fitting CDM. The 95 per cent confidence bounds for the estimated potency are (−0.114, −0.057). The minimum axis to maximum axis and intermediate axis to maximum axis ratios for the best-fitting CDM are |$\frac{a}{c}=0.171$| and |$\frac{b}{c}=0.693$|⁠, respectively.

Figure 9.

The results of the best-fitting CDM for the 2015 Calbuco eruption. The first row shows the InSAR line-of-sight (LOS) displacements, the CDM displacements and the residuals for the ascending acquisition. The second row corresponds to the descending acquisition. The stars in the residual panels indicate the centroid of the best-fitting CDM. The summit of Calbuco is located on the origin of the coordinate system.

Figure 10.

Results from the 5000 genetic algorithm inversions of the InSAR data using the compound dislocation model (CDM). (a) The 1-D and 2-D marginal distributions for estimated parameters of the CDM as well as the two phase bias terms (see Section 4) are shown on the diagonal and off-diagonal panels, respectively. The black distributions represent the parameters that the CDM and point CDM have in common. The grey distributions represent the parameters that determine the source shape and potency. Vertical red lines and black points indicate parameter values for the best-fitting CDM, respectively. The exact numerical values of the estimated parameters and their 95 per cent confidence bounds are given in Table 1. (b) Illustration of the south-north and east-west cross-sections and the top view of the best-fitting CDM in the local coordinate system as introduced in Fig. 9. Geodetic coordinates are provided in Table 1. The colours and labelling of the RDs are the same as in Fig. 3.

Table 1.

The parameter values and 95 per cent confidence bounds of the best-fitting CDM associated with the 2015 Calbuco eruption. The λ₀ and ϕ₀ are the geodetic longitude and latitude of the CDM centroid, respectively, and d is the depth to the CDM centroid from the free surface in the model. The ω_X, ω_Y and ω_Z are the angles of rotation about the X, Y and Z axes, respectively. The a, b and c are the semi-axes of the CDM that, before applying the rotations, align with the X, Y and Z axes, respectively. The bias_asc and bias_dsc are the ascending and descending bias terms. The coordinates of the best-fitting CDM centroid in the local coordinate system in Fig. 9 are (0.263, −0.161) km. The 95 per cent confidence bounds for these local coordinates from the inversions are (−2.445, 1.337) and (−2.710, 1.070), respectively.

Parameter	Value	95 per cent confidence bounds	Parameter	Value	95 per cent confidence bounds
λ0 (°)	−72.611	(−72.643, −72.598)	a (km)	0.409	(0.200, 0.731)
ϕ0 (°)	−41.327	(−41.350, −41.316)	b (km)	1.660	(0.936, 1.844)
d (km)	8.206	(6.557, 9.534)	c (km)	2.395	(1.805, 3.305)
ωX (°)	1.316	(−12.494, 12.669)	u (m)	−4.398	(−6.307, −2.591)
ωY (°)	−4.023	(−8.653, 11.631)	biasasc (cm)	−0.100	(−0.974, 0.340)
ωZ (°)	159.500	(144.886, 178.928)	biasdsc (cm)	1.163	(0.542, 1.698)

In our second modelling attempt, we utilize the point CDM as the source model and follow the same inversion procedure that we mentioned above. The predicted displacements and residuals from the best-fitting point CDM are almost the same as those of the CDM (see Fig. A2). The marginal distributions of the estimated point CDM parameters from the inversions, and also the force dipoles that form the best-fitting point CDM are depicted in Fig. 11. The corresponding estimated parameters and 95 per cent confidence bounds are summarized in Table 2. We use the estimated potencies for the individual point dislocations in ΔV = ΔV_A + ΔV_B + ΔV_C and calculate the total point CDM potency of |$\widehat{\Delta V}=-0.091$| km³ with (−0.176, −0.087) as its 95 per cent confidence bounds. We also calculate |$\frac{da}{dc}=\frac{\Delta V_{C}}{\Delta V_{A}}=0.151$| and |$\frac{db}{dc}=\frac{\Delta V_{C}}{\Delta V_{B}}=0.706$|⁠, where the da, db and dc are the differential axes of the point CDM (see Section 2.2), and |$\frac{da}{dc}$| and |$\frac{db}{dc}$| are the ratios of minor axis to major axis and intermediate axis to major axis, respectively.

Figure 11.

Results from the 5000 genetic algorithm inversions of the InSAR data using the point CDM. (a) The 1-D and 2-D marginal distributions for estimated parameters of the point CDM and the two phase bias terms (see Section 4). The grey distributions represent parameters that determine the source potency. Vertical red lines and black circles indicate parameter values for the best-fitting point CDM. The exact numerical values of the estimated parameters and their 95 per cent confidence bounds are given in Table 2. (b) Cross-sections and the top view of force dipoles of the best-fitting point CDM. The dipole magnitudes are calculated from eq. (3). The orientation of each dipole shows the opening direction of a point dislocation. The colour and labelling allow for comparison to the RDs in Figs 3 and 10. Note that the dipole with the largest magnitude corresponds to the dislocation with the largest potency, which in case of a CDM is the RD with the largest area.

Table 2.

The parameter values and 95 per cent confidence bounds of the best-fitting point CDM associated with the 2015 Calbuco eruption. The ΔV_A, ΔV_B and ΔV_C are the potencies of the point dislocations that, before applying the rotations, are normal to the X, Y and Z axes, respectively. The other parameters are the same as the CDM (see Table 1). The coordinates of the best-fitting point CDM in the local coordinate system in Fig. 9 are (0.521, −0.607) km. The 95 per cent confidence bounds for these local coordinates from the inversions are (−2.445, 1.337) and (−2.710, 1.070), respectively.

Parameter	Value	95 per cent confidence bounds	Parameter	Value	95 per cent confidence bounds
λ0 (°)	−72.608	(−72.630, −72.576)	ΔVA (km3)	− 0.067	(−0.124, −0.055)
ϕ0 (°)	−41.332	(−41.351, −41.309)	ΔVB (km3)	− 0.014	(−0.040, −0.005)
d (km)	8.271	(7.732, 10.953)	ΔVC (km3)	− 0.010	(−0.030, −0.007)
ωX (°)	2.086	(−12.632, 13.110)	biasasc (cm)	0.009	(−0.383, 1.087)
ωY (°)	1.115	(−11.732, 8.229)	biasdsc (cm)	1.119	(0.340, 1.808)
ωZ (°)	160.259	(144.930, 178.616)

4.1 Implications of the 2015 Calbuco eruption modelling

The inversions of the 2015 Calbuco co-eruptive displacements that we detailed in the previous section, do not rely on any a priori information or constraints on the deformation source geometry and mechanism. This is a major advantage, particularly because no prior knowledge was available regarding the depth and type of the magma source. The parameters of the best-fitting CDM in Table 1 clearly show that the deformation source is a rather thin, vertically elongated source that may somewhat resemble a dike. From the best-fitting CDM parameters, a characteristic ratio of 0.3 can be estimated for the deformation source. Given this small characteristic ratio, one can reasonably deduce that the far-field approximation can be applied to the deformation source, implying that the point CDM can replace the CDM in this case. We tested this by repeating the inversion procedure using the point CDM as the source model. As expected, the depth and potency, as well as the aspect ratios and the surface displacement field associated with the best-fitting point CDM, are equivalent to those of the best-fitting CDM (see Section 4).

The potency of the largest RD of the best-fitting CDM, however, is only 71 per cent of the total CDM potency. This result shows that the deviation of the deformation source from a planar tensile dislocation is significant. The large differences between the estimated aspect ratios of the optimal CDM in Section 4 show that the deformation source is also different from the McTigue (1987) and Mogi (1958) isotropic source models. Using eq. (2), we calculated M₁, M₂ and M₃, which are the eigenvalues of the optimal CDM moment tensor. The eigenvalue ratios M₂/M₁ = 0.559 and M₃/M₁ = 0.515 correspond to a point near the lower boundary of the domain of ellipsoids in Fig. 4, again showing that the deformation source has a lens-shaped geometry, different from a plane and substantially different from a sphere. To check this result numerically, we performed other inversions using the Mogi (1958) and Okada (1985) models (see Appendix A). The best-fitting Mogi model is located at 4.5 km to the west and 2.9 km to the south of the volcano at a depth of 8.3 km, and its potency is −0.058 km³. The estimated depth is comparable to that of the best-fitting CDM in Section 4; however, the potency is approximately 41 per cent smaller. Although the model fit to the ascending data in this case is good, the Mogi model completely fails to model the descending data (see Fig. A4). The centre of the best-fitting Okada (1985) model, however, is located at 1.5 km to the east and 0.4 km to the north of the volcano at a depth of 8.1 km. The strike and dip angles of the best-fitting Okada (1985) model are 157° and 88°, respectively. The horizontal and along-dip dimensions of this best-fitting RD are 3.8 km and 2.7 km, and its potency is −0.066 km³. In this case, the depth to the centre also compares well to the CDM depth, but the potency is approximately 33 per cent smaller than the best-fitting CDM potency. As expected from the moment tensor ratios of the optimal CDM, the Okada (1985) model performs better than the Mogi (1958) model. However, the Okada (1985) model residuals are greater than those of the best-fitting CDM (see Fig. A3), and given that the best-fitting Okada (1985) model is almost vertical, the deviation of its centre from the volcano centre cannot be realistic.

Considering the results of the 2015 Calbuco co-eruptive displacement modelling in Figs 9 and A1, the CDM displacements show a very good agreement with both sets of the ascending and descending observations with residuals that are generally smaller than 1 cm. The residuals may partly be due to atmospheric artefacts in the InSAR data. The small systematic residuals near the summit area, however, can be related to the shallow part of the plumbing system, which became active during the eruption. Given the very small magnitude of the residuals, we speculate that any contribution from other shallow processes cannot be significant. We note that each InSAR interferogram reflects the total displacements over the period of time spanned by the two master and slave images. However, these periods for the ascending and descending acquisitions are not the same (see Section 4). The possible differences between the temporal variations of the displacements in the ascending and descending interferograms might be reflected in the results as the systematic residuals. Because both ascending and descending InSAR data can be well explained by a single model, we conjecture that the deformation episode was fully covered by the data.

5 DISCUSSION

The CDM that we presented in this paper is a generalization of the Bonafede & Ferrari (2009) dislocation model. Under the far-field approximation, the CDM represents a pressurized ‘rectangular box’ of arbitrary aspect ratios and orientation (Bonafede & Ferrari 2009; Ferrari et al. 2015), which in some cases is equivalent to the Davis (1986) point ellipsoidal cavity in a half-space (see Fig. 4). The aspect ratios and the pressure boundary condition of the latter case can be determined through the Eshelby (1957) theory. However, in the near-field case, the CDM can represent finite planar intrusions, namely, dikes and sills, and more volumetric variations thereof. In either case, given the simple configuration of the CDM as a kinematic model (see Section 2.2), the potency associated with the underlying deformation source can be easily estimated through simple arithmetic operations, which utilize the CDM dimensions and opening as well as the Lamé parameters. Central to the implementation of the CDM are the new solutions for the displacement and stress fields of the RDs that we developed in Section 2.1.

5.1 The limitations of the CDM

Like all other analytical solutions used for volcano deformation modelling, the CDM assumes that the Earth's crust is an isotropic, uniform elastic half-space. Therefore, to apply the CDM to case studies that involve a rough topography or substantial layering and material heterogeneity, one would need to be very careful about possible contributions of these factors. Moreover, the CDM can only represent the geometry of the deformation that a subsurface fluid reservoir may undergo. In the absence of further physical and geological evidences the inferred deformation geometry does not necessarily represent the real shape of the reservoir.

Regarding the results of the evaluation tests in the near field (see Section 3), for pressurized cigar-shaped cavities, the CDM and the point CDM perform alike. However, for the other cases in the near field, the more oblate the cavities are, the better the CDM performs over the point source models. This is very distinctive particularly for sheet-like intrusions such as dikes and sills in the near field (see Fig. 7). Considering the geometrical structure of the CDM in Fig. 3, it is clear that in the vicinity of pressurized cavities, the CDM cannot accurately simulate the stress field. In such cases and for more detailed modelling of the displacement and stress fields associated with the deformation sources, numerical models provide considerably better results.

In contrast to the source potency, the pressure that a deformation source applies to the surrounding medium does not appear in the CDM formulation. Therefore, the problems that require pressure boundary conditions associated with the cavities cannot be investigated through applying the CDM only. In such cases and for estimating the pressure associated with deformation sources, the CDM relies on the Eshelby (1957) solution. However, because of this simple geometrical and mathematical structure, the CDM and point CDM are versatile models for kinematic source modelling purposes.

5.2 Method improvements and comparison to earlier works

The Okada (1985, 1992) solutions for RDs in a half-space have the limitation that the upper edge of the RDs is constrained to be parallel to the free surface. We addressed this issue by adding the ‘plunge angle’ to the new RD solutions, which have full rotational degrees of freedom (see Figs 2 and 3b). Note that such an RD solution can also be formed as a superposition of two triangular dislocations (TDs), as detailed in Nikkhoo & Walter (2015). However, such an approach is not only computationally more expensive but also due to the diagonal dislocation line in its configuration, it does not allow for calculations at the centroid of the RD. The latter in particular is an important drawback for implementing the RD solutions in computer routines based on the boundary element method (BEM).

Another problem in the Okada (1985, 1992) solutions is the imperfection in dealing with the artefact singularities and numerical instabilities along the edges and underneath and above the vertices of the RDs. For the calculation points within a predefined neighbourhood of the artefact singularities, Okada (1992) replaced the singular and numerically instable terms by their mathematical limits. This approximative approach, however, makes the Okada (1992) solution scale dependent, meaning that depending on how small the RDs in a model are or how close the calculation points to the artefact singularities are, the accuracy of the results may be degraded. Bradley & Segall (2012) considered this problem in a high-resolution BEM simulation based on RDs and addressed the instabilities along the RD edges. However, the Bradley & Segall (2012) work does not remark on the other artefacts and instabilities associated with the Okada (1992) solution. Following the Nikkhoo & Walter (2015) approach, we addressed the problem of the artefact singularities in the RD solutions in both full-space and half-space (Fig. 1).

Bonafede & Ferrari (2009) showed the analogy between the Mogi model and a dislocation model composed of three mutually orthogonal square tensile dislocations. They interpreted the dislocation model as a pressurized box in the far field and detailed its implications in the internal and external volume changes associated with a magma chamber. A generalization of the Bonafede & Ferrari (2009) model in its primary form that is based on the Okada (1985, 1992) solutions could only simulate pressurized horizontal and vertical ‘boxes’. This is because of the geometrical limitation of the Okada (1985, 1992) RD solutions. Thanks to the new RD solutions that we discussed above, the CDM can take any arbitrary orientation in the space.

Ferrari et al. (2015) developed a numerical model for pressurized rectangular ‘boxes’ by using the Okada (1992) solution in a BEM scheme. Each side of the rectangular box in this model can be composed of any number of RDs. A uniform traction boundary condition is fulfilled at the centroid of the RDs that form the box. This numerical model is a good approximation of rectangular boxes in both the far field and near field. The moment tensor representation of the Ferrari et al. (2015) model and the CDM are identical, and their corresponding eigenvalue ratios cover the same area in the moment tensor spectrum (see Fig. 4). Nevertheless, due to the same geometrical limitation in the Okada (1992) RD solution, the Ferrari et al. (2015) model can only simulate pressurized horizontal and vertical rectangular ‘boxes’. Clearly, this issue can be addressed by using the RD solutions that we developed in Section 2.1.

Under the far-field approximation, a second-order moment tensor that corresponds to a system of force couples located at a point can represent generalized point sources, including uniformly pressurized cavities (Aki & Richards 2002). Although the inversion of geodetic data sets for the location and components of a moment tensor does not require any a priori information about the geometry of the deformation source, the interpretation of the source geometry and mechanism through the inferred moment tensor is not unique (Trasatti et al. 2009).

The point CDM, as a special case of a moment tensor (see Sections 1 and 3.1), can only represent volumetric changes due to expansion or contraction in all directions in space. This is because in the point CDM formulation, potencies of the point dislocations that form the point CDM are constrained to have the same sign. As a result, the point CDM does not cover the whole domain of universal moment tensors as illustrated in Fig. 4. Including arbitrary positive and negative potencies for the point dislocations in the point CDM formulation allows for modelling of sources that randomly expand or contract in any of the three main directions. However, in the present paper we do not consider this case.

Based on the moment tensor inversion concept, Davis (1986) developed a multi-stage inversion scheme for estimating the unknown parameters of a generalized point ellipsoidal cavity in the far field. After inverting the surface displacements for the components of a generic moment tensor, Davis (1986) estimated the source parameters by using the moment tensor eigenvalues and eigenvectors as well as the Eshelby (1957) tensor. However, the estimated moment tensor in the first step of the Davis (1986) method might not represent a pressurized ellipsoidal cavity (Trasatti et al. 2009). Consequently, for some cases, the Davis (1986) inversion scheme does not converge to an optimal solution. As another disadvantage, the Davis (1986) method, because of its special multi-stage inversion scheme, cannot be implemented in computer codes that use forward-model-based heuristic search algorithms, such as genetic algorithm (GA) or Markov chain Monte Carlo (MCMC) methods. These problems have kept the Davis (1986) generalized point ellipsoidal cavity model from replacing its isotropic counterpart, that is, the Mogi model. Considering the details in Section 2.2, it is clear that these issues do not apply to the CDM. Using the relationship between the point CDM and point ellipsoidal cavity in Section 3, one can easily perform the forward modelling of the point ellipsoidal cavity. In this way, the convergence issue can be avoided. The point CDM covers the entire domain of pressurized cavities in the far field (see Fig. 4), while its application in practice is as convenient as the Mogi model. The latter is because the point CDM does not rely on the Eshelby (1957) method for calculating the displacement and stress fields and the source volume change. Moreover, similar to other dislocation models, the point CDM directly allows for the forward model calculations, and it is straightforward to integrate this source model in various inversion routines. Therefore, the point CDM is an appropriate alternative for generalized point-source cavity models in the far field.

In addition to the planar and ellipsoidal deformation sources, the CDM can also be used to represent other source geometries. Bonaccorso & Davis (1999) presented approximate analytical models for closed and open vertical volcanic conduits or pipes. The closed pipe model is the limiting case of a pressurized prolate ellipsoid with a negligible short axis. The open pipe model is the limiting case of a tensile cylindrical dislocation with a negligible radius. Both models are effectively uniform distributions of moment tensors along their vertical axis. The moment tensors in the closed pipe model generate lateral and upward pressure, but in the open pipe model, they only exert lateral pressure to the surrounding medium. A narrow vertical CDM has a similar configuration to the open pipe model. However, the CDM generally cannot produce the same amount of upward push that exists in the closed pipe model. The characteristic ratio of the closed pipe in Fig. 12 is 0.5, but the CDM approximation is not perfect. Therefore, the CDM is not an ideal substitute for the closed pipe model, except for considerably smaller characteristic ratios, which may not be realistic. However, the limiting case of a very narrow CDM, which coincides with the pipe's vertical axis, can perfectly simulate the open pipe model (see Fig. 12). Moreover, the CDM in this case can be applied to more general non-vertical conduit models with a non-circular base.

Figure 12.

Comparison of the CDM and the Bonaccorso & Davis (1999) pipe models. In both cases, the conduit is a vertical cylinder with a height of 1 km, a radius of 0.1 km and a centroid depth of 1.5 km. The horizontal and vertical surface displacements of the closed pipe model (left) and the open pipe model (right) are normalized by the maximum vertical displacement and maximum horizontal displacement, respectively. The solid and dashed lines represent the pipe and the CDM, respectively.

In the near field, the finite size of a deformation source must be included in the source model, and rather than a point-source model, a finite source or a multipole point-source model should be used (Davis 1986). In fact, the explicit equations in the Eshelby (1957) theory first showed that the displacement field of a pressurized ellipsoidal cavity in a full-space is exactly equivalent to that of a uniform distribution of three orthogonal force dipoles inside the ellipsoid (Segall 2010). The full-space Green's functions (Kelvin 1848) associated with these force dipoles can be replaced by the half-space Green's functions (Mindlin 1936) to achieve an approximate solution in the half-space (Segall 2010). These latter points better explain the performance of the CDM in modelling the finite cavities in the near field that we discussed in Section 3.2. Each of the RDs that form the CDM has a uniform distribution of orthogonal force dipoles with a moment tensor proportional to those in eq. (1). Therefore, in contrast to point-source models, the CDM can represent a portion of the uniform distribution of some force dipoles in a cavity. This property of the CDM can generally improve the performance of the CDM compared to the point-source models (see Fig. 6). However, this improvement for planar intrusions and thin oblate cavities in the near field can be significant (see Fig. 7). Clearly, when the length of one of the CDM edges tends to zero, the CDM will be reduced to a single RD (see Fig. 3). Therefore, similar to finite dislocation models that in many cases are good source models for dikes and sills (Davis 1983; Segall 2010), the CDM can also adequately model the planar intrusions in the near field.

6 CONCLUSIONS

In this paper, we developed analytical artefact-free solutions for RDs in the full-space and half-space. These RD solutions possess full rotational degrees of freedom, meaning that their upper edge is not constrained to be parallel to the free surface. Using these RD solutions, we constructed a CDM in terms of three mutually orthogonal tensile RDs. Under the far-field approximation, that is, for small source dimension to depth ratios, the CDM can represent any triaxial ellipsoidal or planar sources. The point CDM is a generalized point source that covers a wider spectrum of moment tensors compared to the Davis (1986) point ellipsoidal cavity. In the near field, however, the CDM can simulate planar intrusions, namely, dikes and sills, as well as thin oblate ellipsoidal sources. The potency of the CDM is the sum of the potencies of the individual RDs that it comprises. Regardless of the shape of the deformation sources, their potency can be uniquely determined by using the CDM.

The convenience of integrating the CDM and point CDM in various optimization schemes compares to the Mogi model as the simplest analytical point-source model. Therefore, the CDM can be used as a first model for rapid co-eruptive displacement modelling. This is particularly of great relevance for early warning and rapid response systems, which shortly after or during a crisis require the volcanic source parameters without any preliminary information on the shape and depth of the source. As the first option in these cases, using the point CDM can provide an even higher performance. If this model could not succeed in simulating the ground deformation, it implies that most probably a finite source in the near field is causing the deformation. As the second option, the CDM can then be used.

We applied the CDM to the 2015 Calbuco co-eruptive ground deformation observed by the Sentinel-1A InSAR satellite. The inferred deformation source is a vertically extended lens-shaped magma body at a depth of 8.2 km. The estimated potency of the source is −0.1 km³. We showed that applying the isotropic Mogi (1958) model or the Okada (1985) model to this case study produces erroneous estimations of the deformation source parameters.

We provide the 2015 Calbuco co-eruptive subsampled InSAR data. The MATLAB functions that accompany this paper allow for easy calculations using the new RD solutions, as well as the CDM and the point CDM.

Acknowledgments

We appreciate constructive reviews and comments from Sylvain Barbot, an anonymous reviewer and the editor. This is a contribution to VOLCAPSE, a research project funded by the European Research Council under the European Union's H2020 Programme/ERC consolidator grant No. [ERC-CoG 646858]. We acknowledge the financial support from ERC project ‘TOPOMOD’ (contract No. 264517), as well as from the InSARAP project (ESA Contract 4000110587/14/I-BG). The subsampled InSAR data and the MATLAB functions associated with the RD and CDM calculations in this paper are available upon request from the corresponding author and can also be downloaded under the following link: www.volcanodeformation.com.

APPENDIX A: COMPARISON OF RESULTS FROM INVERSIONS OF THE 2015 CALBUCO CO-ERUPTIVE DISPLACEMENTS USING VARIOUS SOURCE MODELS

The Mogi (1958) and Okada (1985) models are the most well-known and widely used analytical models in volcano deformation studies. Due to the simplicity of the mathematical formulation of the Mogi (1958) model in particular, using it in practice is very convenient. However, as also shown by Davis (1986), the application of the Mogi model may result in large errors in some of the estimated source parameters. To determine how this is relevant in the 2015 Calbuco co-eruptive deformation, in addition to the CDM and point CDM, we also perform separate inversions using the Mogi model and Okada (1985) model, respectively. In the following figures, we show the subsampled InSAR data that we used in the inversions. We use the local XY coordinate system that we introduced in Fig. 9. We show the results associated with the CDM in Fig. A1. In this case study, the CDM results are almost equivalent to those of the point CDM that are shown in Fig. A2. We also show the results of the modelling for the Okada (1985) model and the Mogi model in Figs A3 and A4, respectively. In this case, the Mogi model can simulate the ascending data very well, but it completely fails to simulate the descending data (see Fig. A4). The inferred deformation source derived using the Mogi model is 5.4 km off-centred (see Section 4.1). The deviation of the best-fitting Okada (1985) model from the Calbuco centre is 1.5 km. The Okada (1985) model simulates the general pattern of the ascending and descending displacements better than the Mogi model. Nevertheless, neither of these two models performs as well as the CDM and point CDM. This can be better seen in the root mean square error (RMSE) values of the residuals of the four models in Table A1. It should be noted that the Mogi model can be considered as the special case of a CDM that has three equal axes. Also, the Okada (1985) model is the special case of a CDM that has only two non-zero axes. Therefore, in general the CDM performs better than the Mogi (1958) and Okada (1985) models and in special cases performs the same as these two models.

Figure A1.

The results of the CDM for the 2015 Calbuco eruption. The first row shows the subsampled InSAR observations, the CDM displacements and the residuals for the ascending acquisition. The second row corresponds to the descending interferogram. The summit of Calbuco is located on the origin of the coordinate system. The stars in the residual panels indicate the centroid of the best fitting CDM.

Figure A2.

The results of the point CDM for the 2015 Calbuco eruption. The first row shows the subsampled InSAR observations, the point CDM displacements and the residuals for the ascending acquisition. The second row corresponds to the descending interferogram. The summit of Calbuco is located on the origin of the coordinate system. The stars in the residual panels indicate the centroid of the best fitting point CDM.

Figure A3.

The results of the Okada (1985) model for the 2015 Calbuco eruption. The first row shows the subsampled InSAR observations, the Okada (1985) model displacements and the residuals for the ascending acquisition. The second row corresponds to the descending interferogram. The summit of Calbuco is located on the origin of the coordinate system. The stars in the residual panels indicate the centroid of the best-fitting Okada (1985) model.

Figure A4.

The results of the Mogi model for the 2015 Calbuco eruption. The first row shows the subsampled InSAR observations, the Mogi model displacements and the residuals for the ascending acquisition. The second row corresponds to the descending interferogram. The summit of Calbuco is located on the origin of the coordinate system. The stars in the residual panels indicate the location of the best-fitting Mogi model.

Table A1.

Quantitative comparison of the residuals of the CDM, the point CDM, the Okada (1985) model and the Mogi model for the 2015 Calbuco eruption. The min, max and RMSE in the table stand for the minimum residual, maximum residual and the root mean square error of the residuals, respectively. All the numerical values are given in (mm).

	Asc. data			Desc. data			Asc. & Desc. data
Model	Min	Max	RMSE	Min	Max	RMSE	Min	Max	RMSE
CDM	− 9.8	9.2	3.5	− 12.5	12.4	5.2	− 12.5	12.4	3.9
Point CDM	− 10.9	8.7	3.6	− 11.9	13.3	5.4	− 11.9	13.3	4.0
Okada (1985)	− 19.2	10.9	4.6	− 15.4	12.9	4.4	− 19.2	12.9	4.6
Mogi (1958)	− 13.2	18.6	5.5	− 14.2	53.4	15.2	− 14.2	53.4	8.2