Analytical Solutions for Gravity Changes Caused by Triaxial Volumetric Sources

Volcanic crises are often associated with magmatic intrusions or the pressurization of magma chambers of various shapes. These volumetric sources deform the country rocks, changing their density, and cause surface uplift. Both the net mass of intruding magmatic fluids and these deformation effects contribute to surface gravity changes. Thus, to estimate the intrusion mass from gravity changes, the deformation effects must be accounted for. We develop analytical solutions and computer codes for the gravity changes caused by triaxial sources of expansion. This establishes coupled solutions for joint inversions of deformation and gravity changes. Such inversions can constrain both the intrusion mass and the deformation source parameters more accurately.

In this study, we use the Okubo (1991) expressions to derive analytical solutions for the gravity changes associated with the point CDM. We show how gravity changes due to point and finite ellipsoidal sources can be calculated by using the point CDM. We compare the point CDM gravity changes with the Hagiwara (1977) and Trasatti and Bonafede (2008) solutions. Finally, we elaborate on the potential of the model for coupled inversions of surface displacements and gravity changes.

Methods
Deformation-induced gravity changes are usually expressed as the sum of contributions due to deformation in the source region and country rocks and the surface uplift. Here, we adopt a decomposition scheme compatible with the point CDM formulation. We assume a homogeneous, isotropic elastic half-space as a model for the Earth's crust. We denote the Poisson's ratio, shear modulus, and bulk modulus in the medium by ν, μ, and K, respectively. We adopt a right-handed xyz Cartesian coordinate system with the origin at the free surface and the z axis pointing upward. By "gravity change" we refer to the change in the absolute value of the gravity vector's z component. Thus, a positive mass change (mass increase) below a gravimeter leads to a positive gravity change (gravity increase).

Gravity Changes Caused by Magma Chamber Pressurization
As an example, suppose that magma degassing pressurizes a magma chamber (Figure 1). We assume that the exsolved gases all gather at the interface between the chamber walls and the degassed magma, forming a  Equation 1). The country rocks are subjected to positive dilatation/density decrease (light gray and white contours) and negative dilatation/density increase (dark gray and black contours). Thick black contour marks zero dilatation. The gravity station (black triangle) has been subjected to gravity change δg and vertical displacement u v .

of 12
shell-shaped cavity. The outward expansion of the chamber walls and inward compression of the magma lead to the oppositely signed chamber volume change, δV c , and magma volume change, δV m , respectively. The total volume created by the expansion-compression process-namely the interface volume change, ΔV int -is given by where V c = V + δV c and V m = V + δV m are the chamber volume and magma volume in the deformed state, respectively, and V represents both the chamber volume and magma volume in the undeformed state. The chamber expansion also uplifts the surface and generates a strain field, ϵ ij , in the surrounding rocks. This changes the density of the rocks by δρ r = −ρ r ϵ kk , where ρ r is the rock density in the undeformed state and ϵ kk = ϵ 11 + ϵ 22 + ϵ 33 is the volumetric strain or dilatation-a positive dilatation reduces the density (see Figure 1). Similarly, the magma density change, δρ m , due to the compression is related to the magma compressibility, β m , through δρ m = ρ m β m δp, where ρ m is the magma density in the undeformed state and δp is the pressure change in the chamber (Rivalta & Segall, 2008). Provided that β m and δp are known, we have Since we consider the created volume ΔV int as void, the density change in the δV c and δV m portions is −ρ r and −ρ m , respectively. Similarly, uplift, or subsidence, at the Earth's surface will either fill void space or create a void space. So, the other zone of substantial density change is the Earth's surface, where areas of uplift and subsidence are subjected to density changes +ρ r and −ρ r , respectively.
The same deformation-induced density changes exist if instead of exsolved gases, the interface cavity is formed by, and filled with, the intrusion of some external fluids. In such case, the interface cavity is filled with a net mass where ρ int is the intrusion density.
The magma chamber expansion leads to a vertical displacement, u v , and the following gravity change contributions for each observation point at the surface: 1. Δg β , due to density change δρ m in the magma volume in the deformed state, V m , 2. Δ m , due to density change −ρ m within the δV m volume, 3. Δ c , due to density change −ρ r within the δV c volume, 4. Δ , due to density changes δρ r throughout the country rocks, 5. Δg SM , due to the presence of the displaced surface mass layer with density +ρ r , 6. Δg FA , due to the free air change in gravity associated with u v , 7. Δg ΔM , due to the added intrusion mass ΔM that leads to density change ρ int within the interface cavity, for a total surface gravity change of Δg ΔM , also known as residual gravity, can be used to constrain ΔM (see Battaglia et al., 2008). However, this requires all the other terms in Equation 5 to be quantified first. At each station, δg and u v can be determined through repeated gravity and deformation measurements, respectively. Then we have where γ ≃ −0.3086 mGal/m is the free air gradient, and where G is the gravitational constant. Note that Equation 7 uses the Bouguer plate approximation and is valid for flat topographies. The other terms in Equation 5 can be estimated only by using a deformation model for the chamber pressurization. Note that Equation 5 is valid for sources both in the near field and the far field. In the following, we first introduce an analytical point-source model, which can be applied to sources in the far field, and show that in this case Equation 5 can be simplified. Next, we present a semi-analytical finite-source solution and elaborate on the issues that may limit its applicability to near-field problems.

The Far Field Approximations
The far field gravity changes caused by the intruded fluid mass can be calculated through a point-mass approximation as where d is the depth to the center of the chamber and r is the distance between the center of the chamber and the surface observation point. This approximation can also be applied to the far field gravity changes caused by the other density changes in the chamber as The conservation of the initial magma mass in the chamber implies δρ m V m = −ρ m δV m , which together with Equation 9 yields Note that for shallow finite sources Equation 10 does not necessarily hold, as mass redistribution within the chamber may lead to measurable gravity changes. The far field form of Equation 5 can now be written as which expresses the surface gravity changes associated with a deep pressurized chamber as the sum of contributions due to displaced mass at the chamber walls ( Δ c ) , volumetric strain in the host rocks ( Δ ) , displaced mass at the Earth's surface (Δg SM ), and the vertical displacement of gravity stations (Δg FA ), superimposed on the mass change contribution (Δg ΔM ).
Note that Equations 1-11 hold for any chamber shape and boundary conditions on the chamber walls.

Gravity Changes Caused by the Point CDM
The point CDM represents the far field of triaxial sources of expansion with arbitrary spatial orientations . The point CDM is composed of three mutually orthogonal point tensile dislocations (see Figure 2a) constrained to either expand or contract together. The strength of each point tensile dislocation is determined by its potency, defined as the product of dislocation surface area and opening (Aki & Richards, 2002;Nikkhoo et al., 2017, see also Appendix A). The point CDM has 10 parameters: three location coordinates, three rotation angles, three potencies specifying the expansion intensity along the three principal axes of the source, and Poisson's ratio, ν. The total potency of the point CDM, denoted by ΔV, is the sum of the three potencies. ΔV has the units of volume but it is not a physical quantity. Rather, it is a measure of the source strength and it holds ΔV = ΔV int , provided that K m = K, where K m = 1/β m is the bulk modulus of magma.
Triaxial sources of differing shapes, but identical far field deformation, have the same point CDM representation and thus, the same ΔV. However, in order to have the same δV c , these sources must also have identical shapes (except for ν = 0.5, which leads to ΔV = δV c ). For example, the uniformly pressurized cuboidal and ellipsoidal chambers in Figure 2 have the same potencies but their volume changes are different. Analytical expressions relating ΔV to δV c are available for ellipsoidal sources from Eshelby (1957). For uniformly pressurized ellipsoids, we have : Recalling that = 2 (1+ ) 3(1−2 ) and that for a spherical source of radius a the total volume and volume change are Sph = 4 3 3 and Sph c = 3 , respectively, Equation 12 becomes which for ν = 0.25 leads to Δ Sph = 1.8 Sph c (see also Aki & Richards, 2002;Bonafede & Ferrari, 2009;Ichihara et al., 2016). Gravity changes caused by point tensile dislocations can be calculated through the Okubo (1991) analytical expressions (Appendix A). By superimposing the gravity changes associated with three mutually orthogonal point dislocations (Equation A1), we derive the analytical gravity changes associated with the point CDM as where Δg ΔV is the interface cavity contribution (white space in Figures 2b and 2c) and Δg MD is the contribution due to the medium dilatation both inside and outside the source (gray space in Figures 2b and 2c). Noting that Δ Δ = Δ c + Δ m and Δ MD = Δ + Δ and using Equation 10 we have from which it follows that the δg from Equation 14 and the δg from Equation 11 are equivalent. Therefore, the point CDM can be used to compute the effects of deformation on gravity change and thus estimate the mass change ΔM.

Gravity Changes Caused by Point and Finite Pressurized Ellipsoidal Cavities
For any point ellipsoidal model after Davis (1986), there is an equivalent point CDM related to the elastic parameters of the medium and the ellipsoid semi-axes and pressure change through the Eshelby (1957) tensor (see Nikkhoo et al., 2017). Thus, Equation 14 also holds for point ellipsoidal sources. By calculating δV c for ellipsoidal cavities Δ c (Equation 9) and thus, Δ (Equation 15) can be determined for ellipsoidal sources.
Assume that a point CDM with potencies (Δ , Δ , Δ ) represents the far field of a pressurized ellipsoidal cavity with semi-axes ( ) . Then, a set of N point CDMs with potencies (Δ ∕ Δ ∕ Δ ∕ ) , uniformly distributed within the ellipsoid (see Figure 2d), approximates the near field deformations of the pressurized cavity (Amoruso & Crescentini, 2011;Amoruso et al., 2008;Davis, 1986;Eshelby, 1957;Segall, 2010;Yang et al., 1988). Provided that N → ∞, this procedure leads to an accurate solution, unless the cavity is immediately below the free surface (Amoruso & Crescentini, 2011;Segall, 2010;Yang et al., 1988). Similar accuracies can be achieved by using the finite Ellipsoidal Cavity Model (finite ECM) after Nikkhoo and Rivalta (2022), which uses a smaller number of point sources with depth-dependent spacing and strengths. By incorporating the expressions for the point CDM gravity changes in these configurations, we derive new solutions for the gravity changes caused by a finite pressurized ellipsoidal cavity. While the finite ECM is more accurate than the point CDM in modeling shallow pressurized ellipsoidal cavities, it is still an approximate solution for both deformation and gravity change calculations. Similar to the Yang et al. (1988) solution, the finite ECM provides excellent accuracies in the limit that the source dimensions are small compared to its depth (see Nikkhoo & Rivalta, 2022, for further details). Hagiwara (1977) derived closed-form expressions for the gravity change contributions caused by the Mogi (1958) source, later used to validate analytical (Okubo, 1991) and numerical solutions (Currenti et al., 2007(Currenti et al., , 2008Trasatti & Bonafede, 2008).

Comparisons With Other Gravity Change Solutions
An isotropic point CDM is equivalent to the Mogi (1958) model (Bonafede & Ferrari, 2009 which are equivalent to the Hagiwara (1977) expressions (see also Okubo, 1991;Rundle, 1978;Savage, 1984;Walsh & Rice, 1979). This validates the gravity change solution for the point CDM in the case of point spherical cavities. As proved by Walsh and Rice (1979), the sum of the three terms in each set of Equations 16 and 17 vanishes. Note also that, for any point CDM, if ν = 0.5 then Δ MD = Δ = 0 .
We now show that the gravity change solutions for the point CDM also provide a basis for rigorous benchmarking of numerical solutions. We use the point CDM and the finite ECM to calculate the surface displacements ( Figure 3a) and gravity changes ( Figure 3b) associated with the Trasatti and Bonafede (2008)  There is also a good agreement between the gravity changes from all approaches (Figure 3b). The maximum differences between Δ c , Δ , Δg SM , and Δg from the finite ECM and point CDM are ∼6%, ∼9%, ∼9%, and ∼6%, respectively. Since the cavity in this example is relatively deep, the finite ECM calculations are very accurate. Thus, in this particular case, the subtle differences between the finite ECM and the FEM gravity change contributions mostly reflect the errors in the FEM vertical displacements and cavity volume change. The largest difference between the Trasatti and Bonafede (2008) and the other solutions is slightly above 1 μGal, which is more than double the error that Trasatti and Bonafede (2008) estimated by comparison with Hagiwara (1977). This suggests that comparing numerical models with the solution for spherical cavities alone may lead to underestimated errors for the numerical models. Dieterich and Decker (1975) showed that different source shapes produce almost indistinguishable uplift patterns if the source depths are appropriately adjusted. However, the associated horizontal displacements will be completely different. The implication is that in order to constrain all source parameters reliably, horizontal and vertical displacement data must be inverted together. Similar to horizontal and vertical surface displacements, the deformation-induced gravity changes depend on the deformation source parameters. Thus, gravity changes can potentially help better constrain them (Trasatti & Bonafede, 2008).

Implications for the Retrieval of Deformation Source Parameters
We use the point CDM to simulate the radial and vertical displacements and the gravity changes associated with three different axially symmetric deformation sources: a horizontal sill, an isotropic source, and a prolate source (see Figure 4). For all sources, ΔM = 0. The source depths in Figure 4a lead to similar vertical displacements (Figure 4c) but distinct horizontal displacements ( Figure 4d) and distinct gravity changes (free air contribution removed; Figure 4b). Adjusting the source depths differently (Figure 4e) such that the horizontal displacements match (Figure 4h), leads to distinct vertical displacements ( Figure 4f) and distinct gravity changes (Figure 4g). This implies that, from a theoretical perspective, gravity changes may also help to better constrain the deformation  source parameters, beside the mass changes. In practice, however, if ΔM ≠ 0, gravity changes may be dominated by Δg ΔM and thus, depending on the signal-to-noise ratio of the data, the Δg curves (Figures 4b and 4f) may become indistinguishable.

Discussion
Volcanic gravity changes caused by the net mass of intruding magmatic fluids and the induced host rock deformations may have comparable magnitudes to those of hydrological origin, such as changes in the water table.
Such hydrogravimetric disturbances can be estimated by employing hydrological monitoring and modeling techniques (Battaglia et al., 2003(Battaglia et al., , 2006Creutzfeldt et al., 2010;Kazama et al., 2015;Lien et al., 2014;Van Camp et al., 2010) or by analyzing time-lapse gravity data (Güntner et al., 2017). Thus, the mass of intruding fluids at volcanoes can be inferred reliably once such effects are removed.
New generation, low-cost, and accurate gravimeters might soon provide gravity measurements at an unprecedented spatiotemporal resolution (Carbone et al., 2017(Carbone et al., , 2020. Permanent networks provide opportunities for new insight on magmatic plumbing systems (Battaglia et al., 2008;Carbone et al., 2019). One main challenge associated with these developments is to perform both detailed Bayesian inferences for in-depth understanding of the volcano, and rapid inversions for hazard assessment and early warning.
The available FEM gravity change models can incorporate various chamber shapes (Currenti, 2014;Currenti et al., 2007Currenti et al., , 2008Trasatti & Bonafede, 2008), the Earth's surface topography (Charco et al., 2009;Currenti et al., 2007), crustal density and material heterogeneities (Currenti et al., 2007(Currenti et al., , 2008Trasatti & Bonafede, 2008;Wang et al., 2006), viscoelasticity of the Earth's crust (Currenti, 2018), self-gravitation effects (Charco et al., 2005(Charco et al., , 2006Fernández et al., 2001Fernández et al., , 2005, and magma compressibility (Currenti, 2014). Beside difficulties in implementing the FEM, such as meshing issues, this powerful method is computationally too demanding to be used for detailed inverse modeling. In contrast, the point CDM is a half-space model, but has already proven to be suitable for exploring the parameter space in both detailed Bayesian inferences (see Lundgren et al., 2017) and rapid and unsupervised inversions of deformation data (see Beauducel et al., 2020). The gravity change solutions for the point CDM, which we provide here, extend this potential to joint inversions of surface displacements and gravity changes. Volcanic deformation sources are often deep or far enough from the observation point to be treated as far field sources. The point CDM can provide a first order solution, which can be later improved by more sophisticated numerical models. Some complexities such as layering or viscoelasticity can be accounted for (Amoruso et al., 2008) by using appropriate Green's functions for point dislocations (Okubo, 1993;Sun & Okubo, 1993;Wang et al., 2006). Besides, theory errors, arising from ignoring real Earth complexities, can be estimated in terms of noise covariance matrices within a Bayesian framework (see Duputel et al., 2014;Minson et al., 2013;Vasyura-Bathke et al., 2021).
Finite pressurized ellipsoidal cavities can be approximated by a set of point CDMs uniformly distributed in the cavity volumes. With a high number of point CDMs, this approach can be used for benchmarking numerical models. An alternative solution is the finite ECM after Nikkhoo and Rivalta (2022), which provides comparable accuracies for a lesser number of point CDMs. The finite ECM is very fast, and thus, provides a practical way for performing coupled inversions of surface displacements and gravity changes.
It is important to recall that for ellipsoidal deformation models in the half-space, including the finite ECM and the Yang et al. (1988) spheroid, the full-space expressions are used to calculate δV c (Amoruso & Crescentini, 2009). While this approximation may often be acceptable for deformation studies, it may lead to large errors in gravity change calculations involving shallow finite sources. This warrants future systematic comparisons with numerical models in order to quantify the associated error.
Deformation-induced gravity changes may be substantial (see Figure 3b) and should be accounted for in joint inversions of surface displacements and gravity changes. Provided that coupled models are employed for such inversions, the gravity changes may be exploited to better constrain the deformation source parameters besides the mass change. How practical this may be, depends on the observation uncertainties and the signal-to-noise ratio. We will explore this feature in future studies.
Coupled inversions of surface displacements and gravity changes constrain the deformation source parameters and the intrusion mass without making any assumption on the properties of the intruding fluid. The intrusion density can be estimated from the inferred mass only if the interface volume change, ΔV int , is known (ΔV int should not be mistaken for the chamber volume change δV c ). It can be shown from Equations 2 and 3 that the determination of ΔV int requires knowledge of the fluid compressibility. This shows that unlike mass change estimates, the estimates of the intrusion density are prone to large uncertainties.

Conclusions
1. Surface gravity changes are sensitive to both the intruding fluid mass and the deformation-induced surface uplift (subsidence) and country rock dilatation. Due to this coupling between the gravity changes and host rock deformations, gravity changes can also be used to constrain deformation source parameters, namely the location, spatial orientation, and potency of triaxial source models for expanding reservoirs. 2. We provide analytical solutions and MATLAB codes for surface displacements and gravity changes caused by both the point CDMs, a model for triaxial sources of expansion and the finite ECM, a model for ellipsoidal sources of uniform pressurization. 3. While modeling gravity changes caused by shallow sources, it may be necessary to account for the mass redistribution within the source. This issue and also the inherent error in δV c for half-space solutions may limit the applicability of the finite ECM. 4. The analytical solutions presented here can be used to validate new numerical gravity change models. Such validations should ideally consider various source depths and aspect ratios. 5. By using the point CDM and the finite ECM, coupled inversions of surface displacements and gravity changes can now be performed.