Analytical models of volcanic ellipsoidal expansion sources

Modeling non-double-couple earthquakes and surficial deformation in volcanic and geothermal areas usually involves expansion sources. Given an ensemble of ellipsoidal or tensile expansion sources and double-couple ones, it is straightforward to obtain the equivalent single moment tensor under the far-field approximation. On the contrary, the moment tensor interpretation is by no means unique or unambiguous. If the far-field approximation is unsatisfied, the single moment tensor representation is inappropriate. Here we focus on the volume change estimate in the case of single sources, in particular finite pressurized ellipsoidal sources, presenting the expressions for the computation of the volume change and surficial displacement in a closed analytical form. We discuss the implications of different domains of the moment-tensor eigenvalue ratios in terms of volume change computation. We also discuss how the volume change of each source can be obtained from the isotropic component of the total moment tensor, in few cases of coupled sources where the total volume change is null. The new expressions for the computation of the volume change and surficial displacement in case of finite pressurized ellipsoidal sources should make their use easier with respect to the already published formulations.


Introduction
Modeling non-double-couple (non-DC) earthquakes and surficial deformation in volcanic and geothermal areas usually involves expansion sources, e.g.microcracks, magma chambers, and dikes.
Observations of non-DC earthquakes include earthquakes related to landslides and volcanic eruptions, long-period volcanic earthquakes, short-period volcanic and geothermal earthquakes, earthquakes at mines, deep-focus earthquakes, and other shallow earthquakes [for a review, see e.g.Miller et al. 1998].As regards shallow earthquakes in volcanic or geothermal areas and mines, isotropic components are consistent with failure involving both shear and tensile faulting, which may be facilitated by high-pressure, high-temperature fluids, and with cavity closing in mines.Also microearthquakes induced by hydraulic fracturing [e.g.Sileny et al. 2009] may have significant non-DC parts including a positive isotropic component consistent with opening of cracks oriented close to expected hydraulic fracture orientation, and could be the result of crack opening by fluid injection.For mixed mode earthquakes, which combine shear and tensile dislocations on a planar fault that can possibly be opened or closed during the rupture process, the deviation of the dislocation vector from the fault, i.e. the tensility (or slope) of the source, can be determined uniquely [Vavryčuk 2011]; volume change of the source ΔV is given by m/(m+2n/3), where m is the isotropic (ISO) component of the moment tensor (one third its trace), and m and n are the Lamé parameters.
Seismic moment tensors are not always consistent with combined shear and tensile faulting [see e.g.Foulger et al. 2004].As regards source mechanisms of deep low-frequency earthquakes, the isotropic component is often assumed to represent a spherical magma oscillation [e.g.Julian et al. 1998, Miller et al. 1998], and in this case ΔV is given by m/(m+2n).Very-long-period earthquakes at volcanoes have sometimes been interpreted in terms of a mutual deflation and inflation of two connected magma chambers [e.g.Nishimura et al. 2000] whose shape is not necessarily spherical.
It is common practice to decompose the moment tensor into ISO and trace-free parts, and obtain ΔV from the isotropic component only.On the other hand Wielandt [2003] has demonstrated that ΔV cannot be obtained from ISO without choosing a priori a source model.There is no general relationship between ISO (m) and ΔV, but only between ISO and the volume displaced through any spherical surface enclosing the source (i.e. the surface integral of the normal displacement over the spherical surface), which is given by m/(m+2n).This result is true for any position of the point moment tensor within the given sphere and thus remains valid for distributed sources.For a given ISO, a source model with any volume change, positive or negative, can be hypothesized.Only for a single expansion source (e.g., an ellipsoidal cavity) a positive ISO can only result in a positive ΔV.It also follows that the surface integral of the normal displacement over a surface enclosing the source cannot be used to estimate ΔV, but in case the enclosing surface is exactly the boundary of a single expansion source, like a pressurized cavity or a tensile fault.Unfortunately, displacements on the source boundary computed using the point moment tensor approximation of the source are different from the actual ones because the far-field assumption (implicit in the point moment tensor approximation) is obviously unsatisfied.These fundamental results are often not considered in research work.
Expansion sources are usually modeled as ellipsoids and tensile cracks.Amoruso and Crescentini [2009] have computed exact expressions for the volume change of a pressurized ellipsoidal cavity in an infinite homogeneous elastic medium that are approximate solutions for a homogeneous half-space.Amoruso and Crescentini [2009] have also shown that the moment tensor of earthquakes generated by a sudden magma exchange between two ellipsoidal cavities (e.g., because of the breaking of a barrier) can include all the three (ISO, DC, Compensated Linear Vertical Dipole -CLVD) components, and the volume change obtained from the isotropic component only can be much smaller than that really involved in the magma exchange process.From the analysis of the CLVD and ISO components, it is possible to get clues on the source shapes and estimate their volume change.
As regards surficial deformation, its analysis is a widely used tool for studying volcanoes, and in particular the critical stages prior to eruptions.The source is often modeled as a single pressurized cavity [e.g.Amoruso et al. 2007, Amoruso et al. 2008], and the density of the material entering the source may be inferred from gravity data, if the source geometry and volume change are known; this approach may allow discriminating between deformation episodes due to the hydrothermal system and those ascribable to magma intrusion.If the source dimensions are small with respect to depth, surficial deformation is the same as from a single (point) moment tensor, and what previously outlined for non-DC earthquakes still holds.If the source dimensions are not small with respect to depth, a further complication arises because of the small number of analytical solutions for finite sources [for a review, see e.g.Lisowski 2006].The point moment tensor representation is equivalent to the monopole term in the multipole expansion of any extended deformation source, but the resulting far-field approximation may be inappropriate.
Presently the library of available analytical (approximate or exact) solutions for surface displacements includes few finite cavity shapes [spheres, McTigue 1987;prolate spheroids, Yang et al. 1988;fluid-pressurized closed pipes or dislocating open pipes, Bonaccorso and Davis 1999;circular horizontal cracks, Fialko et al. 2001; finite triaxial ellipsoids, Amoruso and Crescentini 2011] and very small (with respect to depth) cavities [spheres, Mogi 1958; generic triaxial ellipsoids, through its representation in terms of a single moment tensor, Davis 1986].Assuming homogenous elastic conditions, general geometrical configuration of pressure sources can be obtained from surficial displacement data [e.g.Vasco et al. 2002, Camacho et al. 2011].These sources are described as an aggregate of pressure point sources, which are able to model hydrothermal systems [e.g.Rinaldi et al. 2010] but not expanding non-spherical cavities or aggregates of cracks.
In this work we give closed analytical expressions for computing the volume change of finite pressurized triaxial ellipsoidal sources.We detail that, if a single source is assumed and the moment tensor eigenvalues are consistent with an ellipsoidal source or a mixed mode dislocation, the volume change is proportional to the moment tensor trace, but the proportionality constant depends on the source shape.We also discuss three examples of coupled cavities.
At last, we give explicit analytical expressions for computing surficial displacements from a pressurized ellipsoidal cavity, under the quadrupole approximation.

Estimation of the volume change of expansion sources
In this section, we summarize some results about the computation of ΔV and the moment tensor representation of pressurized ellipsoidal cavities, and give explicit expressions some of which still unpublished.Later on, we discuss how the volume change of each source can be obtained from the ISO component of the total moment tensor, in few cases of coupled sources where the total volume change is null.

Ellipsoidal cavity
The elastic field due to an ellipsoidal inclusion in a homogeneous half-space has been treated by Davis [1986], following Eshelby's [1957] approach but using the point force solution for a half-space rather than a full-space (Kelvin force) as the fundamental Green's function.The solution satisfies boundary conditions exactly on the free surface but approximately on the ellipsoid, and is reasonably accurate if the depth to the ellipsoid center is larger than twice its dimension.Davis' far-field solution is the sum of displacements from three co-located double forces acting along the axes of the ellipsoid.Far-field deformation from a pressurized ellipsoidal cavity is thus the same as from a moment tensor, whose eigenvalues are proportional to the product of the source pressure and volume PV and depend on the axis ratios; the normalized eigenvectors of the moment tensor represent the directions of the el-lipsoid axes.Although ΔV is of great importance, e.g., for estimating the intrusion density from gravity data, its link to PV is not trivial.If the reference axes are aligned with the ellipsoid axes (a, b, c), Davis [1986] showed that far-field deformation from a pressurized ellipsoidal cavity is the same as deformation from the moment tensor and gave expressions for computing P a T /P, P b T /P, P c T /P as a function of the ellipsoid axis ratios, in terms of ellip-tic integrals.Amoruso and Crescentini [2009] have shown that for an ellipsoidal inclusion where o is the Poisson's ratio, p T = P a T + P b T + P c T is rotational invariant, being proportional to the trace, and can be calculated for any orientation of the cavity axes, and P is positive for overpressure.In terms of ΔV, the moment tensor can be written as Thus, ΔV is proportional to the ISO component (m) of the moment tensor that represents the far-field displacements from the ellipsoidal source, but the pro-portionality factor depends on the shape of the ellipsoidal cavity: For a given m, ΔV can vary up to a factor of about two, depending on the moment tensor eigenvalue ratios.
If two of the axes of the ellipsoid are equally long (i.e., the ellipsoid is a spheroid), the eigenvalues of the moment tensor can be computed in terms of elementary functions (for details, see the work of Eshelby [1957]).Amoruso and Crescentini [2009] obtained approximate results for very small or very large axis ratio using appropriate Taylor expansions of equations (3.15) and (3.16) in the work by Eshelby [1957].
Here we give the details of the computations, in particular explicit expressions for P a T /P, P b T /P, P c T /P in terms of elliptic integrals and in algebraic form for a couple of limit cases.These expressions are not new, but are given in a form that might make their use easier with respect to previously published ones.We stress again that they are valid exactly in an infinite homogeneous elastic space and approximately in a half-space.Following Eshelby [1957] and Davis [1986], P a T /P, P b T /P, P c T /P can be obtained by solving the linear system (2) where x = − P a T /P, y = − P b T /P and z = − P c T /P, is the Carlson elliptic integral of the second kind (numerically evaluable, e.g.Press et al. [1992]).By the way, there is a typo in both the Equation (13) and the system of Equations ( 14) in the work by Davis [1986], where the right-hand side of the equation and the known term of the system of equations miss a division by 3.For a spheroid (a = b ≠ c), there are only two unknowns and the previous linear system becomes (since x = y) and as in Davis [1986] where Putting e = c/a, for a = b > c (e < 1, oblate spheroid) and a = b < c (e > 1, prolate spheroid), following Eshelby [1957] we have Taylor expansions of the expressions for We have seen that if the moment tensor obtained from seismic and/or deformation data represents a mixed mode dislocation on a fault plane or the expansion of an ellipsoidal cavity, ΔV is proportional to the ISO component of the moment tensor, but the proportionality factor is not a constant.Figure 1 shows the ratio (R iso ) of ΔV of the real source to the volume change (ΔV iso ) of an isotropic source sharing the same value of m.
Figure 1 also shows domains of possible ratios of eigenvalues for mixed mode (shear and tensile) dislocations on a planar fault (from pure compressive source, slope a = −90°, to pure shear source, slope a = 0°, to pure extensive source, slope a = 90°; Vavryčuk [2011]) and ellipsoidal pressurized cavities (light gray area), embedded in an elastic medium with o = 0.25.
In case of mixed mode dislocation, Again it is possible to compute both ΔV and the moment tensor in terms of ΔV: It is thus possible to compute both ΔV and the moment tensor in terms of ΔV: where

Coupled cavities
When dealing with expansion sources, volumecompensating processes involving mass movements should often be taken into account.For example, in volcanic environments sudden magma movement between coupled cavities can originate seismic waves.Amoruso and Crescentini [2009] have already shown that source properties retrieved from far-field deformation data for interconnected equally pressurized ellipsoidal cavities are less reliable than for single sources.The case of a source region including ellipsoidal cavi-ties interconnected by narrow fluid-filled conduits (whose contribution to the total volume change is negligible) can be treated under the assumption that the ellipsoidal cavities share the same pressure, and that the distance between cavity surfaces is generally larger than twice the cavity size (thus an ensemble of equally pressurized ellipsoidal cavities can be treated as if each cavity were isolated).Under the far-field approximation, moment tensors associated with the cavities inside the volume source can be considered co-located, and the overall volume change is given by It is interesting to note that R iso for ellipsoidal sources approaches R iso for tensile faults as the ellipsoid degenerates into a penny-shaped crack.Unfortunately, as previ-ously mentioned, if the moment tensor cannot represent an ellipsoidal pressurized cavity or a mixed mode dislocation, no univocal estimation of ΔV is possible.
as the ratio of ΔV of the real source to the volume change of an isotropic source sharing the same value of the product of pressure and volume.We get depends on the eigenvalue ratio and consequently on the ellipsoid shape.Amoruso and Crescentini [2009] obtained is constant and equal to 1.8 for o = 0.25, independently of the slope of the source.
In case of ellipsoidal pressurized cavities . (58) where is also the trace of the moment tensor that represents the far-field displacements from the ellipsoidal sources [Amoruso and Crescentini 2009]; is the bulk modulus.Here we just recall that, if the cavities have different shape and orientation, it is impossible to infer the volume change of an ensemble of ellipsoidal cavities from the equivalent moment tensor, if, as usual, PV cannot be estimated independently.
The total moment tensor representation, its decomposition into ISO, DC and CLVD force systems, and its trace, have been widely used not only to interpret surface deformations but also seismic records in volcanic areas.The ISO component is usually attributed to a spherical magma oscillation, the DC component to shear faulting, and the CLVD to rapid movement of magmatic fluid [e.g.Julian et al. 1998].
In this framework, Amoruso and Crescentini [2009] have considered the case of two coupled spheroidal cavities, and assumed that for some reason (e.g., fluid migration) volume of both cavities change and total volume change can be different from zero.They considered the case in which the symmetry axes of both cavities are supposed parallel to the z axis of the coordinate system and the case in which the symmetry axes of the two cavities are parallel to the z axis and y axis of the coordinate system respectively.The decomposition of the moment tensor into ISO, DC and CLVD force systems were performed using the method of Knopoff and Randall [1970], which makes the major axis of the CLVD coincide with the corresponding axis of the DC.Amoruso and Crescentini [2009] have shown that the moment tensor of earthquakes generated by a sudden magma exchange between two spheroidal cavities (e.g., because of the breaking of a barrier) can include all the three (ISO, DC, and CLVD) components.In particular, while the ISO component is given by the same expression (Equation 21 in Amoruso and Crescentini [2009]), the DC component is null if the symmetry axes of the two cavities are parallel to the z axis and not null in the other case.Although the interpretation of the moment tensor decomposition is by no means univocal, results in Amoruso and Crescentini [2009] can be used to interpret non-DC earthquakes.As an example, they re-examined three deep low-frequency earthquakes occurred beneath Iwate volcano, Japan, in 1998 and 1999 [Nakamichi et al. 2003], and showed that volume-change values obtained from the isotropic component only can be much smaller than those really involved in the magma exchange process.
Similar computations can be performed also for the moment tensor obtained from the inversion of surficial deformation data and different coupled sources, like spheroidal magma chambers and cracks, but their application to real cases is a case-to-case problem and a-priori assumptions about the involved sources are always necessary.
Here we focus on the simpler problem of how the volume change of each source can be obtained from the ISO component of the total moment tensor, in few cases of coupled sources where the total volume change is null.In other words, the volume change ΔV 1 of the first cavity (here always an ellipsoid) and the volume change ΔV 2 of the second cavity (a very thin oblate spheroid, a very elongated prolate spheroid, or a sphere) are such that ΔV 1 = −ΔV 2 = ΔV.If m 1 and m 2 are the ISO components of the two sources, the apparent total volume change (ΔV app ) given by the total ISO component is (59) independently of the orientations of the two cavities.In all the described cases, (60) We will see that ΔV app is not zero, its sign can be positive or negative, and it is proportional to ΔV.

Ellipsoid and very thin oblate spheroid
This case schematizes magma transfer between a magma chamber and a dike.The ISO component of the second source is (61) thus (62) Figure 2 shows ΔV app /ΔV for different shapes of the ellipsoidal chamber.In this case ΔV app /ΔV is always positive.

Ellipsoid and very thin prolate spheroid
This case schematizes magma transfer between a magma chamber and a closed conduit.The ISO component of the second source is (63) thus (64) Figure 3 shows ΔV app /ΔV for different shapes of the ellipsoidal chamber.In this case ΔV app /ΔV can be either positive or negative.

Ellipsoid and sphere
This case schematizes magma transfer between a magma chamber and a secondary chamber.The ISO component of the second source is (65) thus (66) Figure 4 shows ΔV app /ΔV for different shapes of the ellipsoidal chamber.In this case ΔV app /ΔV is always negative.

Surficial displacement generated by finite ellipsoidal cavities
In this section we give the expressions for surficial displacement (under the quadrupole approximation) generated by finite ellipsoidal cavities embedded in a homogeneous half-space [Amoruso and Crescentini 2011], after a brief recall of the necessary basic concepts.Amoruso and Crescentini [2011] mentioned that explicit expressions (in closed analytical form) for displacements and stresses can be given for homogeneous media but gave a more approximate solution through the combined effects of seven suitable point sources of appropriate location and strength, to make it usable in heterogeneous media.
The external displacement field due to a pressurized ellipsoidal cavity (axes a ≥ b ≥ c, parallel to the coordinate axes) embedded in an infinite elastic medium, is the same as that due to a uniform distribution of moment tensors located inside the cavity itself.In this special coordinate system, moment density is [Eshelby 1957] (67) For any other system it can be found by the usual law for tensor rotation.
The equation is a practically useful approximation also for a semi-infinite medium, provided that the ellipsoid is small with respect to its depth [Davis 1986].Thus, surficial displacement (v x , v y , v z ) of an observation point (x s , y s , 0), due to a pressurized ellipsoidal cavity centered in (x, y, z) and embedded in a semi-infinite elastic medium (z < 0), is given (approximately) by the volume integration over the ellipsoid (68) where is the Green's function giving i-component of displacement due to a very small (point) ellipsoidal source (unitary PV product) located in = = (x, y, z) on a point at = (x s , y s , 0), and the z-axis is perpendicular to the free plane surface.Similar equations also hold for stresses and strains, provided the proper Green's functions are used.Taylor series expansion up to order 2 of about the ellipsoid center is In a homogeneous half-space it is possible to express the displacement Green's functions, and thus their second derivatives , in analytical form.
Here we decompose the moment tensor repre-senting a point ellipsoid into those representing three orthogonal point cracks (tensile faults) whose openings are parallel to the ellipsoid axes.Other decompositions (e.g., in terms of double forces) can be used as well.Crack potency are Thus, since ∫ V (x j − x 0j ) dV = 0 because of symmetry, where ii) crack perpendicular to axis b i) crack perpendicular to axis a where Π a indicates potency of the crack perpendicular to axis a, and so on.We describe the spatial orientation of the ellipsoid by means of the three Euler angles a, b, and c, following the convention of Z-X-Z co-moving axes rotations (x convention in Goldstein et al. [2001]) for moving the (x, y, z) reference frame (here x Eastward, y Northward, z upward) to the (X, Y, Z) referred frame (X parallel to c, Y to b, Z to a).Because of the ellipsoid symmetry, Euler angle ranges can be restricted to −180° ≤ a < 180°, 0° ≤ b ≤ 90°, 0° ≤ c < 180°.
Direct computations show that strike (z, clockwise from North) and dip (d) of the three cracks are related to the Euler angles as follows: where i = x, y, z, Π is the tensile fault potency (hereinafter expressed as ∑Δu, surface times opening), and a i , b i are given in Equations ( 10) and (11) in Okada [1985].
As regards the second derivatives of the displacement, after some tedious algebra we get the following expressions, grouped according to the component of the displacement.
Using this decomposition, surface displacements and their second derivatives are computable from equations 10, 11, and 12 in Okada [1985].For the sake of clarity we use the Cartesian coordinate system by Okada [1985], i.e. right-handed, x-axis parallel to strike and z-axis upward.We put R 2 = (x 0 − x) 2 + (y 0 − y) 2 + z 2 and q = y sin d − z cos d, where d is the dip of the tensile fault.Final expressions can be can be easily implemented in computer codes and results can be transformed in any other coordinate system through proper axis rotation.The three components of the displacement due to a crack share the same structure iii) crack perpendicular to axis c (79) where where where and and and (1) (2) where where and and and where where where and and These closed analytical expressions can be easily implemented in any computation code.They are valid only in the Okada coordinate system, which is different from crack to crack, but can be used to obtain the correct results in any other coordinate system.The procedure is essentially the following: 1. given the Euler angles of the ellipsoid, determine dip and strike of the cracks that are perpendicular to the ellipsoid axes; 2. determine crack potency for unitary PV; 3. for each crack, compute the components of the displacement vector and the elements of the secondderivative third-order tensor in the Okada coordinate system using expressions in our work; 4. for each crack, compute the rotation matrix from the Okada coordinate system to the reference one; 5. for each crack, use the rotation matrix from step 4 to compute the components of the displacement vector and the elements of the second-derivative third-order tensor in the reference coordinate system.
If R is the proper rotation matrix from the Cartesian coordinate system in Okada [1985] to the reference one, the elements of the second-derivative third-order tensor can be transformed by using the expression (134) As previously mentioned, R depends on the Euler angles identifying the ellipsoid orientation and is different from crack to crack.

Conclusions
With respect to already published material, here we have presented the explicit expressions for the computation of the volume change and the generated surface deformation for finite pressurized ellipsoidal sources, in a closed analytical form to make their use easier in comparison to previous formulations.We have also evidenced and discussed a few results from the literature, that are often ignored, on the computation of the volume change for ellipsoidal sources and mixed mode dislocations on a fault; ignoring these findings might bias derived estimates (e.g.intrusion densities for volcanic environments).With the same aims, we have discussed how the volume change of each source can be obtained from the ISO component of the total moment tensor, in few cases of coupled sources where the total volume change is null.
for both cases where Taylor expansions of the expressions for I a = I b in case e > > 1 give

Figure 1
Figure 1 (left).Ratio (R iso ) of ΔV of the real source to the volume change (ΔV iso ) of an isotropic source sharing the same value of the moment-tensor trace; o = 0.25.Here m 1 and m 3 are the absolutely largest and smallest moment-tensor eigenvalues; m 2 is the absolutely intermediate one.The light gray area indicates the domain consistent with ellipsoidal sources.The piecewise solid line indicates the domain consistent with mixed-mode planar faults, from pure compressive source (slope a = −90°) to pure shear source (slope a = 0°) to pure extensive source (slope a = 90°).Note that |m 2 |/|m 1 | and |m 3 |/|m 1 | depend on |a| only.Figure 2 (right).Ratio of the apparent total volume change (ΔV app ) to the volume change (ΔV) of an ellipsoidal chamber in case of mass transfer to a dike, when the actual total volume change is null.Here a>b>c are the axes of the ellipsoid; o = 0:25.

Figure 3
Figure 3 (left).Ratio of the apparent total volume change (ΔV app ) to the volume change (ΔV) of an ellipsoidal chamber in case of mass transfer to a closed conduit, when the actual total volume change is null.Here a>b>c are the axes of the ellipsoid; o = 0.25.Figure 4 (right).Ratio of the apparent total volume change (ΔV app ) to the volume change (ΔV) of an ellipsoidal chamber in case of mass transfer to a spherical chamber, when the actual total volume change is null.Here a>b>c are the axes of the ellipsoid; o = 0.25.

∂
2 u i ∂x j ∂x k whatever = (R) il (R) jm (R) kn ∂ 2 u l ∂x m ∂x n Okada