A New Analysis of Caldera Unrest through the Integration of Geophysical Data and FEM Modeling: The Long Valley Caldera Case Study

: The Long Valley Caldera, located at the eastern edge of the Sierra Nevada range in California, has been in a state of unrest since the late 1970s. Seismic, gravity and geodetic data strongly suggest that the source of unrest is an intrusion beneath the caldera resurgent dome. However, it is not clear yet if the main contribution to the deformation comes from pulses of ascending high-pressure hydrothermal fluids or low viscosity magmatic melts. To characterize the nature of the intrusion, we developed a 3D finite element model which includes topography and crust heterogeneities. We first performed joint numerical inversions of uplift and Electronic Distance Measurement base-line length change data, collected during the period 1985–1999, to infer the deformation-source size, position, and overpressure. Successively, we used this information to refine the source overpressure estimation, compute the gravity potential and infer the intrusion density from the inversion of deformation and gravity data collected in 1982–1998. The deformation source is located beneath the resurgent dome, at a depth of 7.5 ± 0.5 km and a volume change of 0.21 ± 0.04 km 3 . We assumed a rhyolite compressibility of 0.026 ± 0.0011 GPa − 1 (volume fraction of water between 0% and 30%) and estimated a reservoir compressibility of 0.147 ± 0.037 GPa − 1 . We obtained a density of 1856 ± 72 kg/m 3 . This density is consistent with a rhyolite melt, with 20% to 30% of dissolved hydrothermal fluids.


Introduction
The Long Valley Caldera (LVC), located in east-central California on the western edge of the Basin and Range Province and at the base of the Sierra Nevada frontal fault escarpment, is an east-west elongated oval depression formed by the eruption of the Bishop Tuff, 767,100 ± 900 years ago ( Figure 1). Beginning in the late 1970s, the caldera entered a period of unrest, without any eruptions, that continues to the present time (e.g., Figure 3 in [1]). The unrest episodes include recurring earthquake swarms beneath the South Moat Seismic Zone (SMSZ) and the Sierra Nevada (SN) block, accelerated inflation of the central Resurgent Dome (RD), variations in the geothermal system and gas emissions around the flanks of Mammoth Mountain (MM) on the southwest margin of the caldera ( [1] and references therein). The U.S. Geological Survey (USGS) began an intensive effort to monitor the unrest in LVC between 1975 and 1983 with the setup of new leveling lines in 1982, a two-color Electronic Distance Measurement (EDM) network in 1983, trilateration arrays in 1979, a dense seismic network in 1982, and a high-precision gravity network in 1982. Continuous Global Positioning System (GPS) measurements have been made since 1993. Both ground-based and space geodesy (including satellite interferometry) observations reveal a consistent radial and upward deformation pattern, centered at the RD and decreasing with radial distance (e.g., [2][3][4][5][6][7][8]).
The processes driving the unrest at LVC remain unclear, with the main likely source of unrest being either a magmatic intrusion into the upper crust [1], or pulses of highpressure hydrous fluid intrusion into the upper crust [18]. Geologic and petrologic evidence support the hypothesis that the LVC rhyolitic magmatic system is moribund and that the magma body that fed the caldera-forming eruption may now be in the final stages of crystallization. The most recent eruptions along the Inyo Craters/Mono Domes chain and Mammoth Mountain have been fed by a different magmatic system. All of the eruptions inside the LVC have been rhyolitic, with the most recent eruption ~100 ka in the west moat. There has been no eruption on the resurgent dome over the last 500 ka. No significant seismicity and no emission of CO2 or other magmatic gases has been recorded beneath the resurgent dome. Finally, the drilling of the resurgent dome found temperatures of only 100° at a depth of 3 km [18]. On the other hand, several pieces of geophysical evidence point to a possible magma intrusion as the cause of the present unrest. Multiple seismic imaging studies using different techniques (e.g., teleseismic tomography and fullwaveform ambient noise tomography) highlighted large low velocity zones in the middle and lower crust, which have been interpreted as evidence of the presence of a partial melt. Different geodetic data (both ground-and satellite-based) measuring deformation at LVC since 1979, have not recorded any substantial deflation episodes yet. This might instead be expected if the inflation involved the injection of hydrothermal fluids with poroelastic swelling followed by diffusion, as observed at other calderas, such as Yellowstone and Campi Flegrei [1].
While ground deformation can provide insights about volume changes in the underground reservoir, it cannot constrain the mass of the intrusions and therefore discriminate between magma and hydrous fluid intrusion. Combined deformation and gravity measurements can be used to infer the density of the intrusive fluids and better define the source of unrest [19][20][21][22][23][24][25][26][27]. Given the density difference between silicate melts (~2300 kg/m 3 ) and hydrothermal fluids (~800 kg/m 3 , [28]), density estimates can, in principle, be used to distinguish between these two possible sources of caldera unrest.
Gravity measurements at LVC have been conducted yearly between 1980 and 1985 and repeated in 1998 and 1999 [28]. In this period, the RD experienced a quasi-steady uplift, with accelerated phases in 1989-1990 and 1997-1998, when the most rapid deformation occurred (e.g., [4,13]). These data have been analyzed, together with different kinds of deformation records (EDM, leveling, GPS, InSAR) in different studies using analytical models and considering increasing complexities, from point source to tilted finite ellipsoidal source, from homogeneous to vertically layered elastic half-space [12,[28][29][30][31]. The results of these studies suggest that gravity data are more compatible with the addition of a magma intrusion than pulses high-pressure hydrous fluids.
In this paper, we consider the 1982-1999 unrest period. This time has the best and most complete gravity dataset. We perform numerical computations based on the finite element method (FEM), exploring the effect of topography and realistic medium heterogeneities on the parameters (e.g., location, depth, density) of the source of unrest. We first invert EDM and leveling data from 1985 to 1999 in order to constrain the location, depth, and geometry of the unrest source. We then use the inferred source to model the deformation and gravity changes between 1982 and 1999, and to compute the source volume change, and density.
In Section 2, we present the methods including data, model setup and model computations; in Section 3, we show the results; in Section 4, we discuss our findings and conclusions.

Data
We used the data from the Long Valley Caldera GIS Database ( [13]; https://doi.org/10.3133/ds81, accessed on 15 February 2021). The database includes extensive geologic, monitoring, and topographic datasets from studies conducted in Long Valley caldera between 1975 and 2001. The unrest is investigated using three sets of data: baseline length changes (an approximation of horizontal deformation) from two-color EDM, vertical deformation from a combination of GPS and leveling, and gravity changes.
The two-color EDM network consists of two sets of seven baselines. The first set is formed by the sites Hot, Knol, Krak, Mine, Shark, Sher and Till, observed from the central monument CASA (green triangles in Figure 1a). The second set is formed by the sites Bald, Dead, Knob, Krak, Micr, Mike and Sage, observed from the LKT monument (cyan triangles in Figure 1a). Measurements at these baselines span the 1985-1999 inflation period, which is included in the targeted time in this work (Supplementary Materials, Table S1). The methods used to extract the displacement and its error for each of the baselines are described in [11,17]. The EDM deformation data that were used are from [13].
Vertical deformation (uplift) measurements were taken during different leveling surveys along the 65-km-long line along Hwy 365 from Tom's Place to Lee Vining, and along several other routes within the caldera, and are obtained by combining leveling and GPS data ( Figure 1a). Complete leveling of the caldera occurred each summer from 1982 to 1986, and in 1988 and 1992. In 1999, reference [13] occupied 44 leveling benchmarks with GPS to bring up to date the direct measurement of vertical deformation. The data sets employed in this work consist of the 44 benchmarks with leveling and GPS for the period 1985-1999 [13], and 34 benchmarks with leveling and GPS for the period 1982-1999 [28] (red circles respectively in Figure 1a,b; Supplementary Materials, Tables S2 and S3). The benchmark C916, located near Lee Vining (Mono Lake), is the elevation datum for the vertical deformation. The standard error for each elevation difference is calculated according to [13].
The Long Valley caldera gravity monitoring network is centered near Tom's Place (the primary reference station) and extends from the Sierra Nevada west of Lee Vining, CA, southeastward to a station in the White Mountains east of Bishop, CA [32]. Gravity data (gravity changes, noise from the water table and gravity corrected for the water table and free-air effect) are from [28], Supplementary Materials, Table S4. In Section 2.4, we employ the gravity data corrected for water table and free air contribution to estimate the density of the intrusion.

Model Setup
We develop a three-dimensional (3D) numerical model using the finite element method (FEM) and the software COMSOL Multiphysics (www.comsol.com, accessed on 15 February 2021). The geometry refers to a Cartesian reference system and is composed of a domain of 120 km × 130 km. The model is 40 km deep (up to the Moho depth in the area [33,34]), with zero depth corresponding to the sea level. The chosen size represents a crust portion which includes the LVC and a significant part of its surroundings ( Figure  2a).
Inside the domain, we assume the existence of an internally pressurized ellipsoidal prolate cavity that we invert for its location, dimensions, and overpressure ( Figure 2b  Pressurized cavities of simple geometry can mimic/approximate the crustal stress field produced by the actual source. None of these geometries reproduced an actual source. The actual deformation source, beneath the resurgent dome, is probably a network of fractures filled with fluids (or magma) ascending from the crystallizing Pleistocene pluton below [18].
We explore three different crust configurations. A homogeneous elastic domain with flat stress-free top surface (labeled HF) representing the average altitude of the area (~2300 m a.s.l.), a homogeneous elastic domain with topography (labeled HT), and a fully heterogeneous elastic domain with topography (labeled HeT). The topographic surface is generated by using the STRM digital elevation model (DEM) from the USGS Earth Explorer [35], resampled at 600 m resolution. Material heterogeneities (density, bulk modulus, shear modulus and Poisson's ratio; Figure 2c-f) are obtained from pressure (VP) and shear (VS) wave velocity distributions. Shear wave velocities VS are from [36] while pressure velocities VP are calculated from VS using a VP/VS ratio of 1.75 for the SN block [37] and of 1.79 elsewhere, with a gradient d(VP/VS) of 3% [33,38]. VP, VS velocities are converted into Poisson's ratio (ν), density (ρ) and dynamic Young's Modulus (E) using the equations in [39]: However, to properly represent the medium strain rate in a quasi-static condition, we need to refer to quasi-static mechanical properties. Laboratory tests [40,41] show, in fact, that for lithostatic pressures in the range 1-3 kbar (3.8 to 11.5 km depth), the ratio between quasi-static and dynamic bulk modulus Ks/Kd for granite is different from 1 and can vary between 0.5 (at 0.09 kbar-0.4 km depth) to 0.9 (at 3 kbar-11.5 km depth). For the range 0-3 kbar of lithostatic pressure (equivalent to the distance between the top surface and 11.5 km depth), we calculate the dynamic bulk modulus from VP, VS values and multiply it by the Ks/Kd ratio values from [40] at the corresponding lithostatic pressure (depth) level to estimate the quasi-static bulk modulus. The relation between quasi-static and dynamic mechanical properties is empirical and depends on several factors including stress state and stress history [40,41], however our approach leads to a better characterization of the material response with respect to what can be obtained using pure dynamical properties. An interpolation function guarantees a smooth transition between different lithostatic pressure levels. At a greater depth, where the lithostatic pressure is higher than 3 kbar, we assume a Ks/Kd ratio of 1. From the quasi-static bulk modulus, we can calculate the quasi-static shear modulus. Since Poisson's ratio is not expected to change significantly [42], we can retain its dynamic value. Crust properties are summarized in Table 1 and represented in Figure 2c-f. In terms of boundary conditions, the model bottom is fixed, the top surface is stressfree while at the lateral boundaries we apply a roller condition (no displacement in the direction normal to the boundary). An infinite element condition, set at the lateral and bottom boundaries, simulates the far-field, and guarantees that the displacement vanishes at a very far distance from the original geometric size, thus avoiding any boundary effects. We prescribe a parametrized overpressure on the boundaries of the ellipsoidal cavity. The model domain is meshed with tetrahedral elements while the source boundaries and the top surfaces are discretized with triangular elements. Automatic adaptive mesh refinement tests are carried until an optimal performance is found without further variation of the output. The mesh for the whole domain is shown in Figure 3. We validate the FEM model for the deformation by performing a benchmark calculation for a vertical prolate ellipsoid in a flat homogeneous domain and comparing the numerical results against the analytical solution by [43], see Appendix A.

Inverse Modeling of Deformation Data
We use the FEM model described in Section 2.2 to perform numerical inversions of EDM baseline changes and uplift (leveling-GPS). Inversions are performed in two steps.
In the first step, we jointly invert the EDM and leveling-GPS data for the period 1985-1999 [13] and for each model configuration (with/without topography or with topography and heterogeneities) we infer the best-fit source dimensions, position, and overpressure. In the second step, we keep the deformation source stable in size and location and further optimize only for the source overpressure by performing a second inversion of the leveling-GPS data for the period 1982-1999 [28]. This second step provides the source parameters needed to model the gravity changes.
Inverse modeling is performed using the Nelder-Mead solver [44] and coupling the structural mechanics and the optimization module in COMSOL Multiphysics. The link between the data to invert and the source parameters is built by setting up the objective function [45] and = Mi are the modeled data, Di the observed data, Ei the observation error, Wi the weights and the index i relates to each benchmark. The inversion goal is to minimize the least weighted squares Equation (4).
The inversion of EDM and leveling-GPS datasets for the period 1985-1999, is made considering the following seven parameters of the ellipsoidal source: the semiaxes (Ea, Eb and Ec), oriented along the cartesian x, y and z axis respectively; the horizontal position along the x and y direction (Ex, Ey) with respect to a reference point located at longitude −119°W, latitude 37.5°N; the source vertical position (Ez), with Ez = 0 m at the sea level; the overpressure (ΔP) applied at the source internal boundaries. The first three parameters control the source geometry, the second three control the source position and the last one controls the source overpressure. The source is assumed to be vertical, with a plunge of 90°. For each parameter we assign an initial value (used in the first iteration), and lower and upper bounds to search for reasonable values during the computation ( Table 2). Initial values and ranges of parameters are based on results from previous studies (e.g., reference [7] and references therein). In particular, because of a poor data coverage in the south caldera rim, we constrained the north-south source position Ey to fall within the resurgent dome area. The solver is set to perform a maximum number of 400 model evaluations. This threshold has been chosen by looking at the convergence rate and considering that further evaluations no longer have any significant impact on the value of the objective function. According to [29], the inflation source could be slightly tilted with a dip angle between 91 and 105 degrees. To check whether this was the case for our source, we performed preliminary tests, including additional parameters for a source rotation of ±10 degrees around each of the three cartesian axes. Results showed that in all cases (with/without topography or heterogeneities) the optimal rotation is minimal, <1° around x and y axis and <2° degrees around the z axis. For this reason, we did not include these parameters in our inversions.

Computation of Gravity Changes
The total gravity change recorded at a benchmark during unrest episodes contains the effect of different contributions: (i) the free-air effect, due to the vertical displacement of the benchmarks at the ground surface during unrest; (ii) the water table effect, proportional to the water table level change in the area; (iii) the deformation effect, due to the coupling between elastic deformation and gravity; and (iv) the residual gravity, which depends on the density change related to the introduction of the new mass into the pressurized volume (e.g., [20][21][22]). Furthermore, the estimation of gravity variation is sensitive to model complexities, such as volumetric source geometry, topography, material heterogeneities and fluid compressibility (e.g., [22][23][24][25][26][27]).
The best fit parameters of the ellipsoidal source from the three different crust configurations (HF, HT, HeT; Figure 2), are used to compute the gravity change at the free surface. Following the methodology in [46], we first compute the displacement field from the previously estimated best source models (cf. Section 2.3), and we then solve the Poisson's equation relating the gravity potential (φg) to the change in density distribution (Δρ) caused by the subsurface mass redistribution ∇ = −4 ∆ ( , , ), where G is the gravitational constant. The gravity change can be then computed as δg = −∂ φg/∂ z. In particular, the relation between the gravity potential φg and the density variations can be expressed by the following contributions: where u is the inflation-related displacement field, ρ0 the embedding medium density and ρin is the density of the intruding fluid. Equation (7) gives the gravity contribution δg1 due to the displacement of density boundaries in heterogeneous media, corresponding to the Bouguer correction at the surface in case of flat homogeneous models. Equation (8) gives the gravity contribution δg2ΔV due the displacement of the source boundaries, which implies replacement of the surrounding mass. Equation (9) gives the term δg3, which considers the effect of dilatational/compressional strains in the host rock, while Equation (10) gives the term δg2V which considers the input of material (of density ρin) into the source volume [21]. Equations (7)- (9) can be used to compute the massless deformation contribution to the gravity changes while Equation (10) represents the contribution due to the source mass change.
To numerically solve the Poisson's equations, we modify the model geometry by adding an additional domain with same size of the rock domain (Figure 2b), but made of air (assuming E = 1 Pa, ρ = 1 kg/m 3 , ν = 0.25; Table 1). Furthermore, solving for all contributions to the gravity potential requires the embedded source to be a domain and not a cavity, as done during the inversion of displacements. Poisson's Equations (7) and (8) are solved on the stress-free surface and on the source boundaries, respectively. Poisson's Equation (9) is solved on the domains surrounding the source, and (10) is solved on the source domain.
We validate the FEM model by performing a benchmark calculation for a vertical prolate ellipsoid in a flat homogeneous domain and comparing the numerical results to the analytical solutions by [47] (see Appendix A).
When estimating the gravity changes due to reservoir inflation, it is important to consider that the volume change accommodating for the input of new mass could arise, not only from the expansion of the source wall that deforms the surrounding medium, but also from the compression of the material stored in the reservoir (e.g., [25,[48][49][50]). The relation between the actual volume of the mass intrusion, ∆ , and the volume change from the inversion of deformation data (geodetic volume, cf. Section 3.2), ∆ , is (e.g., [48,51]) where = is the compressibility of the material stored in the reservoir, = is the compressibility of the reservoir due to medium elasticity and reservoir shape, and are the bulk moduli, rV is the volume ratio, and βm is a function of several parameters, like pressure, gas volume fraction, temperature, phenocryst content and source depth (e.g., Table 3 from [52]). Finite element calculations of reservoir compressibility βc as a function of the source geometric aspect ratio = indicates that in our case ≈ , where μ is the shear modulus (see Figure 5 in [53], and Tables 1 and 3). can also be computed as [48]: where ΔP is the overpressure, and V is the source volume before the application of the overpressure (Table 4). Finally, the density corrected for the effect of compressibility (ρcmp) can be computed as:  1 Note that here we indicate the source depth with respect to the stress-free surface and not the sea level, i.e., accounting for the average elevation of LVC area (~2300 m a.s.l.).

Deformation: Best Fit Source
We find that the three different crust configurations (homogeneous, flat elastic halfspace HF; homogeneous elastic domain with topography, HT; heterogeneous elastic domain with topography, HeT) give similar results for the position (Ex, Ey), depth (Ez), geometric aspect ratio (A), and pressure change (ΔP) of the source ( Table 3).
The nonlinearity of the inverse problem makes the evaluation of uncertainties difficult; nonlinear error propagation is a difficult problem to address, COMSOL does not have a feature that allows extraction of the covariance matrix, and the model covariance matrix may not give a good estimate of the uncertainties [54]. A solution could be to employ a Monte Carlo method. Unfortunately, this method requires each model to be run thousands of times. We employ the result from the inversions (350 to 400 runs) to mimic a Monte Carlo method and obtain an estimate of the uncertainties of the source parameters. We then propagate the errors to the density results (see Appendix B).
The source is moved about 1 km eastwards (Ex), 2.4 km southwards (Ey) and 600 m deeper (Ez), with respect to its starting position and starting depth. No major differences can be seen between a homogeneous flat crust (HF), homogeneous crust with topography (HT) or heterogeneous crust with topography (HeT). The source size is similar for the HF and HT crust models, while we can observe in the HeT crust model a significant reduction of about 1/3 of all semiaxes. The source shapes for each crust model before and after inversion are showed in Figure 4. Although the main source parameters are similar for the three different models (Table 3), material heterogeneities make a difference, especially in the estimate of the absolute volume of the source (see Table 4).  Figure 5 shows the total displacement (combination of vertical and horizontal displacements) at the free surface for each model configuration. For the HF and HT cases we can clearly observe two lobes with higher displacement northwards and southwards of the source with the topography playing a minor damping role and a slight clockwise rotation of the northern lobe. When heterogeneities are introduced (HeT), the southern lobe disappears while the northern lobe further rotates clockwise and shows an area of maximum total displacement.   Tables S1 and S2 in Supplementary Materials. Results for EDM show a good agreement between models and observations ( Figure  6a) except for LKT-MIKE baseline which is slightly underestimated in all three crustal models. However, the difference between the models' results and the measurement for LKT-MIKE is within the data uncertainty (thin solid black lines). From the EDM residuals (Figure 6b) we can observe that the fit improves when we add topography (χ 2 decrease by 19% from HF to HT, red and blue arrows) and material heterogeneities (χ 2 decrease by 21% from HF to HeT, red and yellow arrows).  The inversion results for the leveling data ( Figure 7a) show a good agreement between the observed and modeled data at the benchmarks located inside the caldera border (black dashed line). However, the model underestimates the observed uplift by 5-10 cm at the benchmarks located outside the southeastern (SE) caldera border (near Crowley Lake) and by 2-5 cm for the benchmarks located at the northwestern (NW) caldera border (Figure 8). This discrepancy influences the χ 2 value (Figure 7b), which slightly decreases when topography is included (5% decrease in χ 2 from HF to HT) but increases when we add the material heterogeneities (16% increase in χ 2 from HF to HeT). This is probably because, in the HeT model, the material outside the caldera area is stiffer than the material inside it (Figure 2). 3). In this case, we reach a good agreement between the models and uplift for all three crustal models, with residuals within the observation error (Figures 9b and 10). Some discrepancy is still observed for the benchmarks southeast of the caldera border, HeT since the uplift in this area is probably controlled by the deformation along the Sierra Nevada block [55]. In this case, the χ 2 value is reduced by 36% when we add topography (from HF to HT) and by a further 10% when we also add the heterogeneities (46% total decrease χ 2 in from HF to HeT). Observed values and numerical results for the leveling data over the period 1982-1999 are reported in Table S4, Supplementary Materials. The large residuals observed for three benchmarks close to the center of the resurgent dome are from the exploitation of the hydrothermal aquifers by the local geothermal power plant [13]. Other discrepancies are because of motion along faults in the caldera South Moat [4] or the Sierra Block (e.g., Figures 7 and 9) [8,55]. The heterogeneous model (HeT) can better fit the uplift for 1982-1999 than the two homogeneous models (HF and HT; Figure 9).

Density of the Intrusion
The observed gravity change, after free-air and water-table correction, shows a positive anomaly centered on the resurgent dome, with peak amplitude of about 60 μGal ( [28]  HeT and black circles in Figure 11), that suggests mass intrusion into the sub-caldera crust beneath the resurgent dome.
The estimate of the density of the intrusion requires three steps. First, we calculate the gravity changes associated with the deformation of the source and of the surrounding crust (δg1 + δg2ΔV + δg3) by solving (7)-(9), the so-called "deformation effects", see Figure 11. The gravity variations due to deformation effects are substantial in the near-field of the source location, with the highest values at the source-tip location (up to −60 μGal, i.e., comparable, in magnitude, to the observed gravity change after free-air and water-table correction) and decaying to magnitudes <10 μGal at a distance of ~10 km. Figure 11. Comparison between observed gravity changes (1982-1999, Table S4, Supplementary Materials) after the removal of free air and water table effect (black circles with 1-sigma error bars) and the total deformation contribution to gravity changes (δg1 + δg2ΔV + δg3) from the solution of equations (7)-(9) corresponding to different model configurations (colored circles), ordered according to the horizontal distance from the resurgent dome center. HF: homogeneous crust, red circle; HT: homogeneous crust with topography, blue circle; HeT: heterogeneous crust with topography.
We then subtract the deformation effects from the observed gravity changes in order to calculate the residual gravity (black circle with error bars in Figure 12). It is worth HeT noting that the observed gravity changes had been previously corrected for water table noise and the free-air effect (Table S4, Supplementary Materials) (details in [28]). The residual gravity, δg2V, depends on the mass change accompanying the deformation. Following the methodology described in Section 2.3, we solve the inverse problem with Poisson's Equation (10) to obtain the density value for the intruding fluid which best matches the observed residual gravity change δg2V. Numerical results for modeled residual gravity change agree within the measurement errors with most of the residual gravity observations ( Figure 12). Figure 12 shows a good fit to the observed residual gravity for the FEM model of a homogeneous crust with topography (HT). Adding additional information about heterogeneities in the crust does not significantly improve the fit.  Table 4. HF: homogeneous crust, red circle; HT: homogeneous crust with topography, blue circle; HeT: heterogeneous crust with topography, orange circle. Table 4 shows the resulting density values of the intrusion, assuming that the mass change is due to either an incompressible ( ) or compressible ( ) fluid intrusion. The relation between the density of incompressible ( ) or compressible ( ) fluid intrusion is given in Equation (13). Reference [53] and Equation (12) allow for computing the compressibility of the crust surrounding the intrusion from quasi-static elastic properties (Table 1) and the results of the FEM models ( Table 4).
The density of the intrusion depends on magma and reservoir compressibility-Equations (11)- (13). According to [56], the isothermal compressibility for a rhyolite, with a volume fraction of water between 0% and 30%, is 0.026 ± 0.0011 GPa −1 . We estimated a crust compressibility of 0.147 ± 0.037 GPa −1 for the heterogeneous model. Using these values, we obtained a density of 1856 ± 72 kg/m 3 . This density is consistent with a rhyolite melt (no crystals) with 20% to 30% of dissolved hydrothermal fluids.

Summary
We build a 3D finite element model to investigate the source of observed displacements and gravity changes at Long Valley Caldera for 1982-1999. Using the geodetic data available in Long Valley Caldera GIS Database (https://doi.org/10.3133/ds81, accessed on 15 February 2021)-EDM baseline change, uplift, and gravity change measurements-we explore different model configurations starting from a flat and fully homogeneous domain and then adding additional complexities, such as topography and medium heterogeneities. Limits on the coverage of the existing geodetic monitoring network, ambiguities on the interpretation of subsurface distribution of the crust elastic properties, and the nature of non-linear inversion make our models and solutions non-unique. We can improve the bounds on the parameters of the deformation source by employing an appropriate modeling approach.
Since the joint inversion of horizontal and vertical deformation better constrain the geometry of the deformation source, we first invert EDM and leveling data for the period 1985-1999 to infer the size and location of an ellipsoidal source under the caldera. We then optimize the estimate of the source volume and mass change by performing a second set of inversions of leveling and gravity data for 1982-1999. This is the period with the best signal to noise ratio for gravity data.
One advantage of our approach is, not only the inclusion of topography, but also of full heterogeneities (3D). The influence of mechanical heterogeneities in LVC has been considered in other works [30,[57][58][59], but while past works relied on simplified models in terms of geometry (e.g., 2D) or the material heterogeneities distribution (e.g., only vertical), in this study, similarly to [8], we implement both lateral and vertical material heterogeneities. Furthermore, we account for the difference between static and dynamic mechanical properties, since the use of a dynamic bulk modulus for the overall domain would overestimate the medium rigidity at shallow depths (lithostatic pressure < 3 kbar).
To estimate the subsurface mass change of the deformation source, we first estimate the so called "deformation effects" given by (7)-(9); see Figure 11. We then subtract the "deformation effects" from the gravity changes to obtain the residual gravity, since the residual gravity, δg2V, depends on the mass change accompanying the deformation. Finally, we solve the inverse problem with Poisson's Equation (10) to obtain the density value for the intruding fluid which best matches the observed residual gravity change ( Figure 12).
The density of the intrusion will change if the fluid is either incompressible ( ) or compressible ( ), since the density of the intrusion depends on the ratio ⁄ (Equation (13), Table 4). We can estimate either from the shear modulus and source aspect ratio (Table 1; [53]) or from the absolute volume of the source inferred from our numerical analysis (Table 4). Both approaches calculate the same value of for a homogeneous medium. We find the second approach more appropriate for our heterogeneous model.

Conclusions
Gravity data are usually noisier than deformation data (e.g., [28]) but are essential for estimating the density of intrusion, because changes in the gravity potential are related to the changes in density distribution caused by the subsurface mass redistribution. Without gravity data, we cannot obtain information about the nature of the deformation. In this specific case, the major ambiguity is not coming from the errors in the gravity data but from the uncertainty about the appropriate value of magma compressibility. Reference [56] present experimental values of the isothermal compressibility of rhyolite, andesite, and basalt glasses as a function of the volume fraction of water (see Table 4 and Figure 5 in [56]). We assume here the compressibility values for rhyolite, i.e., the main component of erupted magma in LVC [18]. According to [56], the isothermal compressibility for a rhyolite, with a volume fraction of water between 0% and 30%, is 0.026 ± 0.001 GPa −1 . Using our inversion results for source parameters, for the heterogeneous case we estimated a reservoir compressibility of 0.147 ± 0.037 GPa −1 ( Table 4). We therefore obtained a density of 1856 ± 72 kg/m 3 . This density is consistent with a rhyolite melt with 20% to 30% of dissolved hydrothermal fluids.

Acknowledgments:
The authors wish to thank James Hickey (Bristol University) and Gilda Currenti (INGV) for the discussion about the model implementation. Mehdi Nikkoo (GFZ) for supporting the benchmark studies with analytical solutions and for the helpful conversations. The authors thank the Reviewers for the thorough and insightful comments which helped to improve the manuscript and Emily Montgomery-Brown (USGS) for internal review of the manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Validation of Numerical vs. Analytical Solution
To assess the accuracy of the finite element calculations, the absence of edge effects and ensure that the geometry and mesh adopted yield sufficient sensitivity, we validated the FEM, the numerical solution, for a vertical prolate ellipsoid in a homogeneous model without topography against the analytical solutions from [43]. For this comparison, we used a prolate ellipsoid with horizontal semiaxes a = b = 1.5 km and vertical semiaxis c = 3 km. The ellipsoid depth is 7 km and an outward (inflating) overpressure of 70 MPa is applied. We first compare the horizontal and vertical components of displacement on the stress-free surface on ±40 km EW profile ( Figure A1a). Successively, we calculate the gravity changes due to the source deformation and to an intruding fluid with a density of 2700 kg/m 3 using the method described in the main text (Section 2.4) and compare them with the analytical solutions from [47] (Figure A1b,c). For the gravity contribution from the mass change (δg2V), we also checked with the solutions from [60]. We observe a good agreement between the analytical and FEM solutions for both the displacements and the gravity components. Minor differences visible over the first 5 km distance from the source are inside the ranges of uncertainty, which are ~1 mm for horizontal displacements, ~1 cm for vertical displacements and ~10 μGal for gravity changes.