Constraining Magma Reservoir Conditions by Integrating Thermodynamic Petrological Models and Bulk Resistivity From Magnetotellurics

Magnetotelluric (MT) data image the bulk resistivity of the subsurface which can be used to infer magma reservoir conditions beneath volcanoes. The bulk resistivity of magma depends primarily on the melt volume fraction, temperature, and water content. These variables are coupled thermodynamically, yet mixing relations for bulk resistivity implicitly treat them as independent. Here, we use a parameterization of the rhyolite‐MELTS thermodynamic model to constrain relationships between melt fraction, temperature, dissolved water content and bulk resistivity for rhyolitic magmas. This method results in MT interpretations which are (a) thermodynamically consistent at near‐equilibrium conditions, (b) independent of temperature and water content estimates derived from erupted products, and (c) able to consider saturated melts containing a volatile (i.e., aqueous fluid) phase. The utility of the method is demonstrated with three case studies of silicic systems: Mono Basin, Newberry volcano and the Laguna del Maule Volcanic Field (LdMVF). The moderately conductive feature at Mono Basin can be explained by under‐saturated partial melt (6–15 vol%) at <775°C, indicating relatively stable magma storage conditions since the last eruption. However, a relatively resistive feature at Newberry Volcano requires lower temperatures (<750°C) than previous estimates, suggesting that the system has cooled since the last eruption. A conductive feature at the LdMVF cannot be explained by saturated or under‐saturated rhyolitic melt and requires additional conductive phases. These results demonstrate the potential of this new method to reduce uncertainty in MT interpretations and highlight the need for additional coupling strategies between petrology, geophysics, and thermo‐mechanical models to better understand magmatic systems.

The MT method images the bulk resistivity structure of the subsurface and, in volcanic settings, this can be related to melt fraction using laboratory empirical relations and mixing relations. Laboratory studies measure the melt resistivity to develop empirical relations based on parameters such as temperature, volatile content, and composition (e.g., Gaillard, 2004;Guo et al., 2016Guo et al., , 2017Laumonier et al., 2015). The melt resistivity and imaged bulk resistivity are then used to infer melt fraction using a mixing relation (e.g., Glover, 2010;Glover et al., 2000;Hashin & Shtrikman, 1963;Schilling et al., 1997). However, over the last decade, three important challenges have emerged in using MT to derive magma storage conditions beneath volcanoes: (a) nearly all MT interpretations implicitly assume that melt fraction is independent of temperature and volatile content (e.g., Bowles-Martinez & Schultz, 2020;Cordell et al., 2018;Comeau et al., 2016;Hübert et al., 2018;Samrock et al., 2018, etc.), ignoring that these variables are coupled through thermodynamic phase equilibria and saturation relationships (Gualda et al., 2012;Liu et al., 2005); (b) current laboratory studies only examine two-phase undersaturated melts containing crystals and melts and thus interpretations exclude the possibility of a three phase system containing crystals, melt, and a magmatic volatile phase (MVP) in saturated magmatic systems; and (c) MT interpretations generally require constraints on temperature and dissolved water content from petrological analyses of deposits from past eruptions and as such, inherently assume current and past magma reservoir conditions are similar.
To improve our understanding of the relationship between bulk resistivity and in-situ magma storage parameters in shallow rhyolitic systems, we use the thermodynamic software rhyolite-MELTS v. 1.1.0 (Gualda & Ghiorso, 2015;Gualda et al., 2012) to develop a parameterization for the melt volume fraction as a function of temperature, pressure and water content, and couple this to empirical relationships for bulk resistivity. Our method focuses on providing information on the dependencies of these parameters in rhyolitic magma reservoirs at or near closed-system thermodynamic equilibrium. In order to reduce complexities resulting from broader compositional variety, we restrict our interpretation to rhyolite magma reservoirs, although in principle, the method could be expanded to other compositions. Prior to this study, Samrock et al. (2021) have linked rhyolite-MELTS to bulk resistivity based on closed-system isobaric fractional crystallization of a parental basalt in a specific setting. Here, we focus more broadly on how melt fraction is coupled to P-T-H 2 O conditions in rhyolitic magma reservoirs, without a priori assumptions about the origin of the rhyolite magma or the temporal evolution of the reservoir.
Our method allows one to reduce the possible parameter space for melt fraction, temperature, and volatile content estimates at rhyolitic systems, and ensures these variables are internally consistent with thermodynamic melting relationships (i.e., higher temperatures correlate with higher melt fractions). Using the rhyolite-MELTS parameterization, we are also able to expand the bulk resistivity parameter space to consider water-saturated systems with a free MVP composed of aqueous fluids coexisting with the melt. Finally, the coupling provides a way to assess magma reservoir conditions independent of temperature and dissolved water content estimates from petrological analyses of previously erupted volcanic deposits, which could bias the interpretation of present-day in-situ magma storage conditions. We apply the method to three recent MT studies of silicic systems characterized by distinct resistivity anomalies and hence differing interpretations: (a) Mono Basin, California, USA where a moderately conductive (3-10 Ωm) anomaly was identified at 10 km depth (Peacock et al., 2015); (b) Newberry Volcano, Oregon, USA where a moderately resistive (25-50 Ωm) anomaly was identified at 3-4 km depth (Bowles-Martinez & Schultz, 2020); and (c) the Laguna del Maule Volcanic Field (LdMVF), central Chile where a highly conductive (0.3 Ωm) anomaly was identified at 4 km depth (Cordell et al., 2018).

Motivation for the Problem
Estimating magma storage conditions from bulk resistivity requires knowledge about the resistivity of the individual phases in the partial melt. Once the resistivity of each phase in the mixture is estimated, then mixing relations can be used to determine the volume fraction of the individual phases given an observed bulk resistivity. This section reviews the common approach that is widely used in the literature for estimating melt fraction and highlights its potential pitfalls and the need for additional constraints.
For under-saturated melts, the magmatic system contains two phases: the melt and the crystals (or solid minerals). Here, we use the laboratory-derived empirical relation of Guo et al. (2016) to calculate the rhyolite melt resistivity as a function of temperature and dissolved water content, where pressure has a second order effect (see Supporting Information S1). Figure 1a shows the parameter space for Guo et al. (2016) where higher temperatures and water contents lower the melt resistivity due to increased sodium ion mobility (Guo et al., 2016). Next, we apply the two-phase Modified Archie's Law (MAL) mixing relation from Glover et al. (2000) to estimate the bulk resistivity as a function of melt resistivity, crystal resistivity, melt fraction, and a user-defined parameter which reflects the connectivity of the melt phase (see Supporting Information S1). As shown in Supporting Information S1, the resistivity of the crystals has little impact in shallow silicic systems since it is at least two orders of magnitude more resistive than the melt phase, and the connectivity parameter can be held fixed between 1.0 and 1.5 for silicic partial melts (Gaillard & Marziano, 2005;Glover et al., 2000;Miller et al., 2015;Roberts & Tyburczy, 1999). Figure 1b shows the parameter space for MAL at a fixed crystal resistivity and connectivity parameter. We call this method "thermodynamically unconstrained" in that it does not take into account thermodynamic melting relationships.
Figures 1a and 1b describe a 4-D parameter space for the variables temperature, dissolved water content, melt fraction, and bulk resistivity. As such, if three of the four parameters are known (or can be estimated), then the fourth variable can be uniquely determined. Similarly, if one of the four parameters is known (or held fixed), then the other three can be plotted as a 3-D parameter space. Figure 1c shows the 3-D parameter space for water content, temperature, and melt fraction for a fixed 2 Ωm bulk resistivity. In other words, Figure 1c shows all the possible interpretations of magma storage conditions for a 2 Ωm anomaly using the empirical laboratory relation of Guo et al. (2016) and the MAL mixing relation. Note that with this approach, there is no thermodynamic constraint and thus lower melt fraction interpretations require hotter, wetter conditions while higher melt fraction interpretations require cooler, drier conditions. This is because higher temperatures and water contents result in lower melt resistivity values (Figure 1a), and thus less of the conductive melt phase is needed to explain a given bulk resistivity (Figure 1b). For example, as shown in Figure 1c, at 1000°C and 6 wt% H 2 O, the empirical relationships predict <20 vol% melt fraction, whereas intuition would suggest that rhyolite magma under these conditions is likely to be near-liquidus or super-liquidus. Similarly, at 700°C and 0 wt% H 2 O, the empirical relationships predict 100 vol% melt, yet an evolved magma under these conditions should be near-solidus, or even sub-solidus. This indicates that some parts of the thermodynamically unconstrained 4-D parameter space are not necessarily realistic for evolved silicic systems unless such systems are in a state of extreme disequilibrium. However, the empirical relationships in themselves do not provide any information about which parts of the parameter space might be unrealistic. Therefore, to identify interpretations that are potentially unrealistic or require conditions far from equilibrium, it is critical to further constrain the parameter space to a subspace which is consistent with thermodynamics. Thermodynamic modeling can be used to couple the four variables such that only thermodynamically consistent interpretations of magma storage conditions are considered. As shown in the next section, our new method constrains the 3-D parameter space to the white line in Figure 1c.  Guo et al. (2016). (b) Bulk resistivity using Modified Archie's Law from Glover et al. (2000) with 1,000 Ωm crystal resistivity and connectivity parameter of 1.5. (c) The thermodynamically unconstrained parameter space for a fixed bulk resistivity which combines panels (a and b). The white line on panel (c) shows the new thermodynamic constraint using rhyolite-MELTS which we develop in this study. All estimates made at 100 MPa.

Integrating Thermodynamic Constraints
Thermodynamic modeling software such as rhyolite-MELTS (Gualda & Ghiorso, 2015;Gualda et al., 2012) are commonly used to compute the thermodynamically stable equilibrium melt and mineral assemblage for a given bulk rock composition, temperature, pressure, and total water content. The outputs from rhyolite-MELTS include melt fraction and melt composition, the dissolved volatile content in the melt phase, and-in the case of saturated conditions-the amount of a water-dominated free MVP that would be present at equilibrium conditions. We assume an aqueous volatile phase (i.e., H 2 O) and ignore other volatiles (e.g., CO 2 ) which have a second-order effect on the estimated melt fractions for shallow rhyolitic systems.
To explore the parameter space over a range of temperatures, pressures, and total water contents, we constructed a parameterization of rhyolite-MELTS v. 1.1.0 using an average Andean rhyolite glass composition from the GEOROC database (Sarbas et al., 1999; see Supporting Information S1). Closed-system equilibrium crystallization was simulated from 1200°C to 660°C, with total system water contents between 0.5 wt% and 6 wt% (in 0.25 wt% increments) and pressures between 100 and 300 MPa (in 50 MPa increments). Crystal volume fraction ( ), melt volume fraction ( melt ), and MVP volume fraction ( ) were calculated and parameterized using the MATLAB curve fitting tool. The amount of dissolved water in the melt cannot exceed the water solubility limit of the melt, which is defined as a function of pressure and temperature from Liu et al. (2005). See Supporting Information S1 for more details about the parameterization.
The parameterization described above provides melt fraction, volatile fraction, and dissolved water content for an average Andean rhyolite as outputs, given inputs of temperature, pressure, and total system water content. Figure 2 shows the parameter space at 100 MPa. Melt fraction increases as a function of temperature and total system water content from cool colors (blue) to warm colors (yellow). The saturation curve on Figure 2 separates the parameter space into an under-saturated regime and a saturated regime. The saturation curve is near-vertical at 4 wt% total system water content related to first boiling (controlled by water solubility in the melt) and near-horizontal at 760°C related to second boiling (controlled by crystallization-induced degassing of the melt). The gray dashed contours in the under-saturated regime are the dissolved water content outputs from the MELTS parameterization and parallel the saturation curve, with decreasing dissolved water content at higher temperatures for sub-liquidus conditions. Note that the dissolved water content is measured in wt% of the melt, while all other measures of wt% and vol% are in reference to the total system. Once saturation occurs, the dissolved water content remains roughly constant, so no gray dashed contours appear in the saturated regime. However, exsolved volatiles form a third phase (the red dashed lines) which in this model co-exists with the other phases in equilibrium. Note that the liquidus line is here defined as anything >99 vol% melt, since the MELTS parameterization never reaches precisely 100 vol% (see Supporting Information S1). Also note that the dissolved water content contours in the super-liquidus under-saturated regime are vertical because the dissolved water content in the melt equals the total system water content.
The inputs and outputs of the parameterization include three of the four variables in the 4-D parameter space discussed in Section 2.1: temperature, dissolved water content, and melt fraction. The only remaining variable is bulk resistivity. Using the input temperature and output dissolved water content from the MELTS parameterization, the melt resistivity can be calculated using the empirical relation of Guo et al. (2016) at a fixed pressure. Then, the calculated melt resistivity and the output melt fraction from the MELTS parameterization can be used with the MAL mixing relation to compute the bulk resistivity, shown as solid black contours on Figure 2. The solid black contours are subparallel to the melt fraction, with the deviation due to a trade-off between temperature and dissolved water content on the melt resistivity. At higher dissolved water contents approaching the solubility curve, the bulk resistivity contours steepen whereas the melt fraction contours flatten. This is because the melt phase becomes substantially more conductive at the same time as the melt fraction becomes insensitive to the total system water (i.e., the melt is water-saturated), and thus less melt is required to maintain a fixed bulk resistivity.
Using this method, the bulk resistivity is no longer an unknown free parameter in a thermodynamically unconstrained 4-D parameter space (e.g., Figure 1) but is instead determined from the three MELTS-defined variables (temperature, dissolved water content, melt fraction) in a 3-D parameter space. This coupling provides bulk resistivity estimates which are consistent with undersaturated rhyolitic magma storage conditions based on equilibrium crystallization in MELTS. Extracting the temperature, dissolved water content, and melt fraction values along the 2 Ωm contour (white line in Figure 2) enables us to project the thermodynamic constraints onto the unconstrained parameter space (white line in Figure 1c). This restricts the available interpretations of a 2 Ωm anomaly to those that are thermodynamically consistent. Interpretations can still theoretically be made anywhere in the parameter space in Figure 1c, with the knowledge that as interpretations stray farther from the white line, they also require increasing levels of disequilibrium, which may be unrealistic at the reservoir scale.

Assumptions of the Method
We emphasize that nearly all estimates of melt fraction from MT (as well as seismic) data in the literature do not impose any thermodynamic constraints. As such, it is not clear whether previous interpretations of magmatic systems imply equilibrium or non-equilibrium conditions (e.g., Bowles-Martinez & Schultz, 2020;Cordell et al., 2018;Comeau et al., 2016;Hübert et al., 2018;Peacock et al., 2015;Samrock et al., 2018, etc.). Two exceptions are Blatter et al. (2022) and Samrock et al. (2021). The former was concerned with mantle melts in the asthenosphere and coupled a parameterization for peridotite melting with MT data to infer melt fraction and dissolved water content. The latter used isobaric fractional crystallization of parental basalt magma to infer compositional differences in subvolcanic magma systems from bulk resistivity. Here, we consider rhyolitic composition in an evolved system because we are primarily concerned with the conditions of long-lived silicic subvolcanic systems that are responsible for the largest and most hazardous eruptions.
We focus on determining the relationships of intensive parameters in rhyolitic reservoirs at or near in-situ closed-system equilibrium without an assumption of the origin of the rhyolite magma itself. Different processes can influence the generation of rhyolitic magmas, such as crystal fractionation, assimilation and partial melting of hot and deep crustal components (e.g., Moyen et al., 2021). In the present study, we consider the current state of an active shallow rhyolitic subvolcanic reservoir because we aim to image and identify zones of accumulation of rhyolitic magmas rather than the processes that led them to accumulate there. Note that we do not necessarily assume the entire trans-crustal magmatic system is rhyolitic, but only that the shallowest reservoir is so. As we assume a rhyolite composition, the limited compositional contrast between melt and crystal assemblage results in a negligible effect on melt fraction or dissolved water content when comparing closed and fractionated systems (see Supporting Information S1). Dynamic (i.e., non-equilibrium) processes of melt segregation (e.g., crystal settling) can change the geometry and distribution of the melt within a reservoir . For example, the melt may not be homogeneously distributed in the system but may instead collect in melt pockets or lenses. This highlights the distinction between "melt porosity" determined by the distribution and geometry of the melt within the reservoir and "melt fraction" defined thermodynamically. The MT method does not have sufficient resolution to image the specific distribution of melt within the reservoir and instead better estimates the total melt (i.e., melt fraction) within the reservoir (Cordell et al., 2022). While our method does not preclude the possibility of small zones of melt segregation in disequilibrium (below the resolution limits of the MT method), we therefore focus here on the large-scale state of the reservoir at near-equilibrium P-T-H 2 O conditions.

Water-Saturated Melts
The bulk resistivity in Figure 2 is only calculated in the under-saturated regime, and the interpretation in Figure 1 assumes under-saturated melts. This is because the existing empirical relations derived from laboratory experiments only deal with undersaturated melts (e.g., Guo et al., 2016Guo et al., , 2017Laumonier et al., 2015). As a result, nearly all MT interpretations at volcanic systems ignore the possibility of volatile saturated conditions which include a third volatile phase (some exceptions include Afanasyev et al., 2018;Samrock et al., 2021). However, the MELTS parameterization predicts the formation of an aqueous MVP if the water content in the system exceeds the solubility of water in the silicate melt (see Supporting Information S1). The MVP volume fraction is shown as dashed red lines on Figure 2. While the closed-system MELTS model assumes that all MVP remains in the system, we note that the efficiency of MVP migration and escape out of shallow magma bodies remains an open question. The buoyant migration of the MVP depends on the crystallinity of the magma and the volume fraction of MVP present Parmigiani et al., 2016). These studies show that within silicic magmas, bubbles can remain trapped within the melt for prolonged time; it is therefore not unrealistic to expect a finite MVP fraction present in most shallow water-saturated magma reservoirs.
Several complications arise when estimating the bulk resistivity of three-phase magmatic systems with an exsolved volatile phase since the salinity (i.e., resistivity) of the MVP is unknown and MELTS assumes a pure H 2 O (or H 2 O-CO 2 ) fluid. However, there will likely be interactions between the MVP and the surrounding crystals, melt, and/or wallrock which results in an aqueous solution with dissolved ions (Afanasyev et al., 2018). We can estimate end-member values for the MVP resistivity to place bounds on the possible bulk resistivity values in the saturated regime. If we assume the MVP is a saline supercritical aqueous fluid, then the MVP resistivity can be estimated as a function of temperature, pressure, and NaCl-equivalent concentration using the empirical relation derived from the laboratory experiments of Sinmyo and Keppler (2017) (see Supporting Information S1 for more information). We consider two end-member cases: a well-connected, high-salinity MVP representing a coherent brine layer, and a poorly connected, low-salinity MVP representing disconnected gas bubbles. The resistive, low salinity end-member considers 1 wt% NaCl-equivalent solution while the conductive, high salinity end-member uses 16 wt% NaCl-equivalent which represents an upper limit (Afanasyev et al., 2018). We use the estimated resistivity of the MVP (see Supporting Information S1) and the output MVP volume fraction-alongside the inputs/outputs from Section 2.2-to compute the bulk resistivity of the saturated melt using a three-phase MAL for each end-member scenario (Glover, 2010). See the Supporting Information S1 for more details about the three-phase mixing relation. As shown by the case studies in the next section, including the saturated regime when interpreting MT results can provide crucial information about magma storage conditions at silicic systems.

Case Studies
The method described in Section 2 can be applied to any system where an estimate of the bulk resistivity has been made for an inferred rhyolitic magma reservoir. The method could be expanded to other compositions (e.g., andesite, dacite) but would require a different MELTS parameterization. For case studies, we restrict our focus to silicic systems where previous estimates of melt fraction were made for shallow rhyolitic melt and where explicit details were provided by the authors about their assumptions. Based on these criteria, we chose to investigate Mono Basin, California, USA (Peacock et al., 2015), Newberry Volcano, Oregon, USA (Bowles-Martinez & Schultz, 2020), and the LdMVF, central Chile (Cordell et al., 2018). For each case study, we compare the thermodynamically unconstrained parameter space (i.e., Figure 1, termed "unconstrained"), to the "MELTS-constrained interpretation" (i.e., Figure 2).

Mono Basin, California, USA
Mono Basin is located in eastern California, USA and hosts the Mono-Inyo Craters volcanic chain, the youngest set of silicic eruptions in the Long Valley volcanic system. Petrologic data suggests that Mono Basin is underlain by a shallow rhyolitic magma system with primarily crystal-poor, highly evolved erupted products (Hildreth, 2004). MT data were previously collected throughout Mono Basin and inverted to produce a 3-D bulk resistivity model (Peacock et al., 2015). The model contained two primary features which were interpreted to be of magmatic origin (labeled C1 and C3; Peacock et al., 2015). Both are located at approximately 10 km depth with a bulk resistivity between 3 and 10 Ωm. Peacock et al. (2015) used the empirical relation of Gaillard (2004) to calculate rhyolite melt resistivity at a range of temperatures (750°C-850°C) and dissolved water contents (4-5 wt%), and applied the two-phase MAL to arrive at a melt fraction estimate between 10 and 20 vol%. A summary of temperature, water content, and other variables from Peacock et al. (2015) is shown in Table 1.

Thermodynamically Unconstrained Approach
We use the same input parameters as Peacock et al. (2015), but we replace the empirical relation for melt resistivity from Gaillard (2004) with that of Guo et al. (2016). Guo et al. (2016) is preferred since it does not require extrapolation past 3 wt% dissolved water content (see Supporting Information S1). Note that although Peacock et al. (2015) did not explicitly state the pressure or background resistivity, both variables have little impact on the melt fraction estimates (see Supporting Information S1). Here, we assume lithostatic pressure of 250 MPa for a crustal density of 2,500 kg/m 3 . The crystal resistivity is fixed at 1,000 Ωm which corresponds to leucogranite at 800°C (see Supporting Information S1; Hashim et al., 2013). Given the temperature and water content constraints (Table 1), melt resistivity is bounded between 0.4 and 1.0 Ωm (Figure 3a). This results in a melt fraction estimate of 7 vol% (red dot) to 17 vol% (blue dot) for a bulk resistivity of 6.5 Ωm (Figures 3b and 3c). This agrees well with the original melt fraction estimate of 10-20 vol% from Peacock et al. (2015). As discussed in Section 2.1, this thermodynamically unconstrained method results in hotter temperatures (red dots on Figure 3) corresponding to lower melt fractions and cooler temperatures (blue dots on Figure 3) corresponding to higher melt fractions. This is counter-intuitive to the expectations from thermodynamic melting relationships and better constraining this parameter space to be thermodynamically consistent is the motivation for our current study. If we expand the interpretation to include the range of bulk resistivity values from 3 to 10 Ωm, then the range of melt fraction estimates is much broader between 5 and 35 vol% (Figure 3b).

MELTS-Constrained Approach
The MELTS parameter space is shown in Figure 4 for a fixed pressure of 250 MPa with contours of bulk resistivity calculated in the under-saturated regime. The available interpretations of melt fraction, dissolved water Note. The "unconstrained" melt fraction estimates ignore thermodynamics and treat melt resistivity and melt fraction as independent. The new "MELTS-constrained" melt fraction estimates use MELTS to couple thermodynamics to bulk resistivity. MAL = Modified Archie's Law from Glover et al. (2000). content, temperature, and total system water content for under-saturated magma are constrained along the 6.5 Ωm contour. Without any additional constraint on temperature or dissolved water content, this results in a melt fraction estimate between 6 and 24 vol%. Unlike the thermodynamically unconstrained approach, the highest melt fraction estimate from the MELTS-constrained approach corresponds to the highest temperatures (∼850°C) and  lowest dissolved water contents (∼1 wt%), while the lowest melt fraction corresponds to the lowest temperatures (<750°C) and the highest water contents (>6 wt%), which are thermodynamically consistent estimates of magma storage conditions inferred from a 6.5 Ωm anomaly. Additional information that can be taken from Figure 4 is that the total system water content must be less than ∼0.5 wt% to be consistent with a 6.5 Ωm bulk resistivity. Note that this interpretation does not fix the dissolved water content at 4 to 5 wt% as was done in the unconstrained approach, but instead leaves dissolved water content as a free parameter. The white line on Figure 3c shows the thermodynamically constrained parameter space along the 6.5 Ωm contour in Figure 4 projected onto the unconstrained parameter space. This provides a new thermodynamic constraint to the parameter space and suggests that interpretations that are far away from the white line would be in a state of disequilibrium. If we expand the range of bulk resistivity values (3-10 Ωm), then the melt fraction estimates cover a broader range from 5 to 45 vol% with corresponding relationships to temperature and dissolved water content, and total system water content less than ∼1.0 wt% (Figure 4). In this case study, we do not need to invoke interpretations involving saturated melts to explain the bulk resistivity.

Newberry Volcano, Oregon, USA
Newberry Volcano is a broad shield volcano located in central Oregon, USA with a rhyolitic caldera at its summit, situated approximately 50 km east of the main Cascades volcanic arc. The volcano has had several intracaldera eruptions throughout the Holocene, all of which were rhyolitic (Macleod & Sherrod, 1988). MT data were collected inside Newberry caldera in the 1980s (Fitterman et al., 1988) and more recently on the west flank (Waibel et al., 2014). Bowles-Martinez and Schultz (2020) inverted the previous MT datasets to produce a 3-D bulk resistivity model. The model contained a moderately resistive feature 3-4 km beneath the caldera floor with a bulk resistivity of 25-50 Ωm. Bowles-Martinez and Schultz (2020) investigated various scenarios of melt resistivity, connectivity parameters, temperature and water content bounds using Guo et al. (2016) and the original Archie's Law (Archie, 1942) to explain the feature in terms of partial melt at 100 MPa. Their primary conclusion suggested a magma reservoir with 8-11 vol% partial melt at 850°C and 1.5 wt% dissolved water content. Their interpretation agrees with petrological estimates of temperature and water content from past eruptions, as well as melt fraction estimates of 8-12 vol% inferred from seismic results (Heath et al., 2015).

Thermodynamically Unconstrained Approach
For our interpretation to allow direct comparisons, we use the same input parameters as Bowles-Martinez and Schultz (2020) as shown in Table 1. Note that their study did not define a crystal resistivity because they used the original Archie's law (1,942) which assumes an infinitely resistive background resistivity. Here, we set the crystal resistivity to 1,000 Ωm as in Mono Basin. The choice of crystal resistivity has little impact on the result (see Supporting Information S1). Figure 5c shows the parameter space for melt fraction, temperature, and dissolved water content for a fixed bulk resistivity of 38 Ωm (i.e., midpoint of the 25-50 Ωm range) and fixed pressure of 100 MPa. At 850°C and 1.5 wt% dissolved water content, the estimated melt fraction is 9 vol%. If the range of bulk resistivity values encompasses the full range from 25 to 50 Ωm, then the range of melt fractions is 7-12 vol%, which is in good agreement with the original estimate of 8-11 vol% from Bowles-Martinez and Schultz (2020).

MELTS-Constrained Approach
According to the MELTS parameterization at 100 MPa shown in Figure 6, the magma storage conditions used in the thermodynamically unconstrained approach (i.e., 850°C and 1.5 wt% dissolved water content) require melt fractions >40 vol%, which in turn require a bulk resistivity of <4 Ωm. This is significantly more conductive than the 25-50 Ωm feature observed in the resistivity model of Bowles-Martinez and Schultz (2020). Indeed, the 38 Ωm and 50 Ωm contours do not appear in the under-saturated regime at all and the minimum resistivity bound (25 Ωm) narrowly appears in the under-saturated regime near 0 wt% total system water and <800°C ( Figure 6). This also means that it is not possible to plot the thermodynamically consistent part of the parameter space on Figure 5c using 38 Ωm. At face value, this suggests that the anomaly cannot be explained by an under-saturated magma reservoir if at (or near) thermodynamic equilibrium.
One possibility is that the magmatic system is over-saturated with a resistive, low salinity MVP. Figure 6 shows the bulk resistivity of the saturated regime, which was calculated using the three phase MAL for the resistive end-member MVP with 1 wt% NaCl (see Section 2.4). In this case, because the MVP is more resistive and disconnected than the melt phase, displacing melt with MVP increases the bulk resistivity. As a result, the bulk  . New thermodynamically constrained interpretation of Newberry using the MELTS-parameterization developed in this study with melt fraction as a function of temperature and total system water content at 100 MPa. The gray dashed lines are the dissolved water content contours, the thick black line is the saturation curve, the thick dashed line is the liquidus, and the red dashed lines are the magmatic volatile phase volume fraction contours. Thermodynamically constrained bulk resistivity is shown as thin black contours. The thick white line is the 38 Ωm bulk resistivity contour, the dashed white lines are the minimum (25 Ωm) and maximum (50 Ωm) contours estimated for the feature from Bowles-Martinez and Schultz (2020). Note that the white line only appears in the over-saturated regime. The white arrow shows potential cooling of the system from the magmatic conditions of past eruptions (at 850°C) to the observed bulk resistivity of the current system. resistivity contours from 25 to 50 Ωm appear in the over-saturated regime and extend over a range of temperatures <750°C, total system water contents <3 wt%, and melt fractions <12 vol%. Under such conditions, the MVP volume fraction is often comparable to (or greater than) the melt volume fraction component with a maximum MVP fraction of 25 vol% at 680°C for a 25 Ωm bulk resistivity.

Laguna del Maule Volcanic Field, Central Chile
The LdMVF is a silicic system in central Chile composed of a series of vents surrounding an alpine lake located approximately 30 km east of the main Andean volcanic arc. The LdMVF has erupted a variety of lava compositions prior to 25 ka, but since then has erupted nearly exclusively rhyolite (Andersen et al., 2018;Klug et al., 2020). Most notably, the LdMVF has been showing signs of renewed activity including rapid and ongoing ground deformation and seismicity (Cardona et al., 2018;Le Mével et al., 2020). MT data collected at the LdMVF in 2015 and 2016 were inverted to produce a 3-D bulk resistivity model (Cordell et al., 2018). The model contains several conductive features interpreted to be of magmatic origin including a feature labeled C3 near the center of ground deformation at a depth of 4 km (∼100 MPa) with a minimum bulk resistivity of approximately 0.3 Ωm (Cordell et al., 2018). Cordell et al. (2018) used the empirical relation of Guo et al. (2016) to estimate melt resistivity at 4-5 wt% dissolved water content and 760°C-1,000°C, based on petrological constraints from Andersen et al. (2018), then applied the MAL mixing relation to arrive at an estimated melt fraction >75 vol%. The authors concluded that such a high melt fraction may be unrealistic and suggested the presence of a brine lens. A summary of temperature, water content, and other variables from Cordell et al. (2018) is listed in Table 1.

Thermodynamically Unconstrained Approach
For our interpretation to allow direct comparisons, we use the same input parameters as Cordell et al. (2018) as shown in Table 1. Figure 7 shows the unconstrained parameter space for melt fraction, temperature, and dissolved water content for a fixed bulk resistivity of 0.3 Ωm and fixed pressure of 100 MPa. At 760°C-1000°C and 4-5 wt% dissolved water content, the melt fraction estimates range from 75 to 100 vol%, in agreement with Cordell et al. (2018). As discussed in Section 2.1, this thermodynamically unconstrained method results in hotter temperatures (red dots on Figure 7) corresponding to lower melt fractions and cooler temperatures (blue dots on Figure 7) corresponding to higher melt fractions. This is counter-intuitive to the expectations from thermodynamic melting relationships and further constraints are needed, which provides the motivation for our study.

MELTS-Constrained Approach
As can be seen in the MELTS parameter space in Figure 8, the 0.3 Ωm contour requires super-liquidus conditions in the under-saturated regime with temperatures >1000°C. This is well above the maximum pre-eruptive temperature estimates at the LdMVF (Andersen et al., 2018). This suggests that the interpretation of Cordell et al. (2018) with a melt fraction of 75 vol% using the thermodynamically unconstrained approach is not realistic under equilibrium conditions. At face value, these results would preclude an under-saturated magmatic origin for the 0.3 Ωm anomaly.
Like Newberry, one possibility is that the magmatic system is over-saturated with an MVP. In this case, we consider the conductive end-member MVP when calculating the bulk resistivity using the three phase MAL (see Section 2.4). Intuitively, a solute-rich, 16 wt% NaCl-equivalent aqueous solution should be highly conductive with bulk resistivity contours decreasing as the MVP volume fraction increases. However, the resistivity of the MVP increases as a function of temperature which means that, at high temperatures (e.g., >850°C), the super-critical aqueous MVP is more resistive than the melt phase. Thus, the bulk resistivity is relatively insensitive to the MVP at high temperatures ( Figure 8). As a result, the 0.3 Ωm anomaly still requires >1000°C in the saturated regime, significantly hotter than thermometry estimates for most rhyolites (Ellis et al., 2013). This suggests that the anomaly at LdMVF cannot be explained by saturated or under-saturated rhyolite at equilibrium conditions.

Discussion
The three case studies examined above represent very different MT features that span three orders of magnitude in bulk resistivity: a relatively resistive feature (25-50 Ωm) at Newberry volcano, a moderately conductive feature (3-10 Ωm) at Mono Basin, and a highly conductive feature (0.3 Ωm) at the LdMVF. These cases also illustrate several important benefits of the MELTS-constrained approach: excluding potentially unrealistic interpretations, expanding interpretations into the saturated regime, and providing independent insight into current reservoir conditions compared to estimates of temperature and dissolved water content from past erupted products.

Comparison of Thermodynamically Unconstrained Method and MELTS-Constrained Method
In all three cases, the thermodynamically unconstrained approach relied solely on interpretations of under-saturated partial melt and used petrological analyses of previously erupted volcanic products to restrict the range of temperature and dissolved water being considered. Using the MELTS-constrained bulk resistivity reveals that only specific interpretations are compatible with a system at (or near) equilibrium for an evolved silicic magmatic system. The best agreement between the two methods is at Mono Basin where 7-17 vol% melt was estimated using the unconstrained approach and 6-24 vol% was estimated using the MELTS-constrained approach, both of which are consistent with the original estimate from Peacock et al. (2015). However, simply comparing the final melt fraction values belies the differences between the approaches. First, the unconstrained approach uses the dissolved water content and temperature estimated from petrological analyses of erupted lavas, whereas the MELTS-constrained approach leaves both as free parameters; the 6.5 Ωm contour happens to cover the same temperature range (750°C-850°C) as previous petrological estimates. This need not be the case and is instead a powerful independent verification of the inferred magma reservoir conditions. A second important difference is that the MELTS-constrained approach leaves the dissolved water content as a free parameter, which is expected to be between 1 and 6 wt% based on the bulk resistivity. If the dissolved water content in erupted lavas is taken as a minimum bound (i.e., >4 wt%) in the MELTS-constrained approach, then the melt fraction range is narrowed to 6-15 vol% with temperatures restricted to <775°C. This agrees with lower temperature estimates at Mono Basin from Marcaida et al. (2014). Another important difference between the methods is that increasing temperature using the MELTS-constrained approach results in higher melt fractions and lower dissolved water contents: a system at 750°C has 6 vol% melt and >6 wt% dissolved water while a system at 850°C has 24 vol% melt and ∼1 wt% dissolved water. In contrast, the unconstrained approach suggests 7 vol% melt at 850°C and 5 wt% dissolved water (Figure 3c), but, as shown in Figure 4, such temperature and water content conditions would result in near-liquidus melts, and a bulk resistivity much less than 6.5 Ωm. Although the unconstrained interpretation is indeed possible, it requires that the bulk system is significantly out of equilibrium at the reservoir scale. This example at Mono Basin shows that, even when the melt fraction estimates derived from the thermodynamically unconstrained and MELTS-constrained approaches align well with one another, the MELTS-constrained approach provides additional information and avoids interpretations which are incompatible with thermodynamic constraints.
In contrast to Mono Basin, both Newberry and the LdMVF result in a mismatch between the thermodynamically unconstrained approach and the MELTS-constrained approach, which ultimately lead to divergent interpretations. In the case of Newberry, the relatively resistive 25-50 Ωm anomaly cannot be explained with 8-12 vol% undersaturated melt at 850°C and 1.5 wt% dissolved water when considering closed-system thermodynamic equilibrium, because an evolved magma under these conditions would have a melt fraction >40 vol% and a bulk resistivity near 2 Ωm ( Figure 6). The unconstrained interpretation is not impossible but would again require that the system is far from thermodynamic equilibrium at the reservoir scale. Knowing whether an MT interpretation is far from equilibrium is valuable in and of itself and raises the question of what the processes are that might generate such disequilibrium. Similarly, at LdMVF, the highly conductive 0.3 Ωm anomaly cannot be explained with 75 vol% undersaturated melt because a 0.3 Ωm bulk resistivity requires super-liquidus melt at temperatures >1000°C ( Figure 8). Both case studies indicate that alternative interpretations are required to explain these anomalies within a thermodynamically consistent framework.

Saturated Melts and Total System Water Content
The MELTS-constrained approach allows us to investigate whether water saturation can explain the bulk resistivity observed beneath silicic volcanoes. Interpretations also provide constraints on the total water content in the system, in addition to the dissolved water content of the melt phase. At Newberry, only low melt fraction (<12 vol%) saturated magmas at temperatures <750°C can explain the anomaly if and only if the MVP is poorly connected and/or low in salinity ( Figure 6). This could indicate, for example, that the magmatic system at Newberry resembles an MVP-rich fractured pluton in the final stages of solidification and degassing, driving circulation of geothermal fluids similar to those encountered during drilling campaigns (Fitterman et al., 1988). This follows interpretations from earlier MT studies (Fitterman et al., 1988), as well as seismic attenuation tomography (Zucca & Evans, 1992), which concluded that the caldera is underlain by a fractured pluton based on the velocity and seismic attenuation that was observed. More recent seismic work by Heath et al. (2015) and Beachly et al. (2012) used a denser seismic array and argued for a magmatic origin with 8-12 vol% melt due to a low velocity anomaly. However, seismic methods are insensitive to fluid compositions and cannot easily distinguish between MVP and melt. As a result, the low seismic velocity could be explained by a saturated system with a total melt + MVP of 8-12 vol%, in agreement with the MELTS-constrained MT interpretation. It is worth noting that the transition from an MVP-rich, near-solidus mush to a sub-solidus fractured pluton is gradual and thus defining Newberry as one or the other is ambiguous. The less ambiguous conclusion is that the temperature is likely less than former pre-eruptive temperatures of 850°C, and the system must be over-saturated with an MVP present in order to explain the 25-50 Ωm bulk resistivity.
At the LdMVF, neither the thermodynamically unconstrained approach nor the MELTS-constrained approach can explain the low resistivity anomaly using rhyolitic melt. Even when saturated partial melt with a solute-rich MVP is included in the MELTS-constrained interpretation, the 0.3 Ωm anomaly cannot be explained using magma reservoir conditions inferred from past eruptions. This is due to the fact that the MVP is slightly more resistive than the melt phase at temperatures >750°C, whereby adding more water to the water-saturated system does not significantly influence the bulk resistivity ( Figure 8). However, as the temperature decreases below 750°C, the conductivity of the MVP increases due to changes in charge mobility ( Figure 8). This results in non-linear behavior where the bulk resistivity in the saturated regime switches from decreasing as a function of temperature above 750°C, to increasing as a function of temperature below 750°C (Figure 8). Because of this non-linearity, the MT interpretation in the saturated regime is highly sensitive to the bulk resistivity. This does not change the interpretation of the 0.3 Ωm anomaly at the LdMVF, which cannot be explained with anything other than 100 vol% rhyolitic melt, even if the system is saturated. However, if the modeled feature was somewhat more resistive (e.g., 0.7 Ωm), then a wide range of saturated magma reservoir conditions can explain the anomaly within the temperature ranges suggested by petrological thermometry estimates, from melt-rich magmas with minimal MVP, to melt-poor magmas with MVP volume fraction >30% (Figure 8). Unlike the case studies from Newberry and Mono Basin, the work by Cordell et al. (2018) did not offer a range of bulk resistivity estimates and further work is needed to examine the uncertainty in the bulk resistivity of the LdMVF anomaly. A more likely possibility is that the feature at LdMVF contains an additional, more conductive phase which we did not consider. This could include cooler and hence highly conductive saline hydrothermal brines, or well-connected conductive minerals such as metal precipitates. Given the fact that the LdMVF is showing signs of significant unrest, another possibility is that the LdMVF is in a restless, dynamic state of thermodynamic disequilibrium. Hot, mafic recharge likely flushes volatiles and heat into the rhyolitic reservoir and has been linked to triggering of previous rhyolitic eruptions (Klug et al., 2020). Such mafic melts that are out of thermal and chemical equilibrium with cooler and/or more evolved reservoir areas at LdMVF could thus significantly reduce the conductivity. The MT model of Cordell et al. (2018) also images a deeper feature at ∼10 km depth which may be a less silicic reservoir that is providing material to the shallower rhyolitic reservoir. Further work is needed to elucidate these possibilities, which is beyond the scope of this paper.
It is worth noting that, at Mono Basin, there is no need to invoke a more complicated three-phase saturated system to explain the bulk resistivity. Although this result does not necessarily preclude the possibility of water saturation, undersaturation is consistent with the fact that the feature at Mono Basin is deeper than the features at Newberry and the LdMVF. As such, the magmatic system at Mono Basin is less likely to have reached the stalling level where saturation occurs Plank et al., 2013;Rasmussen et al., 2022).
These results also provide a constraint on the total system water content which is difficult to estimate through petrological methods that tend to focus on the dissolved water content in the melt. In the case of Mono Basin and Newberry, the total amount of water in the current system is estimated to be <3 wt% whereas all saturated interpretations at LdMVF suggest larger amounts of water in the system (>3 wt%). A rather dry system at Newberry could explain why geothermal fluids there are primarily of meteoric origin with limited magmatic signature (Sammel et al., 1988). Similarly, the Baños Campanario hot springs near the more hydrous LdMVF have a strong magmatic signature and some of the highest total dissolved solids observed at springs in central Chile, which is consistent with more magmatic water (Benavente et al., 2016). In the case that the LdMVF anomaly is not due solely to rhyolitic melt, the presence of saline brines and/or metal precipitates would still necessitate a fluid-rich magmatic system. The total system water content may thus provide insights into the larger-scale water buffering capacity of the crust and the amount of water available for exsolution in each system, with implications for example, volcanic hazard assessment and metallogenesis.

Past and Present Magma Reservoir Conditions
In all three case studies, hundreds (or thousands) of years have elapsed since the last eruption, which raises the question of whether past erupted products accurately reflect current reservoir conditions. The thermodynamically unconstrained approach (e.g., Figures 3, 5 and 7) uses temperature and dissolved water content estimates from past eruptions to infer current melt fraction, which may be biased if the magma system is no longer in the same state as it was prior to the last eruption. The MELTS-constrained approach provides an independent method for evaluating magmatic reservoir conditions since it uses thermodynamic modeling and bulk resistivity to constrain temperature, dissolved water content, and melt fraction independent of estimates of temperature and dissolved water content from petrological analyses of erupted products that record the reservoir conditions at the time of past eruptions.
At Mono Basin, the 6.5 Ωm contour falls within the temperature ranges and dissolved water content ranges inferred from past erupted products, suggesting that the past and current magma reservoir conditions are similar. In contrast to Mono Basin, Newberry volcano shows a mismatch between the pre-eruptive temperature and dissolved water contents of erupted lavas at Newberry compared to those inferred from the observed bulk resistivity. The MELTS-coupled bulk resistivity estimates indicate that past reservoir conditions at 850°C and 1.5 wt% dissolved water content require >40 vol% under-saturated melt with a bulk resistivity of ∼2 Ωm. The observed bulk resistivity of 25-50 Ωm suggests that either the system is very far from equilibrium in the undersaturated regime (e.g., <10 vol% melt at 850°C) or that the system is near to equilibrium in the saturated regime and has cooled significantly since the last eruption to <750°C (white arrow on Figure 6). The latter interpretation suggests that total system water content has remained roughly constant as a result of degassed MVP remaining trapped in the pore space of the crystallizing pluton. A trapped MVP could further explain why geothermal fluids at Newberry have very little magmatic signature since it implies a limited hydrologic connection between the cooling, degassing pluton and shallower meteoric hydrothermal reservoirs.
The fact that Mono Basin shows similar magma reservoir conditions to past erupted products while Newberry and the LdMVF do not, could be related to the depth of the features. Mono Basin is the deepest of the three features at 10 km depth, perhaps reflecting a long-term storage regime in the mid-crust with relatively stable temperatures and much higher water solubility limits. This could result in the melt remaining water-undersaturated, lending the system a greater capacity to retain its water rather than lose it during escape of exsolved MVP. In contrast, the features at Newberry and LdMVF are located at shallower depths (<4 km), where it is more difficult for the magmatic system to maintain a long-lived stable magma reservoir over the interval of thousands of years since the last eruption without significant heat flux and where loss of water through MVP exsolution is more likely.

Conclusions
This study couples bulk resistivity to thermodynamic constraints from a rhyolite-MELTS parameterization to improve estimates of magma storage conditions derived from MT data. Resulting interpretations from bulk resistivity models are thermodynamically consistent and avoid unrealistic tradeoffs (e.g., high temperatures being used to infer low melt fractions). The MELTS-constrained bulk resistivity not only constrains melt fraction, tempera ture and dissolved water content, but also places bounds on saturation conditions and total system water content in the present-day system. The fact that the MELTS-constrained bulk resistivity is independent of estimates of temperature and dissolved water content from petrological analyses of erupted products at a particular volcano allows for a comparison between past magma reservoir conditions as recorded in eruptive products, and present magma reservoir conditions as imaged with MT. Finally, the MELTS-constrained approach expands interpretation to water-saturated conditions where an MVP is present, although more work is needed to better constrain the volume fraction, salinity and fluid-rock interactions of the MVP, the network topology of the three-phase mixture, and the resulting bulk resistivity.
We demonstrated the benefits of the MELTS-constrained bulk resistivity in interpretation via three case studies with three different anomalies. Mono Basin can be explained by under-saturated partial melts with melt fractions between 6 and 15 vol%, dissolved water contents >4 wt%, and temperatures <775°C. These conditions are in good agreement with estimates of temperature and dissolved water content from petrological analysis of past eruptive products, suggesting relatively stable magma reservoir conditions. In contrast, the relatively resistive feature at Newberry Volcano cannot be explained by under-saturated partial melt at temperatures >750°C, and may instead reflect a cooler (<750°C), water-saturated melt-poor mush that has cooled significantly since the last eruptions.
Interpretations at 850°C and <10 vol% melt are very far from equilibrium which raises the question of whether such interpretations are realistic at the reservoir scale. Finally, at the LdMVF, the MELTS-constrained bulk resistivity suggests that neither saturated nor under-saturated rhyolitic magma reservoir conditions can adequately explain the MT anomaly even when considering a high-salinity, well-connected volatile phase. Instead, our analysis suggests that the anomaly at the LdMVF requires additional conductive phases (e.g., sub-critical solute-rich hydrothermal brines, metal precipitates, or hot mafic material). These examples illustrate the wide range of cases where MT bulk resistivity observations coupled with thermodynamic constraints provide crucial additional information about current magma reservoir conditions.