DensityX : A program for calculating the densities of magmatic liquids up to 1,627 °C and 30 kbar

Here we present a standalone program, DensityX , to calculate the densities of dry and hydrous silicate melts given pressures, temperatures, and major oxide compositions in wt% in the 10-component system SiO 2 TiO 2 Al 2 O 3 Fe 2 O 3 FeO MgO CaO Na 2 O K 2 O H 2 O. Here we use DensityX to analyze over 3,000 melt inclusions over a wide compositional range to visualize the distribution of natural silicate liquid densities in the Earth’s crust. The program is open-source, written in Python (as a library hosted in PyPi or as a downloadable Python script), and can be accessed and run via an online interface through a web browser at https://densityx.herokuapp.com or by downloading and running the code from a github repository. A compan-ion Excel spreadsheet can also be used to run density calculations identical to those in the Python script but only for one sample at a time. In another example application, we demonstrate how DensityX can be used to constrain density-driven convective cycling in the phonolitic lava lake of Erebus volcano, Antarctica.


Introduction
The density of a silicate melt affects a myriad of physical and magmatic processes, including magma mixing [Grove and Baker 1983;Jull and Kelemen 2001], melt migration dynamics [Stolper et al. 1981;Hack and Thompson 2011], and crustal storage [Sparks and Huppert 1984;Walker 1989;Chaussard and Amelung 2014]. Together, these processes lead to the diversity of magma types observed on Earth. Melt density is also an important parameter in planetary differentiation, where density-driven stratification segregates a planetary body into a crust, mantle, and core [Agee and Walker 1988;Ohtani 1985;Circone and Agee 1996]. A significant amount of effort has been put into modeling silicate melt densities, which relies upon determination of thermodynamic parameters such as partial molar volumes of major oxide from experiments conducted at a range of pressure (P ) and temperature (T ) conditions and expressions describing changes in the volume of these components as a function of P and T (e.g., related to thermal expansion and compressibility of ions). Much of this work was pioneered by R.A. Lange, V.C. Kress, and coworkers who produced a series of papers describing the densities of the major components of silicate liquids [Lange and Carmichael 1987], thermodynamic relations needed to calculate silicate melt density at elevated P and T [Lange and Carmichael 1990;Kress and Carmichael 1991;Lange 1997], and the important addition of thermodynamic data related to the density of H 2 O in silicate melts [Ochs III and Lange 1999;Richet et al. 2000]. This work resulted in an equation of state (EOS) for silicate melts in the 10-component system SiO 2 TiO 2 Al 2 O 3 Fe 2 O 3 FeO MgO CaO-Na 2 O-K 2 O H 2 O [Lange and Carmichael 1990]. The thermodynamic values underpinning these models (e.g. liquid partial molar volumes and melt compressibilities) were recalibrated by Ghiorso [2004], who increased the size of the database used for earlier calibrations and provided a more rigorous treatment of the effect of the oxidation state of Fe on the densities of melts. In the Ghiorso-Kress approach, Fe-bearing systems are speciated into three (rather than two) components: FeO, Fe 2 O 3 , and FeO 1.3 . This new treatment and expanded database resulted in an improved model recovery of density measurements in Fe-bearing systems but requires knowledge of the melt oxygen fugacity.
Over the last four decades, research in this area has produced a large knowledge base of thermodynamic relations as well as the formulae needed to interpret geologic data. The resulting formulae for density calculations are typically presented as thermodynamic equations in the literature or integrated into larger thermodynamic modeling programs such as MELTS (first published as Ghiorso and Sack [1995]), Perple_X [Connolly 2005], Theriak_D [Duesterhoeft and de Capitani 2013], and the model of Fluegel [2007], which rely on the above mentioned published datasets. However, Iacovino & Till, 2019 no standalone program for calculating the density of a silicate melt as a function of P and T exists, making the calculation cumbersome. Automated calculations within currently available thermodynamic modeling programs are useful but require a working knowledge of each program's infrastructure, are time consuming to perform on very large datasets, and the inclusion of the density calculation within the programs means that it is not readily extensible by the end-user.

DensityX
Here we present a standalone program, called Den-sityX, to calculate the densities of hydrous silicate melts given pressures, temperatures, and major oxide compositions in wt% in the 10-component system SiO 2 TiO 2 Al 2 O 3 Fe 2 O 3 FeO MgO CaO Na 2 O− K 2 O H 2 O, using thermodynamic data for these components from Lange and Carmichael [1987], Lange [1997], Kress and Carmichael [1991], and Ochs III and Lange [1999].

Program DensityX
The program DensityX is based on previously published thermodynamic properties of silicate liquids given in Lange and Carmichael [1987], Lange and Carmichael [1990], and Ochs III and Lange [1999] and equations to describe the pressure-and temperaturedependent partial molar volumes of oxide components given in Lange and Carmichael [1990]. The density of a silicate melt depends upon the individual densities of the oxide components of the melt, which in turn depend on their partial molar volumes. Partial molar volumes of the components will depend upon several factors, including overall melt composition, degree of melt polymerization, and changes to coordination of cations with pressure. Here we use an empirical description of melt density, which assumes ideal mixing between components [Lange and Carmichael 1990]: where ρ liq is the density of the silicate liquid, X i is the mole fraction and M i is the molecular weight of oxide component i in the melt, and V liq is the volume of the silicate liquid at a given P and T defined as the sum of the partial molar volumes of each oxide component at P and T and normalized to the concentration of each component, as: where V i is the calculated partial molar volume at given P and T of the oxide component i in the melt.
To calculate melt density at P and T , DensityX first calculates the individual molar volumes of each component oxide at P and T . We use a simplified model equation that incorporates the effects of pressure (dV i /dP ) and temperature (dV i /dT ) on volume, which is sufficient to describe the volume of natural liquids up to 30 kbar [Lange and Carmichael 1990]: where V i (T ref , 1 bar) is the partial molar volume of component i in the liquid at reference temperature T ref and a pressure of 1 bar and dV i /dT and dV i /dP are the temperature and pressure derivatives of the partial molar volume of component i. The pressure dependence of the partial molar volume will be stronger with pressure, and so the second pressure derivative d 2 V i /dP 2 would be required to calculate V i at P > 30 kbar [Lange and Carmichael 1990].
DensityX employs the EOS of Lange and Carmichael [1990], since the updated Ghiorso-Kress EOS requires knowledge of the system's oxygen fugacity, which is often not known (or not reported) for both natural and experimental samples. A comparison of density values calculated with the DensityX and Ghiorso-Kress databases are given in Appendix 1. Values used in the DensityX model are calibrated for silicate melt compositions ranging from: SiO 2 = 37-75 mol%; TiO 2 ≤ 4 mol%; Al 2 O 3 ≤ 27 mol%; Fe 2 O 3 ≤ 15 mol%; MgO ≤ 38 mol%; CaO ≤ 43 mol%; Na 2 O ≤ 33 mol%; K 2 O ≤ 29 mol%; and H 2 O < 19 mol% [Lange 1997;Ochs III and Lange 1999]. The model is valid in the range 1-30 kbar [Lange and Carmichael 1990] and from the glass transition temperature for any given melt up to 1627°C (1900 K; [Lange 1997]). Below the glass transition temperature, rearrangements in the glass structure become kinetically impeded, and so a correction must be applied when calculating the density of a glass below this temperature. This requires knowledge of the compositionally dependent values for the glass transition temperature and coefficient of thermal expansion, and so this correction is not applied in DensityX. Calculations at temperatures below the glass transition will therefore have higher error than reported here.

Description of the program package
DensityX is capable of performing batch calculations on large sample sets (thousands of discrete samples) automatically. The calculation of each sample density is performed iteratively (not in parallel), but the execution is fast. The example input file used here (and available with the model code and online interface) contains over 3,000 individual samples and was processed in <10s during testing (3 GHz Intel Core i7 processor). The code is open source and simple enough to allow the end-user to adjust or update the program with minimal working knowledge of the computer language Python. DensityX is also packaged Table 1 -Partial molar volumes and pressure and temperature derivatives for each oxide component. For batch calculations, an input file with the oxide compositions, P , and T of each sample must be prepared by the user. The file test-data.xlsx is included with the program (and can be downloaded via the web interface) and serves as a template for an input file. Although this file can have any name, it must follow the structure of test-data.xlsx, with columns named as follows (with composition in wt%, pressure in bars, and T in°C): Sample_ID, SiO 2 , TiO 2 , Al 2 O 3 , Fe 2 O 3 , FeO, MgO, CaO, K 2 O, H 2 O, P , T . The order of the columns does not matter, as the program locates data columns in the input file based on the column name given in the first row of the file. If using the densityx Python library, a Pandas DataFrame with the column names as described above serves as input data. On the web interface, or if using the DensityX.py script, the user will be prompted to browse for a .xlsx file.
Upon choosing the input file (or upon passing a Pandas DataFrame), the following operations are done on the data to calculate density. First, necessary constants are defined for all oxide components (Table 1): molecular weights; partial molar volumes at the reference P and T of 1 bar and 1773 K, V i [Lange and Carmichael 1987;Lange 1997;Ochs III and Lange 1999]; thermal expansion coefficients expressed as the temperature derivative of V i , dV i /dT [op. cit.]; and compressibility coefficients expressed as the pressure derivative of V i , dV i /dP [Kress and Carmichael 1991;Ochs III and Lange 1999].
The program then solves Equation 1 for each sample in the input file in two steps: by solving first for the numerator (Equation 2) and then for the denominator (the liquid volume, Equation 3). The data is then written to an excel spreadsheet and saved as [Input File Name]_output.xlsx to the folder where the input file exists. When using the densityx library, a Pandas DataFrame containing density values is returned (see /lib/README.md available in the github repository for more information on how to interact with the densityx library).

Densities of natural silicate liquids
As magma bodies ascend, cool, crystalize, and degas, the density of the liquid component will change due to a number of processes including: decreasing pressure during crustal ascent; decreasing temperature during magma storage and cooling; increasing temperature due to magma mixing or mafic recharge events; changing major element composition during differentiation; and changing H 2 O concentration during magma ascent and/or volatile exsolution.
To illustrate the compositional dependence of the density of natural silicate liquids, geochemical data for over 3,000 crystal-hosted melt inclusions from volcanic settings around the world ranging in SiO 2 concentration from ∼40-80 wt% were compiled from the GEOROC database (Figure 1; data from GEOROC http://georoc.mpch-mainz.gwdg.de/georoc/). These melt inclusions represent residual liquids isolated from crystallizing magma bodies in the Earth's crust. Because the P and T of melt inclusion entrapment is rarely known or reported, all samples were run at P = 5,000 bar and T = 1,100°C, representing magma storage conditions in the mid-crust. DensityX was used to calculate the densities of all samples in this dataset in two model runs: one using anhydrous compositions and one including reported dissolved H 2 O concentrations. The data used here are given test-data.xlsx (Supplementary file 2). Anhydrous melt inclusion compositions range in density from ca. 2.4-3.0 g/cm 3 ( Figure 1B). Of the anhydrous components, SiO 2 and mafic components MgO+FeO correlate most strongly with melt density (R 2 = 0.91, 0.81 respectively). This is likely due in part to the relative abundance of these oxides in silicate melts and to the propensity of SiO 2 concentration to be anticorrelative with MgO+FeO in natural silicate melts. When reported H 2 O concentrations are included in the total Fe as FeO model run, the density range increases significantly to ca. 2.0-3.0 g/cm 3 , with a much higher variability in density for any given magma composition ( Figure 1A). H 2 O exhibits such a strong control on density due to its molecular abundance in natural liquids (H 2 O commonly accounts for up to ∼25 mol% in silicate melts) and its large partial molar volume (V H 2 O = 22.9 ± 0.6 cm 3 /mol; [Ochs III and Lange 1999]) relative to other major oxide components.
During magma ascent and differentiation, pressure and temperature changes will play a strong role on the evolution of the liquid melt density. DensityX was run for a mid-ocean ridge basalt [Dixon et al. 1995] and a high-silica rhyolite [Wallace et al. 1999] with varying H 2 O concentrations from ca. 0-10% over a range of T at P = 5,000 bar and over a range of P at T = 1,000°C (Figure 2).
The dependence of liquid density on temperature increases with increasing H 2 O concentration and/or decreasing SiO 2 concentration. Thus, T has the strongest effect on hydrous basalts and the smallest effect on dry rhyolitic melts. In basalt, an addition of 1 wt% H 2 O has an effect on liquid density equivalent to an increase in T of ∼300°C or a decrease in P of ∼4,000 bar. In rhyolite, that effect is equivalent to an increase in T of ∼500°C or a decrease in P of ∼2,500 bar. This effect is particularly important when calculating densi-  Dixon et al. [1995] using an Fe 2 O 3 /FeO ratio from experiment #7 (op cit.). Anhydrous rhyolite composition: Bishop Tuff sample BT87-2 [Wallace et al. 1999] using an Fe 2 O 3 /FeO ratio calculated at the quartz-fayalite-magnetite buffer using the model of Kress and Carmichael [1991]. ties from melt inclusions, which may lose significant H 2 O due to diffusive re-equilibration during storage in a region shallower than the original inclusion entrapment depth. At constant T , the effect of changing P on melt density is stronger in rhyolite than in basalt, due to the compressibility of SiO 2 . The effect is slightly more pronounced at higher H 2 O concentrations. At 1,000°C and 10 wt% H 2 O, basalt densities range from 2.2-2.9 g/cm 3 (dρ/dP = 2.4 × 10 −5 g/cm 3 -bar) while rhyolite densities range from 1.9-2.7 g/cm 3 (dρ/dP = 2.8 × 10 −5 g/cm 3 -bar).

Applications to high-P melts
In the years since the development of the Lange and Carmichael [1990] equation of state and the publication of the partial molar properties used in Densi-tyX (Table 1), significant progress has been made on the properties of melts, including density, over a wide range of compositions and H 2 O concentrations [San-loup 2016;Terasaki and Nishida 2018]. Many of these recent works employed density measurements for various melts up to pressures well exceeding the limits of DensityX ( 100 kbar). It has been illustrated that third order Equations of State become necessary to accurately predict densities of melts at very high pressures due to changes in the coordination of various cations [Sanloup et al. 2013a;Sanloup et al. 2013b]. Figure 3 shows a comparison of density estimates by DensityX to those reported from recent in situ density measurements up to 30 kbar. DensityX reproduces the majority of reported densities to within a 2σ error of 5% (indicated by dashed lines). Thus, the database used for DensityX remains valid up to at least 30 kbar in light of recent advancements regarding very high pressure melts.
6 Example application to density-driven convection in the Mount Erebus lava lake Density-driven convection in volcanic systems is sometimes invoked to explain the growth of large, zoned crystals. This process is exemplified in systems with continuously convecting lava lakes, such as Erebus (Antarctica) and Nyiragongo (Democratic Republic of Congo), which host zoned megacrysts up to 10 cm in size. At Erebus, oscillatory chemical zones in anorthoclase megacrysts erupted from the phonolite lava lake are thought to record magmatic convection cycles reaching as deep as ∼7 km (2,000 bar) with little to no temperature change during convective cycling [Moussallam et al. 2015]. Evidence from crystal zones and observations of cyclicity in lava lake eruptive behavior (e.g. gas bubble bursts) and gas compositions provide good evidence that convection is driven by cycles of degassing at the surface followed by volatile replenishment in a shallow magma reservoir [Oppenheimer et al. 2009;Ilanko et al. 2015].
Using DensityX, we calculated the pressure-density paths of volatile-bearing (0.2 wt% H 2 O) and degassed Erebus phonolite lava between 0-2,500 bar at a constant T of 1000°C (Figure 4). At magmatic temperatures, Erebus anorthoclase crystals have a density of 2.524 g/cm 3 [Moussallam et al. 2015], and as such remain negatively buoyant within anhydrous phono-lite magma at pressures up to about 2,000 bar (illustrated in Figure 4). This is consistent with a scenario in which anorthoclase crystals continuously settle through a volatile-poor magma conduit [e.g. Molina et al. 2012] until they reach ∼2,000 bar. The periodic infiltration of rehydrated magmas, which are positively buoyant relative to the surrounding degassed phonolite, could provide upward momentum sufficient to carry anorthoclase crystals back up toward the lava lake, recorded as cyclic chemical zoning in individual crystals. The exsolution of volatiles and generation of bubbles within upwelling magma will further decrease magma density, propelling it to the surface faster. Density curves in Figure 4 suggest that magma rehydration must occur at a minimum pressure of 2,000 bar, the point of neutral buoyancy for anorthoclase crystals within degassed phonolite melt. The actual rehydration pressure may be greater but would not be recorded by anorthoclase crystals, which are too buoyant to reach such depths.
DensityX does not consider the degassing of CO 2 , which may have a significant effect on carbon-rich melts. CO 2 has a larger partial molar volume than H 2 O at the same temperature [Sakamaki et al. 2011] and so the addition of CO 2 to the melt would decrease the density of volatile-bearing Erebus phonolite (Figure 4). However, despite Erebus's carbon-rich nature, lava lake phonolites are stored very shallowly and so are significantly degassed, containing only about 0.03 wt% CO 2 Figure 4: Density-driven convection within the shallow plumbing of Erebus volcano is likely the result of cycles of near-surface degassing followed by deep rehydration of phonolite melts at a minimum pressure of ∼2,000 bar (horizontal red dotted line). Low-density volatilebearing melts (blue curve) rise to the surface where they degas and sink back into a shallow magma storage zone (red curve) where they are replenished with volatiles. This cycling is recorded in cyclic chemical zoning within anorthoclase megacrysts (black dashed line) that settle to a maximum depth of ∼2,000 bar and are carried upward during pulsatory ascent of hydrated magma batches. [Iacovino 2015]. Thus, it is unlikely the addition of CO 2 to our model would change our interpretations for Erebus drastically.
Other parameters can affect the efficiency of crystal settling and buoyancy in a magma. For example, the high viscosity of Erebus phonolite suggests that there will be a competition between the crystal settling time and the upwelling time. This may result in retrograde cycles, where sinking crystals are pushed upward before having the chance to completely settle to their point of neutral buoyancy. Such complex cyclicity is observed in zoning patterns of natural anorthoclase megacrysts [Moussallam et al. 2015]. High phonolite viscosity may also result in crystals being dragged below their point of neutral buoyancy at 2,000 bars, but it is unlikely that crystals could sink much below this limit before buoyant forces overcome melt viscosity.

Conclusions
DensityX is a powerful and extensible tool for calculating densities of silicate melts at pressures up to 30 kbar and temperatures up to 1,627°C. The program is written to be easily incorporated into existing code written in Python via the densityx PyPi library. In addition, a standalone Python script capable of performing density calculations on 1000s of samples using a provided Microsoft Excel input file and a standalone Microsoft Excel spreadsheet capable of calculating the density of a single sample allows DensityX to be utilized by users of almost any computer skill level. We have demonstrated that DensityX reproduces a wide range of experimentally measured silicate melt densities to within a 2σ error of 5%. DensityX was used to calculate the densitypressure relationships of dry and hydrated phonolite melt from Erebus volcano (Antarctica). A comparison with the known density of anorthoclase megacrysts, which record dehydration-rehydration cycles in the upper Erebus conduit, allowed us to constrain the rehydration depth to a minimum pressure of 2,000 bar. tal. They also thank Mark Ghiorso, Erik Duesterhoeft, Chrystele Sanloup, Charles LeLosq, and an anonymous reviewer for insightful comments on the manuscript and Magdalena Oryaëlle Chevrel for editorial support. K.I. was supported by National Science Foundation Frontiers in Earth System Dynamics (FESD) grant no. 1338810.

Author contributions
KI conceived of the project, wrote the DensityX code and online interface, and authored the manuscript. CBT assisted with the example applications to crustal magma density and in writing the manuscript.

Data availability
The DensityX program demonstrated in this manuscript is open-source and MIT licensed. The full version of this code, which supports multiple sample input, is written in Python and is available in the github repository https://github.com/kaylai/DensityX. Den-