Relations between pressurized triaxial cavities and moment tensor distributions

Pressurized cavities are commonly used to compute ground deformation in volcanic areas: the set of available solutions is limited and in some cases the moment tensors inferred from inversion of geodetic data cannot be associated with any of the available models. Two different source models (pure tensile source, TS and mixed tensile/shear source, MS) are studied using a boundary element approach for rectangular dislocations buried in a homogeneous elastic medium employing a new C/C++ code which provides a new implementation of the dc3d Okada fortran code. Pressurized triaxial cavities are obtained assigning the overpressure in the middle of each boundary element distributed over the cavity surface. The MS model shows a moment domain very similar to triaxial ellipsoidal cavities. The TS and MS models are also compared in terms of the total volume increment limiting the analysis to cubic sources: the observed discrepancy (~10%) is interpreted in terms of the different deformation of the source interior which provides significantly different internal contributions (~30%). Comparing the MS model with a Mogi source with the some volume, the overpressure of the latter must be ~37% greater than the former, in order to obtain the same surface deformation; however the outward expansion and the inner contraction separately differ by ~±10% and the total volume increments differ only by ~2%. Thus, the density estimations for the intrusion extracted from the MS model and the Mogi model are nearly identical.


Introduction
Ground deformation in volcanic areas is generally recognized as a reliable indicator of unrest, possibly resulting from the intrusion of fresh magma within the shallow rock layers.Following this suggestion, deformation sources are generally modeled as pressurized cavities endowed with some prescribed geometrical shape.The most popular of these models is the Mogi [1958] source which describes the deformation due to a spherical cavity with radius much smaller than its depth; McTigue [1987] provided corrections accounting for the finite radius of the spherical cavity.Bonafede and Ferrari [2009] have shown some hidden implications of the well known equivalence between isotropic moment sources, Mogi sources and a set of three orthogonal dislocations: very different estimates of the incremental volume are predicted, because the nearfield displacement is involved.
Solutions for the deformations associated with more general triaxial ellipsoidal cavities were proposed by Davis [1986], under the point-source assumption.Yang et al. [1988] and Cervelli [2013] provided approximate solutions for finite ellipsoids in a homogeneous half-space.Fialko et al. [2001] provided semi-analytic solutions for a finite horizontal circular crack.Bonafede and Ferrari [2009] generalized the Mogi source to a viscoelastic half-space and to a spherical viscoelastic shell embedding the source in an otherwise elastic halfspace, providing an interpretation of a pressurized cavity in terms of dislocation sources and related moment tensor.These solutions are used to infer from geodetic data the depth and the (incremental) volume of the intrusion, which are very important parameters for volcanic risk implications.
However, it is well known [e.g., Backus andMulcahy 1976a,b, Aki andRichards 1980] that any internal source of deformation can be described in terms of a suitable moment tensor, but there is no general scheme of inversion from which the geometrical shape of an "equivalent pressurized cavity" can be inferred: for example the moment tensor domain pertinent to pressurized ellipsoids is very limited [e.g., Trasatti et al. 2011].Furthermore important parameters such as the overpressure within the cavity and the source volume expansion are related to the near-field deformation which can be very different for different source shapes.Overpressure estimates are important to assess the possibility of dyke opening and fracture propagation in the region surrounding a magma reservoir; source volume changes are important to infer density estimates of an intrusion from gravity measurements.In the present paper, we use parallelepiped cavities to extend the correspondence between triaxial pressurized cavities and moment tensors: we show that triaxial parallelepiped cavities are mostly equivalent to ellipsoidal cavities [Davis 1986], when surface deformation is computed in the far-field.Two different models are presented (pure tensile sources and mixed mode sources) which are studied employing a boundary element approach [Crouch and Starfield 1983] and the solutions given by Okada [1992] for the internal displacements and strains due to rectangular dislocations buried in a homogeneous halfspace.Furthermore the boundary element method allows simulating triaxial pressurized cavities with linear dimensions comparable or greater than depth.

Source models
In this section, we introduce two models of deformation sources which are studied using a boundary element approach.The boundary element method is implemented in a C++ code and two different implementations of Okada expressions are used to compute the stress and the displacement induced by finite rectangular dislocations: -the original fortran code (dc3d.f ) written by Okada [1992]; -a new C/C++ code optimized for boundary element methods.
The results obtained employing the new implementation are compared with those obtained with the original code to check their correctness and to test the computation speed improvements.
The parameter N determines the total number N TOT of dislocation elements that are used to discretize the six faces of the parallelepiped source (l 1 , l 2 , l 3 length edges and X i ± unit normal vector in the direction of the x i axis, Figure 1).An algorithm was designed to achieve a uniform distribution of nearly square boundary elements on the cavity boundaries.
Two types of sources are considered: -Tensile source (TS): on each face of the cavity only tensile dislocations operate, with Burgers vector normal to the dislocation surface (Figure 2).The displacement discontinuity components (b T ) associated to each element are found solving the linear system of N TOT equations in N TOT unknowns: (1) where T stands for tensile dislocation and: -(v i T ) 0 are the prescribed boundary tractions, normal to dislocation area; -A ij TT are the influence coefficients obtained from Okada [1992].
In this model, the problem is defined assigning the normal component of traction v T at the center of each  boundary element and then it is solved for the b T component of the Burgers vectors.
-Mixed source (MS): on each face a mixed (tensile, strike-slip and dip-slip) dislocation operates (Figure 2).In this case, the displacement discontinuity components (b T ,b S ,b D ) associated to each element are found solving the linear system of 3•N TOT equations in 3•N TOT unknowns: where T,S,D stand for tensile, strike-slip and dip-slip dislocations, respectively, and: - Okada [1992].
In the MS model, the problem is defined assigning the three components of traction v T , v S , v D at the center of each element and then it is solved for the b T ,b S ,b D components of the Burgers vectors.If a pressurized cavity is assumed, then v i S = v i D = 0 but non vanishing values apply if fluid intrusion takes place in a prestressed medium [see, e.g., Trasatti et al. 2011].In this model strike-slip and dip-slip dislocations operate and so the source is allowed to shift as shown in Figure 3.
Several tests have been performed which show a general consistency between the results obtained using our implementation of Okada expressions and the original code by Okada (dc3d.f ).A major issue is found for the MS model when the depth of the source is comparable with its dimensions: the cavity interior is subjected to a major rigid body translation along the vertical di-rection which implies non-physical results (e.g., matter interpenetration).The displacement in the external domain is univocally determined imposing null displacement at infinity, but the solution in the internal domain is undetermined by a rigid body translation in accordance with the uniqueness Kirchoff theorem.A detailed analysis (see Appendix A for details) has shown that the numerical problem can be solved fixing the displacement of the interior region in order to prevent it from translating or rotating arbitrarily.
Finally, the tests show that significant speed improvements can be achieved using our C++ implementation of Okada expressions (see Appendix B for details) and so in the following we shall refer only to the new code to perform computations.

Point source approximation
Before we proceed, it is necessary to describe how a single moment tensor can be associated with each TS and MS source.If N=1, over each face a rectangular dislocation is considered with surface area A ± i (where ± denotes the two faces normal to x i , see Figure 1 According to Kirchoff uniqueness theorem [e.g., Fung 1965], outside a pressurized parallelepiped the deformation field is identical to that provided by 6 pressurized rectangular cracks over its faces: -according to the TS model, these cracks may be approximated as 6 tensile dislocations and their Burgers vectors b ± i are computed in order that they provide the same overpressure at the mid point of each face (see linear system (1)).It must be stressed that the TS model is not exactly pressurized since non-vanishing shear stress may appear over the cavity boundary; -according to the MS model, these cracks may be represented by 6 mixed dislocations and their Burgers vectors b ± i are computed in order that they provide the same overpressure (v i T ) 0 and vanishing shear tractions (v i S ) 0 , (v i D ) 0 at the mid point of each face (see linear sys-  tem (2)).Accordingly the MS model should provide a better approximation of a pressurized cavity.
In the case of the TS model, the six tensile dislocations are equivalent, in the point-source approximation, to 3 orthogonal tensile dislocations, located in the center of the cavity, with surface areas and Burgers vectors The moment tensor (describing these three orthogonal tensile dislocations) is simply obtained (employing the axes as basis vectors) from the theorem of body force equivalents [Burridge and Knopoff 1964] The point source (PS) so defined is located in the center of the parallelepiped cavity and it is represented by a single moment tensor M ij , Equation (4).
Instead, in the case of the MS model, the six mixed dislocations are equivalent, in the point-source approximation, to 3 orthogonal mixed dislocations, located in the center of the cavity, with surface areas and Burgers vectors The moment tensor (describing these three orthogonal mixed dislocations) is computed as follows according to the theorem of body force equivalents; again, the point source (PS) which we have defined is located in the center of the parallelepiped cavity and it is represented by a single moment tensor M ij , Equation (6).
When N>1, Equation ( 4) is still valid in the far field (point source approximation) for the TS model if we put 4  1 , X ± 2 , X ± 3 , respectively.In the same way, when N>1, Equation (6) remains valid for MS sources in the far field if we put where the mean value of Burgers vectors components are computed with formulas similar to Equations (7).
In this section, we have realized an association between a single moment tensor (PS) and a finite source (FS) which represents instead a moment tensor distribution; this correspondence is certainly acceptable when we want to evaluate the displacement and the stress in the far field.Instead, in the near field, this approximation must be checked in order to evaluate if and where it is acceptable.
For this aim, we consider three different source configurations (Figure 4) and we plot to evaluate the error introduced using a point source (PS) to evaluate the deformation (u z vertical component, u r radial component) at the free surface (z=0) produced by the finite source (FS).We expect that the difference between PS and FS is significant only when the source dimension is comparable to depth.Let h denote the ratio between depth d and the largest edge, h=d/max(l 1 ,l 2 ,l 3 ).First, we consider the case N=1; in this case the models TS and MS give practically the same results.For h=5, the point source approximation is still acceptable (Figure 5); in each source configuration, the absolute value of du z and du r never exceeds 5%.The major discrepancy is observed for |x 1 |/d<0.5,|x 2 |/d<0.5 while, outside this region, du z and du r are negligible.More precisely, above the source at the surface, if we use the point source approximation, u z , u r are underestimated in the case of the cubic source (a), overestimated in the case of the sill (b).Instead, in the case of the dike (c), no significant discrepancy is observed.The sign and the magnitude of the error are dependent on the shape of the source and more precisely on the dimensions and the position of X ± 3 faces.In the case of the cubic source, the error is negative (underestimate) since the point source is located in C(0,0,−d) and so the effect produced on ground deformation by the horizontal face X ± 3 (located at depth z=−d+l 3 /2) is underestimated since the asymmetric opening of the tensile components over X ± 3 faces is not taken into account.In the case of the horizontal sill, the effects produced by X ± 1 and X ± 2 faces are negligible and so under the point source approximation the concentration of the total moment in C(0,0,−d) determines a little overestimate of ground deformation over a small central region; however the error introduced is negligible since the source shape (l 3 << l 1 = l 2 ) and its horizontal orientation makes the source center depth truly representative of the source depth.In the case of the vertical dike, the effects produced by the X ± 3 faces are very little since their area is negligible and so in this case the point source approximation is quite accurate (Figure 5c).
If h is smaller, N>1 must be considered and the models TS and MS are no longer equivalent (Figure 6).The MS models show a better agreement with the point source solution since a significant discrepancy is observed only in the region |x 1 |/d<0.5,|x 2 |/d<0.5 (|du z |>5% and |du r |>5%) and |du z |,|du r | never exceeds 15%.In order to obtain an accurate evaluation of the surface displacement, N may be reasonably low (N=15 in Figure 6); however we shall see that higher values of N must be used if the near-field displacement is needed, such as in the evaluation of source boundary expansion.

Moment tensor domain and pressurized sources
As previously observed, triaxial ellipsoids [Davis 1986] and circular horizontal cracks [Fialko et al. 2001] are generally employed in the inversions of ground deformation and gravity data to infer the geometrical shape of deformation sources in volcanic areas.However these models do not represent the most general internal sources, since they can generate only a limited subset of moment tensors [Trasatti et. al 2011].
In the following we investigate if this subset can be wider when pressurized parallelepiped cavities are considered.We consider the case of a source embedded in a homogeneous space and we fix the side length l 1 , l 2 , l 3 to determine how the eigenvalue ratios M 2 /M 1 , M 3 /M 1 are dependent on the edge ratios l 1 /l 3 , l 2 /l 3 of the cavity.The choice of the vertical edge l 3 as reference is arbitrary.
For both TS and MS models, when N=1, the domain in the plane M 3 /M 1 vs M 2 /M 1 is a triangular area (Figure 7a) which extends considerably (Figure 7c) the domain covered by triaxial ellipsoid cavities (Figure 7b).The moment ratios M 2 /M 1 and M 3 /M 1 must be positive and ratios lower than 1/3 cannot be obtained (in the Poisson approximation).This map has been retrieved PRESSURIZED CAVITIES AND MOMENT SOURCES  In panel (a), we show the non uniform distribution in the moment ratios domain (M 3 /M 1 , M 2 /M 1 ) of the different configurations (correspondent to different edge ratios) of a rectangular parallelepiped cavity.On each line of the graph the value of l 1 is fixed and moving from right to left the value of l 2 is increased by a constant step from the value l 1 to the value fixed for l 3 .In panel (b), the map identifies the distribution of the different source configurations of a rectangular parallelepiped source in terms of the value of the ratio l 1 /l 2 .The top vertex of coordinates (1,1) (blue dot) corresponds to the cubic cavity which yields the same result as a spherical cavity (in the point source approximation).The vertical side of the triangular domain corresponds to the parallelepiped rectangular cavities having l 1 =l 2 <l 3 (prolate sources) and the vertex (1,0.5)(pink dot) is associated to cavities which have l 1 =l 2 <<l 3 .The oblique side, which connects the top vertex with the bottom vertex of coordinates (0.33,0.33) (yellow dot), represents cavities with l 1 <l 2 =l 3 (oblate sources) and indeed for l 1 <<l 2 =l 3 the rectangular parallelepiped cavity can be identified with a tensile square crack.considering a set of parallelepiped cavities characterized by different l 1 /l 3 , l 2 /l 3 ratios.These points are not uniformly distributed over the triangular domain (Figure 8a), with density increasing from prolate sources to oblate sources.This relation between density and source geometry implies that, for an oblate source, minor changes in l 1 /l 3 , l 2 /l 3 ratios determine significant variations of M 3 /M 1 , M 3 /M 1 ratios.If we consider Table 1, we see that M 3 < M 3 < M 1 if the parallelepiped edges are in the reverse order (l 1 < l 2 < l 3 ), which is the same qualitative behavior generally observed for a triaxial ellipsoidal cavity with respect to the order of its axes.This is the reason why l 3 is taken as the reference edge and M 1 as the reference moment.Now the TS and the MS models are studied in terms of the discretization parameter N. First, we compare the results obtained when N=1 and N=15 if we consider a TS parallelepiped cavity (Figure 9a).The domain pertinent to rectangular parallelepiped cavities is still the triangular region (Figure 10a) which extends considerably (Figure 10d) the domain covered by triaxial ellipsoidal cavities (Figure 10c).But for every source configuration (namely for fixed ratios l 1 /l 3 , l 2 /l 3 ) the moment ratios are dependent on the value of N (arrows in Figure 9a).The different configurations are affected in different ways by N since the direction of the arrows does not show a uniform pattern.The axially symmetric configurations (l 1 = l 2 , vertical line and l 3 = l 2 , top oblique line) are quite insensitive to the increase of N.  In the case of the MS model, N has a greater impact with respect to the TS model (Figure 9b).The shift increases progressively as we move from the top vertex (1,1) to the bottom oblique side.Furthermore the arrows point towards the top vertex (1,1) for all the source configurations, so a significant reduction of the domain takes place when N>1 (Figure 10b).MS parallelepiped cavities and ellipsoids have very similar domains in the space of moment ratios when N>>1 (Figure 10e).As shown by Amoruso and Crescentini [2011], the source mechanism can be described by a multipole expansion: the monopole term characterizes the far field and the dipole terms vanish if the source exhibits mirror symmetries.Then, MS parallelepipeds and ellipsoidal cavities differ only for terms higher than the dipole.

Isotropic source: volume increment estimate
The inflation of a magma chamber may be produced by internal differentiation processes or by the intrusion of fresh magma; the volume and the density of the intruded material are very important parameters to discriminate between the two processes and to evaluate the volcanic risk.
In our models represents the total volume increment of the parallelepiped chamber due to the opening of the tensile dislocations distributed over its surface.The total volume increment can be split into the sum of two contributions, DV in and DV out , where DV in represents the increment of volume produced by the contraction of the source interior and DV out is the outward expansion over the source boundary.
In the case of a cubic cavity, we compare the estimation of DV m obtained with the TS and the MS models for increasing values of the discretization parameter N (Figure 11a).
For low values of N, the two models are in very PRESSURIZED CAVITIES AND MOMENT SOURCES good agreement but, for N>3, DV m decreases rapidly, particularly for the TS model, until N~15 and continues to decrease slowly for N>15.The discrepancy between the TS and the MS models increases rapidly up to N~15 and slowly for N>15 (Figure 11b).When N=40 the solutions are not fully stable but it is difficult to study for N>40 since the linear system to be solved has dimension ∝N 2 .The overestimation of DV m (due to the overestimates of Burgers vectors) decreases with the increase of N (Figure 11a) as it is usual employing boundary element methods.
In order to understand the discrepancy between DV m estimates, we now compare the separate estimates of DV out and DV in performed with both models (Figure 12).
A very good agreement between TS and MS is observed for DV out .So the difference between DV m estimates must be attributed mostly to the different deformation of the source interior (Figure 11).Indeed, DV in estimates relative to TS and MS models are significantly different (Figure 12b).
To study the behavior of the source interior, we compare the results proper to the TS and the MS models for two different values of the parameter N=1,15 (Figure 13).value pertinent to a pressurized body of arbitrary shape with bulk modulus K).For N=1, the results relative to TS and the MS models are coincident and this fact is consistent with the results of Figure 5 and with the results of the previous section (coincident moment domain for TS and MS when N=1).When N=1, the internal deformation is far from a state of isotropic compression and this result is not surprising since the overpressure v T i is assigned only in the center of each face.To obtain a state of uniform compression, the boundary conditions must be applied uniformly over the cube boundary and, to approach this configuration, N must be increased.For N=15, on each face the boundary conditions are applied in 225 points and quite uniform values for de are obtained; therefore a state of isotropic compression is quite correctly reproduced for both models.The MS model overestimates the isotropic value e I but the relative difference is everywhere less than +10%, apart from the edges of the cavity which are singular domains.In-stead the TS model predicts a significant minor contraction of the interior since e kk is everywhere less by 20% than the theoretical value e I .
In the case of a pressurized sphere with reference volume V 0 , we know that [Bonafede and Ferrari 2009] and so the total volume increment is equal to If we put n=m (Poisson approximation), we obtain DV in /DV m =4/9.Now we study DV in /DV m in terms of the parame-

PRESSURIZED CAVITIES AND MOMENT SOURCES
Figure 13.To study the internal deformation of the cubic cavity, we plot the variable de which measures the relative difference between the computed value and the theoretical value proper to a uniformly pressurized cavity of arbitrary shape.We compare the results obtained using the TS (first row) and the MS (second row) models for two different values N=1 (left column) and N=15 (rigth column) of the parameter N.
ter N (Figure 14): this ratio rapidly converges to 0.29 (for the TS model) and 0.36 (for the MS model) and this difference is easily explained since DV TS out DV MS out (Figure 12a) and DV TS in <DV MS in (Figure 12b).For both models, the results indicate ratios significantly lower with respect to the case of a pressurized sphere.Thus, for the Mogi source, the two contributions to the total volume increment DV m are similar in size (DV in /DV out =4/5), while the DV out contribution is significantly larger than DV in in the case of TS (DV TS in /DV TS out =0.29/0.71=0.41)and MS (DV MS in /DV MS out =0.36/0.64=0.56)models.Therefore the shape of the source controls the near field displacement implying a different partition of the two contributions to the total volume increment.The different values for the different models are not due to discretization problems (low N), but to the different source geometries (Mogi vs MS) and to the different boundary conditions (TS vs MS).
The analysis performed for a cubic cavity in terms of the volume increment can be extended to the case of a generic triaxial parallelepiped cavity.In this case, the results relative to TS and MS models must be compared with the exact (for a homogeneous elastic medium) ellipsoidal cavity results given by Amoruso and Crescentini [2009].

Residual gravity changes: Mogi source vs cubic source
To show an application of the MS model, we compare the residual gravity computed according to a Mogi source and a cubic source.Very similar results would be obtained employing the TS model.Preliminarily, we fix depth (d=5 km) and side length (h=1 km) of the cubic source and we assign v T i = 185 MPa, v S i = v D i =0 (i=1, ... ,6N) as boundary conditions (pressurized cavity) in order to reproduce the maximum uplift (~1.8 m) observed in the case of the 1982-84 uplift episode at Campi Flegrei (n=m=1 GPa as elastic parameters).Then, we consider a Mogi source with center at the same depth and we adjust its strength (DP•V 0 ) in order to reproduce the maximum uplift determined by the cubic source at the free surface.In this way, we find that the cubic source and the spherical source show very similar DV m expectations, DV MS m =333.9•10 6m 3 and DV Mogi m =342.6•10 6m 3 (relative difference between estimations is 2%).Thus, even if the MS source provides a source expansion +13% greater than the Mogi source, the different DV in compensates mostly this difference.Instead the product (DP•V 0 ) is significantly different between the Mogi source and the MS source.Indeed, when we fix the same source volume, the two models produce the same maximum uplift if the overpressure relative to the Mogi source is increased to the value DP=253.8MPa (DP Mogi /DP MS = 1.37).Equivalently, in the opposite case, when we fix the same overpressure (DP=185 MPa), the volume of the Mogi source must be increased to the value V 0 = 1.37•10 9 m 3 with respect to the volume l 3 = 1•10 9 m 3 of the cubic source.It must be stressed that the ratio (DP Mogi •V 0 Mogi )/(DP MS •V 0 MS ) becomes stable for N>5.Now we compare the two sources in terms of the residual gravity computed at the point of maximum uplift (w=u z (0,0,0)); after we correct for the free air gradient, the residual gravity Dg R can be split in the sum of three contributions: where -Dg S represents the gravity change induced by the intrusion from remote distance of new material with density t m into the cavity.This term is evaluated in the following way -Dg V is the gravity change determined by the deformation of the material (density t) surrounding the source volume The integration domain V is discretized in a grid of N GRID volume elements and the integral is computed The surface integral ( 15) is computed as the sum of the contributions over each boundary element.The integrals ( 13) and ( 15) are computed according to the numerical procedure described by Trasatti and Bonafede [2008].
For a Mogi source in a half-space, Walsh and Rice [1979] proved that Dg V +Dg B vanishes identically and our numerical results employing a spherical source are consistent with this result and the discrepancy between the absolute values of Dg V and Dg B is mostly due to the finite domain (100 km (X 1 ) × 100 km (X 2 ) × 50 km (X 3 )) used for the numerical computation of the integral (13).
For the MS source, we obtain the following results which show a little but significant discrepancy (Dg V + Dg B )/w = 2nGal/m that cannot be ascribed to the finite domain used for the numerical integration.This discrepancy is due to the non-isotropic effects of the free surface, since the top and the bottom faces suffer different inflation; the cubic source and the spherical source differ for quadrupole terms (and higher) as previously observed in the comparison of parallelepiped sources and ellipsoidal cavities.
At Campi Flegrei, in the case of the 1982-84 uplift episode, the residual gravity change measured at the point of maximum uplift was Dg R (obs)/w(obs)=(74±12) nGal/m [Berrino 1994]; using this value, we obtain from Equation ( 11), that t Mogi m =1471kg/m 3 ≃t MS m =1469 kg/m 3 since the difference between DV Mogi m and DV MS m is largely compensated by a non vanishing value of Dg V +Dg B in the case of the MS source.It must be stressed that in many papers the Mogi source is used to infer the volume of the intruded material neglecting the isotropic contraction of the source interior.Since DV in and DV out are comparable in size, the severe underestimation of the total volume increment determines a strong bias in density estimation.For example, in our case, if we consider DV Mogi out as total volume increment, we obtain t Mogi m =2649 kg/m 3 which strongly overestimates the previous value.Instead, in the case of the MS model, the total volume increment DV m is directly computed using the Burgers components of the boundary elements and so the density estimation is valid since the isotropic contraction of the interior is automatically taken into account.

Conclusions
Two models were presented to generalize to triaxial geometries the equivalence between pressurized cavities and moment tensor sources.The models refer to the same simple source shape (rectangular parallelepiped) but differ for the boundary conditions assigned on their surface: the TS model assumes tensile dislocations while the MS model allows for mixed (tensile/shear) dislocations providing an assigned overpressure and vanishing shear tractions in the middle of each boundary element.
The correspondence between finite source models and point sources shows that the TS and the MS models are equivalent only if N=1, where N represents the discretization parameter in the boundary element approach, provided that an internal element is added which constrains the translation of the MS interior (see Appendix A).Instead, if N>1, these models are no longer equivalent, particularly if we consider sources close to the free surface for which the point source approximation fails.
The study of the models in the space of moment tensor eigenvalues confirmed the equivalence if N=1.The domain pertinent to pressurized parallelepiped cavities in a M 3 /M 1 vs M 2 /M 1 plot is a triangular area which extends significantly the domain proper to point like triaxial ellipsoidal cavities (Figure 7).If N>1, the MS domain shrinks towards the domain pertinent to pressurized ellipsoidal cavities (see Figure 10).The domain proper to the TS model is unchanged with respect to N=1: this is ascribed to the generation of shear PRESSURIZED CAVITIES AND MOMENT SOURCES ( stresses over the cavity boundary, a non-physical implication for a fluid filled cavity.The total volume increment DV m available to host new magma within the source is the sum of the outward expansion DV out of the source and the inward contraction DV in of the material already resident within the source.The DV m values predicted by the TS model and the MS model tend to differ substantially if the discretization parameter N increases.Limiting the analysis to a cubic source, the difference between the two DV m estimates (~10%, Figure 11a) must be attributed almost completely to the different internal deformation (~30%, Figure 13b) of the two sources.This different behavior is a direct consequence of the different boundary conditions imposed.Instead the DV out values predicted by the tensile source and the mixed source models are nearly coincident, even if both decrease when the parameter N increases (Figure 13a).
The TS and MS models are compared with the Mogi source in terms of the ratio DV in /DV m (Figure 14).The different source shape implies different near field displacements which in turn determines a different partition of the total volume increment DV m between the external contribution DV out and the internal contribution DV in .The cubic shape determines a larger DV out contribution and a lower DV in contribution with respect to a Mogi source with the same moment.
The application of the Mogi model and the MS model to the case of the 1982-84 uplift episode at Campi Flegrei shows different strengths (DP•V 0 ) for the two sources.The lower value expected for the cubic source with respect to the spherical source implies a minor overpressure over the cubic source if the same reference source volume V 0 is considered.Instead, if the same overpressure is assigned, the Mogi source requires an initial volume +37% greater than the MS model in order to obtain the same surface deformation.The ratio between the source strengths DP Mogi •V 0 Mogi / DP MS •V 0 MS is not related to the value of the discretization parameter N if its value is greater than ~5.
The computation of the residual gravity allows to estimate the density of the intruded material.The density estimated at Campi Flegrei using the Mogi source and the MS model are in good agreement if the contraction of the source interior is properly taken into account.However density estimates are strongly dependent on the source shape if triaxial geometries are considered as shown by Currenti [2014].The present method can be easily employed to generalize the partition of DV m between the internal and the external contributions (DV in and DV out ) and the overpressure estimates DP for triaxial cavities.

A. Discussion about numerical problems
To study the numerical issues related to the TS and MS models, it is convenient to consider the cubic cavity, for which we expect the simplest output, being the most symmetrical configuration.Furthermore, we fix N=1 in order to reduce the number of parameters to -6 components for the TS model, -18 components for the MS model.With X i + , X i − we indicate the faces with unit normal vector in the direction of the x i axis (Figure 1).To perform computations, we fix l 1 =l 2 =l 3 =100 m and we prescribe an overpressure v i T =10 MPa in the middle of each face of the cube described as a dislocation surface.Fixing the depth of the source to the value d = 0.5 km, the computations performed with our boundary element code give the output summarized in Table 2.
Considering the TS model, the results obtained with the two implementations are fully consistent; the asymmetry in the tensile component over the horizontal faces is expected since X 3 + , X 3 − are at different depths.In the case of the MS model, the tensile components of X 1 ± , X 2 ± faces are quite correctly identified.The pertur-bations with respect to the value 0.117925 m (the isotropic value valid for h>>1) are related to a small rigid body translation of the cavity interior along the horizontal directions.Indeed, the S components on X 1 ± , X 2 ± and both the S and D components on X 3 ± (~10 -7 − 10 -8 m, F&B code and ~10 -6 − 10 -7 m, Okada code) are small but significative.Instead the D components on X 1 ± , X 2 ± faces (±1340.19m, F&B code and ±3423.32m, Okada code) show that the cavity interior is subjected to a major rigid body translation along the vertical direction.This solution is not physically acceptable since the cavity interior would cross the cavity boundaries and even the free surface providing an altered surface deformation pattern.The motivations under this non-physical behavior must be understood and a mechanism to prevent it must be introduced.As pointed out by Crouch and Starfield [1983] for 2D problems, when the boundary element approach is applied over a closed contour, the problems involving the interior and the exterior regions are solved together and therefore it is necessary to fix the displacement of the interior region in order to prevent it from translating or rotating.In our 3D problem, we put one additional boundary element inside the cubic source and so the displacements in the interior region are determined relative to this fixed point while, in the exterior region, they are determined relative to null displacements at infinity.Once the problem is solved, it is possible to show that the additional internal element is subjected to certain shear and normal stresses which however do not affect the solution of the external region.The dimensions of the additional boundary element must be comparable (less or equal) to the dimensions of the other boundary elements.The choice of the position and the orientation of the internal boundary element is made considering the symmetry properties of the source.In order to preserve the symmetry of the output TSD components, the first attempt was to put the additional element at the center of the source and parallel to the free surface.However the implementations of Okada expressions obey to the practical rules described in Okada [1992]: in order to have non vanishing displacements for the internal boundary element it is necessary to shift it from the center of the source.After several tests, we found that the convenient choice is to move the additional element only a few centimeters away from the central point of the internal region.
In Table 3, we compare the results pertinent to the MS model when the additional boundary element is added in the interior region.The insertion of a boundary element inside the cavity solves the issue related to the translation along the x 3 axis since the vertical translation is now correctly identified.The Burgers vectors are consistent using the two different implementations of Okada expressions and furthermore, removing the rigid body translation, the tensile components of X 3 + , X 3 − faces X 3 + 0.119533 − 0.001211 = 0.118322 m X 3 + 0.116953 + 0.001211 = 0.118164 m reproduce the value identified according to the TS model (see Table 2).
As a final remark we have to establish why the results relative to the vertical translation of the cavity interior are so different (1340.19m, F&B code and 3423.32 m, Okada code) when the position of the MS source is not fixed (Table 2).The answer to this question is immediate if we consider the reciprocal condition number (RCOND, see Appendix C) of the linear system relative to the MS model.RCOND is ~10 -3 and so this linear system is very ill conditioned.Therefore the discrepancy must be attributed only to the ill conditioned nature of the problem and the two implementations of Okada expressions must be considered consistent with each other.

B. Details about programming routines
In this appendix we summarise the main features of our implementation of Okada expressions (F&B C++ code) which differs in some aspects from the original Okada implementation (dc3d fortran code).Our implementation is presented with two different interfaces which are designed to satisfy different needs.
The first interface reproduces the original interface of the dc3d Okada code:  The only major difference between our code and the one from Okada is related to the implementation of the rules described in the conclusions of Okada [1992] article (but not introduced in his code) which allows in some cases to refer to different subsets of the formulas in order to compute displacement and stress fields: -the internal deformation field (z<0) is computed according to Okada as the sum of four contributions: • two unbounded medium terms f A (p,h,z), f A (p,h, −z) • two regular depth dependent terms f B (p, h, z), f C (p, h, −z) which must be combined in the following way The symbol || represents the notation introduced by Chinnery [1961] where p is the coordinate along the dip direction.
-the surface deformation (z=0) is computed according to the subset (29) This first interface allows programmers to create new programs using a well known interface and gives the opportunity to use our implementation by means of few changes in old programs that uses the dc3d code.
However this interface is not the best choice when boundary element methods (such as the displacement discontinuity method) are implemented.Indeed, when the influence coefficients are calculated, the simultaneous computation of displacement and strain components is unnecessary and simply determines a waste of computer resources (see Table 4).According to this consideration, we make available a second interface which is designed for codes using boundary element methods and it consists of two separate functions that compute displacement (dc3d_FB_displ) and strain (dc3d_FB_strain) components separately: Table 4.In the table we compare the execution time requested to compute: (a) the surface deformation pattern (z=0), (b) the influence coefficient matrices relative to a TS cavity and a MS cavity using our implementation (F&B C++ code) and the original Okada Fortran code.The discretization parameter is N=15 and the source considered is a cubic cavity.Thus the number of the boundary elements is 6 N 2 =1350 and the matrices have the following dimensions: TS model, 1350  The IRET variable gives as output the return code of the functions (=0....NORMAL, =1....SINGULAR).
Both functions return separately the displacement and the strain components which are devoted to the strike-slip (S), the dip-slip (D) and the tensile (T) components.This choice reduces execution time when we have to compute more than one of the S, D, T components.The use of the dc3d Okada code requires instead the separate computation of these influence coefficients which determines a longer execution time (see Table 4).
Some minor speed improvements are also obtained using some small modifications introduced by Becker et al. [2005] in a slightly modified version of Okada fortran routines.
For fortran programmers, further informations about how to invoke our C/C++ interfaces are given directly in the comments of our source code which we provide as supplementary material of the paper.
A version of our interfaces is also available in the GNU Octave language [Eaton 2002] and in the future a Matlab version will also be released.

C. Details about the computation of the reciprocal conditional number
To compute the reciprocal conditional number, we adopt the same approach used in the GNU Octave language [Eaton 2002] by the function rcond().This function computes the 1-norm estimate of the reciprocal conditional number as returned by Lapack [Anderson et al. 1999].
The steps performed by the procedure are the following: 1. the 1-norm of the matrix is computed using the function dlange(), 2. the LU decomposition of the matrix is obtained using the function dgetrf(), 3. the reciprocal norm is found invoking the function dgecon().

Figure 1 .
Figure 1.Orientation of the cracks which cover the faces of the cubic source.The source is observed from two different points of view in order to increase the readability of the illustration.

Figure 2 .
Figure 2. Tensile (T), strike-slip (S) and dip-slip (D) components of a rectangular inclined fault (d, dip angle).In our models the dip angle d assumes only two values: 0°, in the case of the horizontal faces and 90°, in the case of the vertical faces of the parallelepiped cavity.

Figure 3 .
Figure 3. Translation of the source interior determined by the slip of strike-slip and dip-slip dislocations located over the faces of the cubic parallelepiped cavity.

Figure 4 .
Figure 4. We consider three source configurations: (a) a cubic source, (b) a square horizontal sill, (c) a square vertical dike.In each case, the center of the source is located in C(0,0,−d) (d>0).
mean values relative to the three distributions of boundary elements correspondent to the couple of faces X ±

Figure 5 .
Figure 5.We fix N=1 and we consider du z (left column) and du r (right column) relative to three different source configurations: a cubic source (a), a horizontal square sill (b), a vertical square dike (c); the source depth d is fixed in order to have h=d/max(l 1 ,l 2 ,l 3 )=5.In the figure, all the distances are normalized with respect to d.

Figure 6 .
Figure 6.We fix N=15 and we compare du z (first row) and du r (second row) computed according to the TS model (left column) and according to the MS model (right column) in the case of a cubic source with a depth d comparable with its dimensions (h=d/l=2).

Figure 7 .
Figure 7.The domain in the space of moment ratios (M 3 /M 1 , M 2 /M 1 ) proper to parallelepiped rectangular cavities (panel (a), coincident for TS and MS models if N=1) is compared with the domain proper to ellipsoid sources (panel (b)).In panel (c), the area with single line pattern represents the extension provided by rectangular parallelepiped cavities to the ellipsoids domain (crossed lines pattern).

Figure
Figure8.In panel (a), we show the non uniform distribution in the moment ratios domain (M 3 /M 1 , M 2 /M 1 ) of the different configurations (correspondent to different edge ratios) of a rectangular parallelepiped cavity.On each line of the graph the value of l 1 is fixed and moving from right to left the value of l 2 is increased by a constant step from the value l 1 to the value fixed for l 3 .In panel (b), the map identifies the distribution of the different source configurations of a rectangular parallelepiped source in terms of the value of the ratio l 1 /l 2 .The top vertex of coordinates (1,1) (blue dot) corresponds to the cubic cavity which yields the same result as a spherical cavity (in the point source approximation).The vertical side of the triangular domain corresponds to the parallelepiped rectangular cavities having l 1 =l 2 <l 3 (prolate sources) and the vertex (1,0.5)(pink dot) is associated to cavities which have l 1 =l 2 <<l 3 .The oblique side, which connects the top vertex with the bottom vertex of coordinates (0.33,0.33) (yellow dot), represents cavities with l 1 <l 2 =l 3 (oblate sources) and indeed for l 1 <<l 2 =l 3 the rectangular parallelepiped cavity can be identified with a tensile square crack.

Figure 9 .
Figure 9.The vectors show how the moment ratios change for each source configuration increasing the value of N from 1 to 15 for a TS (panel a) and for a MS (panel b) sources.

Figure 10 .
Figure 10.In the space of coordinates (M 3 /M 1 , M 2 /M 1 ), we compare the domain proper to ellipsoid sources (panel (c)) with the domain proper to TS sources (panel (a), N ≥1) and with the domain proper to MS sources (panel (b), N>1).In panel (d), TS model, and in panel (e), MS model, the area with single line pattern represents the extension provided by rectangular parallelepiped cavities to the ellipsoids domain (crossed lines pattern).

Figure 11 .
Figure 11.Total volume increment DV m predicted for a cubic cavity using the TS and the MS models (panel a).In panel (b), we show how DV m MS − DV mTS tends to reach a constant non vanishing value as N increases (~10%).We consider a cubic cavity buried in an elastic space (n=m=4 GPa) with edge l=100 m and overpressure v T i =10 MPa.

Figure 12 .
Figure 12.The external volume increment DV out (panel a) and the internal volume increment DV in (panel b) predicted with the TS and the MS models are shown as function of the parameter N. DV in I = V 0 DV -K = 1500 m 3 represents the theoretical value for the internal volume increment of a cubic cavity buried in an elastic space (n=m=4 GPa) with edge l=100 m and overpressure v T i =10 MPa.In panel (b), with N>>1, the relative difference between DV in MS and DV in TS tends to ~30%.

Figure 14 .
Figure 14.The value of the ratio DV in /DV m predicted by the TS and the MS models for a cubic cavity tends to different constant values as N increases and both are significantly different from the value DV in /DV m =4/9 proper to a spherical pressurized cavity; the difference is mostly due to DV out for MS sources.

-
S_strain[ ] : array with the 6 strain components due to the strike-slip component -D_strain[ ] : array with the 6 strain components due to the dip-slip component -T_strain[ ] : array with the 6 strain components due to the tensile component

Table 1 .
Relations between the edges l 1 ,l 2 ,l 3 of the rectangular parallelepiped cavity and M 3 /M 1 , M 2 /M 1 moment ratios valid for TS and MS sources if N=1.

Table 2 .
Cubic parallelepiped cavity (l 1 =l 2 =l 3 =100 m) in a half-space (n=m=4 GPa) with h=5 (source depth d=0.5 km).An equal overpressure v i T =10 MPa acts in the center of each face.

Table 3 .
The presence of the inner boundary element (Inner Element, IE) permits the correct identification of the translation of the source interior.The results refer to the same values of the model parameters used in Table2.