An equation of state for high pressure-temperature liquids (RTpress) with application to MgSiO3 melt

The thermophysical properties of molten silicates at extreme conditions are crucial for understanding the early evolution of Earth and other massive rocky planets, which is marked by giant impacts capable of producing deep magma oceans. Cooling and crystallization of molten mantles are sensitive to the densities and adiabatic profiles of high-pressure molten silicates, demanding accurate Equation of State (EOS) models to predict the early evolution of planetary interiors. Unfortunately, EOS modeling for liquids at high P-T conditions is difficult due to constantly evolving liquid structure. The Rosenfeld-Tarazona (RT) model provides a physically sensible and accurate description of liquids but is limited to constant volume heating paths (Rosenfeld and Tarazona, 1998). We develop a high P-T EOS for liquids, called RTpress, which uses a generalized Rosenfeld-Tarazona model as a thermal perturbation to isothermal and adiabatic reference compression curves. This approach provides a thermodynamically consistent EOS which remains accurate over a large P-T range and depends on a limited number of physically meaningful parameters that can be determined empirically from either simulated or experimental datasets. As a first application, we model MgSiO3 melt representing a simplified rocky mantle chemistry. The model parameters are fitted to the MD simulations of both Spera et al. (2011) and de Koker and Stixrude (2009), recovering pressures, volumes, and internal energies to within 0.6 GPa, 0.1 Å3, and 6 meV per atom on average (for the higher resolution data set), as well as accurately predicting ∗Corresponding author. Email addresses: aswolf@umich.edu (Aaron S. Wolf ), daniel.bower@csh.unibe.ch (Dan J. Bower) Preprint submitted to Physics of the Earth and Planetary Interiors liquid densities and temperatures from shock-wave experiments on MgSiO3 glass. The fitted EOS is used to determine adiabatic thermal profiles, revealing the approximate thermal structure of a fully molten magma ocean like that of the early Earth. These adiabats, which are in strong agreement for both fitted models, are shown to be sufficiently steep to produce either a center-outwards or bottom-up style of crystallization, depending on the curvature of the mantle melting curve (liquidus), with a high-curvature model yielding crystallization at depths of roughly 80 GPa (Stixrude et al., 2009) whereas a nearly-flat experimentally determined liquidus implies bottom-up crystallization (Andrault et al., 2011).


Abstract
The thermophysical properties of molten silicates at extreme conditions are crucial for understanding the early evolution of Earth and other massive rocky planets, which is marked by giant impacts capable of producing deep magma oceans. Cooling and crystallization of molten mantles are sensitive to the densities and adiabatic profiles of high-pressure molten silicates, demanding accurate Equation of State (EOS) models to predict the early evolution of planetary interiors. Unfortunately, EOS modeling for liquids at high P-T conditions is difficult due to constantly evolving liquid structure. The Rosenfeld-Tarazona (RT) model provides a physically sensible and accurate description of liquids but is limited to constant volume heating paths (Rosenfeld and Tarazona, 1998). We develop a high P-T EOS for liquids, called RTpress, which uses a generalized Rosenfeld-Tarazona model as a thermal perturbation to isothermal and adiabatic reference compression curves. This approach provides a thermodynamically consistent EOS which remains accurate over a large P-T range and depends on a limited number of physically meaningful parameters that can be determined empirically from either simulated or experimental datasets. As a first application, we model MgSiO 3 melt representing a simplified rocky mantle chemistry. The model parameters are fitted to the MD simulations of both Spera et al. (2011) and, recovering pressures, volumes, and internal energies to within 0.6 GPa, 0.1Å 3 , and 6 meV per atom on average (for the higher resolution data set), as well as accurately predicting

Introduction
Accurately describing the physical and thermodynamic properties of liquids is a challenging task, due to their disordered nature. Unlike solids, where atomic bonding is constrained by lattice symmetries, atoms in liquids possess a partially randomized structure that evolves readily in response to changes in pressure, temperature, and composition. In general, liquids only retain time-averaged short-range order, since individual atomic bonds are constantly breaking and reforming in response to bulk stresses and localized atomic diffusion (Stebbins, 1988). Due to these complexities, precise equation of state (EOS) modeling for liquids is quite challenging. The highly symmetric lattice structures of crystalline solids are well-suited to approximate representations of both compressive and thermal properties, like the Mie-Grüneisen-Debye model. Such simple and powerful EOS models are generally not available for liquids, however.
Despite this limitation, understanding the thermodynamic properties of liquids remains an important goal in many physical science applications. In the geological sciences, melts play a primary role in planetary evolution, acting as efficient agents for transporting both heat and chemical species. In particular, MgSiO 3 melt represents the dominant chemical component of the mantles of rocky planets, and thus its properties strongly influence volcanism and the early stages of terrestrial planet formation. High-pressure melts may also play an important role today at the base of Earth's mantle near the coremantle-boundary, potentially explaining the observed low seismic velocities in ultra low velocity zones (Nomura et al., 2011;Williams and Garnero, 1996).
During the late stages of planetary accretion, planetary embryo collisions were a commonplace occurrence, with most Earth-mass planets experiencing one to a few giant impacts (Quintana et al., 2016;Chambers, 2001;Agnor et al., 1999). These giant impacts delivered heat to great depth, with Earth's largest moon-forming impact potentially melting the mantle all the way to the core-mantle boundary (Nakajima and Stevenson, 2015). The thermophysical properties of these high-pressure molten silicates strongly influence the convection and cooling evolution of the resulting deep magma oceans. Properties such as density, heat capacity, thermal expansion, and compressibility are all more sensitive to pressure and temperature for liquids than their corresponding solid phases, due to their continuous structural evolution (Wolf et al., 2015;Jing and Karato, 2011). The evolutionary path of magma oceans rests on the coupled interactions between thermodynamics and fluid dynamics (e.g. Solomatov, 2000). Of particular interest are the relative buoyancy of melts-compared to their mineral counterparts-which controls whether newly-formed crystals sink or float (Stolper et al., 1981), influencing the possible creation of a long-lived basal magma ocean at the core-mantle boundary (Labrosse et al., 2007). Additionally, the thermal gradient of the melt relative to the liquidus (or melting curve) dictates the depth of crystallization ). These and other first-order properties of high pressure silicate melts have yet to be fully explored and have the potential to fundamentally alter the early evolutionary path of rocky planets (e.g. Andrault et al., 2017).
To enable effective study of high-pressure phenomena in silicate melts, we develop a new liquid EOS for high P-T conditions, the RTpress model. This semi-empirical EOS uses a thermal perturbation approach based on the original Rosenfeld-Tarazona (RT) equation (Rosenfeld and Tarazona, 1998), which has been shown to effectively describe a wide array of liquid types (Ingebrigtsen et al., 2013). Unfortunately, the RT model has limited utility for high pressure applications, since its standard form applies only to thermal profiles at constant-volume. To extend the RT model to high pressures, we introduce a reference isothermal compression curve combined with a reference adiabatic thermal profile. Thermal deviations from these reference profiles are described using a generalized form of the Rosenfeld-Tarazona equation. This piecewise modeling approach provides a powerful framework for modeling the thermodynamic properties of liquids at high pressures and temperatures, utilizing many familiar EOS parameters while accurately capturing the compressive and thermal properties unique to liquids. As a first application of the RTpress model, we determine the equation of state of high-pressure MgSiO 3 melt, representing a simplified chemical model for terrestrial planetary mantles. The model is fit to molecular dynamics simulations from Spera et al. (2011) and, and both are nicely consistent with the experimental shock-wave measurements of Mosenfelder et al. (2009) and broadly consistent with one another at high pressures. We finally explore the thermal structure of a deep magma ocean and the possible implications for the depth of crystallization.

Adiabatic Melt Compression
The distinct compressive properties of liquids, as compared to their solid counterparts, arise from their freedom to structurally evolve during compression. The thermodynamic consequences of liquid structure are readily seen in the compression behavior of the thermodynamic Grüneisen parameter, which is a unitless quantity that describes the relation between temperature, T , and volume, V , along an isentrope, or path of constant entropy, S: where the Grüneisen parameter, γ, can be estimated indirectly from measurements of specific heat capacity, C V or C P , density ρ, thermal expansion α, and bulk modulus, K T or K S ; (see e.g., Boehler and Kennedy, 1977;Boehler et al., 1979;Boehler and Ramakrishnan, 1980). Understanding the behavior of the Grüneisen parameter is critical when studying convecting systems, like magma oceans, due to its control over the thermal profile: where P is the pressure. Convecting systems adopt nearly adiabatic (constant entropy) profiles as they carry heat from depth up toward the surface, with the Grüneisen parameter setting the scale factor for the steepness of this thermal profile. Accurate simulations of magma ocean cooling thus rely on correctly modeling the unusual compression properties of silicate melts. For crystalline solids, this thermodynamic quantity relates directly to the vibrational Grüneisen parameter, and the two are approximately equal for many materials (see, Kieffer, 1982;Williams et al., 1993). When combined with the quasi-harmonic approximation for low-amplitude vibrations, this equivalence leads to the common simplifying assumption that the Grüneisen parameter depends only on volume, with negligible temperature dependence. It has been universally observed for solids that compression causes the value of the Grüneisen parameter to decrease, in response to anharmonicity (Anderson, 1974). The popular Mie-Grüneisen equation of state is a direct consequence of these solid approximations, and can reasonably model solid behavior up to moderate temperatures below the melting point, where vibrations remain harmonic (Oganov and Dorogokupets, 2004;Wu and Wentzcovitch, 2009).
In contrast to solids, liquids universally show paradoxical Grüneisenbehavior, as first noted by Shapiro (1969, 1970) who calculated γ values from thermodynamic property tables for both liquid water and mercury over an 80 K temperature range for compression ratios of up to 20%. This study reveals non-solid-like behavior of γ for liquids in terms of both increasing values with compression and considerable temperature dependence. Since then, these same compression and thermal behaviors have been shown to hold for a wide array of liquid types and compositions. These include simplified model systems (Lennard-Jones fluid), cryogenic noble gases (Argon), organic solvents (pentane, ethanol, methanol), metallic liquids (mercury), and a broad array of silicate and oxide melts (Boehler and Ramakrishnan, 1980;Brown et al., 1987;Amoros et al., 1988;Hess et al., 1998;de Koker and Stixrude, 2009;Spera et al., 2011;Asimow, 2012).

RTpress EOS Model
To address the challenge of accurately representing the thermodynamics of liquid compression, we develop the RTpress model, an extension of the Rosenfeld Tarazona EOS to high P-T conditions. Important features of this semi-empirical EOS are its applicability to a wide range of liquids and its physically meaningful parameterization, which enables easy analysis and interpretation. It is designed with a relatively simple analytic form that is readily usable, and the developed analysis code is freely available as open source Python software (see github.com/aswolf/xmeos). Rosenfeld and Tarazona (1998) proposed a simple mathematical form for a thermal equation of state generally appropriate to liquids. They derive this expression using the variational perturbation method for hard spheres (Mansoori and Canfield, 1969). Despite its idealized foundation, this theoretical framework has been shown to apply well to a broad class of liquids (including simplified model liquids, molten salts, hydrocarbon fluids, and liquid water; Ingebrigtsen et al., 2013). Molecular dynamics (MD) simulations of silicate melts, including SiO 2 (Saika- Voivod et al., 2000), Mg 2 SiO 4 (Martin et al., 2009), MgSiO 3 (Spera et al., 2011), and CaAl 2 Si 2 O 8 , have also been shown to be well-represented by the RT-EOS form.

Rosenfeld-Tarazona EOS
The original RT model represents the total potential energy of liquids as a simple shifted power-law: where a and b are free parameters specific to a particular constant-volume (isochoric) heating profile, and must be determined empirically by fitting data from molecular dynamics simulations. An additional kinetic energy term enables the RT model to accurately describe the properties of a wide variety of simulated liquids. Initially, its limitation to constant volume profiles posed a barrier to wide-spread use. To provide a complete thermodynamic description, an EOS must describe the full free-energy surface as a function of both temperature and volume (or pressure). The compression properties of the RT model constants, a and b, cannot be directly constrained by basic physics, except in the overly-simplified case of a Lennard-Jones Fluid (Ingebrigtsen et al., 2013). A practical approach was introduced by Saika-Voivod et al. (2000)-and later extended by Ghiorso and Spera (2010)-to expand the RT model to cover a range of pressures by assuming smooth behavior for the RT coefficients. Saika-Voivod et al. (2000) consider a grid of isochoric paths and assume that the RT coefficients progressively evolve as a function of volume, ensuring a smooth and differentiable free energy surface. High-order polynomials describe the volume-dependence of both coefficients, a(V ) and b(V ), introducing ∼10-14 arbitrary free parameters that must be determined by fitting simulated isochores (in addition to other EOS parameters). Though goodness-of-fit does not directly translate to predictive accuracy for interpolation or extrapolation, it is remarkable that this approach enables recovery of the training data at the 0.1% to 0.01% level in some cases over much of the simulated P-T range (e.g., Saika-Voivod et al., 2000;Ghiorso and Spera, 2010). The unavoidable downside of this method, however, is that it relies on an unphysical parameterization of the compression-dependence, introducing a large number of free parameters whose estimation requires abundant high-precision simulation data, introducing overfitting issues while excluding experimental applications.

RTpress Thermal Perturbation
By building on previous efforts, we extend the Rosenfeld-Tarazona EOS to high pressure by employing a thermal perturbation approach. The pri- There are three independent model components, highlighted in green boxes: 1. isothermal compression curve, 2. reference adiabat, and 3. generalized Rosenfeld-Tarazona thermal model. Descriptions of each of these components are given in Secs. 3.3 and 3.4. The corresponding parameters for each model component are: 1. the compression parameters {V 0 , K 0 , K 0 }, 2. the Grüneisen model parameters {γ 0 , γ 0 }, and 3. the RT parameters {b 0 , ..., b n } that describe the compression dependence of the heat capacity. Note that the RT offset coefficients (a's from Eq. 3) are now implicitly defined by the adiabatic reference curve, and thus are no longer required as fitting parameters (see details in Sec. 3.4). To constrain these model parameters, we use data of three different types, shown in gray boxes: (I) isothermal compression data, (II) high P-T data providing the thermal pressure derivative, and (III) data to constrain the heat capacity.
Heat capacity data for the liquid can come from multiple possible sources, as indicated in Fig 1 (panels IIIa-IIIc). Direct heat capacity information, either at constant volume or pressure, can be obtained either from experimental calorimetry or molecular dynamics simulations (IIIa). Data constraining an adiabatic temperature path can be measured in laser-driven ramp-wave compression experiments (Kraus et al., 2016;Luo et al., 2004a) (IIIb). The temperature effects on the internal energy can be obtained directly from molecular dynamics simulations (IIIc). This model can thus be adapted to take advantage of a range of possible calibration data sources, including both experimental and theoretical techniques, unlike all previous RT-based models. For the example of MgSiO 3 shown later in this paper, we use method IIIc to constrain the heat capacity of the liquid, by simultaneously fitting pressure and internal energy as functions of volume and temperature based on MD simulation data.
The thermophysical properties of the RTpress model are obtained from derivatives of the free energy surface: The corresponding internal energy and entropy surfaces are given by: where each term is rooted in the component pieces of the RTpress model depicted in Fig. 1: 1. F 0T (V ) is the work performed along the isothermal compression curve, 2. T 0S is the thermal path along the reference adiabat (with entropy S 0 ), and 3. ∆E and ∆S are the perturbations in internal energy and entropy, respectively, determined by the Rosenfeld-Tarazona thermal model. The total free energy surface is thus readily calculated using the analytic expressions from each component of RTpress.

Reference Profiles
The reference compression curve can by described using any standard isothermal compression model (e.g., Vinet or Birch-Murnaghan Poirier, 2000). We generally prefer the Vinet model (Vinet et al., 1989) for its superior extrapolation performance to very high pressures (Cohen et al., 2000). However, the fourth-order Birch-Murnaghan is a good alternative for cases where the Vinet is insufficient to capture rapid changes in the bulk modulus and additional degrees of freedom are required. The pressure and energy expressions for the Vinet model are: with, where the model parameters are the familiar zero-pressure volume, V 0 , the isothermal bulk modulus K 0 , and its pressure-derivative K 0 . Since the compression curve represents a constant-temperature path (T = T 0 ), the work integral gives the Helmholtz free energy F 0T (V ) = − V V 0 P 0T (V )dV + E 0 , where the constant of integration, E 0 , defines the energy level at the reference temperature and 0 GPa.
RTpress also relies on a direct description of the reference adiabatic temperature profile T 0S (V ), where S = S 0 . The reference isentropic temperature path can be described in terms of an integral of the thermodynamic Grüneisen parameter (Eq. 1) along the reference adiabat. This reference Grüneisen curve ensures physically meaningful behavior given a sensible representation of γ 0S (V ). Numerous analytic forms are available to describe γ 0S , and we adopt the vibrational finite-strain model of Stixrude and Lithgow-Bertelloni (2005) based on its superior performance at high compression (for more details, see App. A). The isentropic temperature profile is given by: is the finite volumetric strain and the Taylor expansion coefficients, a 1 = 6γ 0 and a 2 = −12γ 0 +36γ 2 0 −18γ 0 , are expressed in terms of the value and derivative of the Grüneisen parameter at zero GPa, γ 0 and γ 0 = V 0 (dγ/dV ) 0 . For more information, see Sec. 5.2 and App. A.

Generalized Rosenfeld-Tarazona EOS
The thermal perturbation for RTpress is given by a generalized form of the original RT model (Eq. 3), adjusted to improve modeling flexibility: defined in terms of the unitless thermal deviation from the reference temperature, f T : Similar to the approach of de Koker and Stixrude (2009), the thermal deviation and its temperature derivative, f T , are defined in terms of the powerlaw exponent m, which is potentially allowed to deviate from its theoretically expected value of m = 3/5 (determined by Rosenfeld and Tarazona, 1998). Following previous efforts, the prefactor for the power-law term, b, is represented as a function of volume. By introducing the thermal deviation which equals zero when T = T 0 , we ensure that the offset term-given by a in the original RT model-now refers directly to the reference isothermal internal energy curve a(V ) = E pot (V, T 0 ). We have thus removed the additional 5-7 free parameters that are required by the polynomial RT model Saika-Voivod et al. (e.g. 2000); Ghiorso and Spera (e.g. 2010) to effectively constrain the temperature-independent energy offset term, a(V ).
We adopt a modified polynomial representation of the thermal coefficients b(V ) (Saika-Voivod et al., 2000;Ghiorso and Spera, 2010): where (V /V 0 − 1) is the unitless volumetric deviation, and b n are the fitted polynomial parameters. This representation simplifies interpretation of the model parameters: b 0 represents the thermal coefficient at the reference volume b(V 0 ) = b 0 , and successive parameters signify the linear, quadratic, and high order dependences of b(V ) on percentage changes in volume. Additionally, all b n parameters retain units of energy, and should generally decay in magnitude as the polynomial order increases.
To enable the thermal perturbation calculation, we determine internal energy differences from the reference isotherm: and entropy differences from the reference adiabat: T (T 0S ) define the thermal perturbation and the thermal perturbation derivative relative to the appropriate reference paths. For more details, including analytic expressions for important quantities, see App. B.

Thermodynamic Properties of the RTpress Model
Expressions for all thermodynamic quantities can be obtained from the appropriate derivatives of the free energy surface (Eq. 4). The pressure, for instance, is given by the volume derivative, P (V, T ) = − dF dV T : where P (V, T 0 ) is the isothermal reference contribution from Eq. 6, and ∆P E and ∆P S are the energetic and entropic contributions from the free energy: where b (V ) is the volume derivative of the RT model coefficient (see Eq. B.2), and C V,0S (V ) and γ 0S (V ) are the heat capacity and Grüneisen parameter evaluated along the reference adiabat (see Eqs. A.3 and B.4). Similar expressions can also be obtained for other quantities like the Grüneisen parameter (discussed in detail in Sec. 5.2), as well as the bulk modulus and thermal expansion.

Application to MgSiO 3 Melt
As a first application of this new liquid EOS model, we will examine the properties of MgSiO 3 melt at high pressure. Geologically, this is an important liquid composition that represents a simple mantle model for rocky planets like the Earth. Molten silicates play a crucial role in the formation of rocky planets, which is marked by giant impacts, leading to the formation of deep magma oceans. Even billions of years after formation, it is still possible that melts may play a geologically significant role near the core-mantle boundary in the form of a long-lived basal magma ocean or as localized pockets of melt (Labrosse et al., 2007;Nomura et al., 2011). Thus, MgSiO 3 represents a good initial melt model to understand geologic evolutionary processes operating throughout rocky planetary mantles.
To determine the EOS of MgSiO 3 melt, we use the high resolution data produced by the classical molecular dynamics simulations of Spera et al. 2011 (yielding model S11). These calculations employ the empirical pair potentials of Oganov et al. (2000), based on the structure and thermodynamic properties of pure endmember bridgmanite. The liquid MD simulations were performed for a set of approximate isotherms from ∼2500 K to 5000 K that span 0 to 150 GPa. Spera et al. (2011) showed that the Oganov potential agrees reasonably well with first principles MD simulations of Stixrude and Karki (2005); de Koker and Stixrude (2009), as well as other classical MD simulations using alternate potentials (Lacks et al., 2007). Importantly, Spera et al. (2011) published a complete set of simulated output data including calculations of pressure and energy over a range of volumes and temperatures. To demonstrate the versatility of our new EOS form, we also fit the first principles MD simulations reported in de Koker and Stixrude 2009 (model dK09), obtaining similar performance and comparable properties at high pressures and temperatures.
To fit these data, we have developed an open-source Python package called xmeos: Xtal-Melt EOS. This software is freely available on the author's github website (github.com/aswolf/xmeos), and includes a sizable library of equation of state functions relevant to solid and liquid high-pressure phases. Additionally, it provides many routines to enable evaluation and fitting of thermodynamic properties for both mineral and melt phases.
In order to determine the EOS parameters for MgSiO 3 melt, we employ weighted least squares (or χ 2 ) minimization, fitting simulated energies and pressures as functions of volume and temperature. To combine these two data types in a simultaneous fitting procedure, we assume that the errors in the data are proportional to the standard deviation of each quantity within the training data set: where f is the uncertainty scale-factor which controls the size of the model parameter uncertainties, and is determined from the residual scatter to the best-fit model (assuming that the model adequately captures the overall trends in the data). The goal of this joint fitting procedure is to obtain a model that simultaneously behaves well in terms of energy, pressure, and volume; due to the high compressibility of liquids, however, small deviations in pressure can introduce large volume (or density) errors near 0 GPa. This causes major issues when using the model to assess near-surface conditions of a magma ocean. To address this shortcoming, we introduce an effective pressure uncertainty term,σ P , that equally weights deviations in both pressure and volume, ensuring that the model maintains accuracy over the full magma ocean pressure range (see App. D for details). The fitting cost function is given by the familiar χ 2 formulation, with one summation term for pressure and one for energy: Prior to fitting, we must fix the reference temperature T 0 to a sensible value, in order to define the isothermal reference path and the potential temperature of the reference adiabat. We select T 0 = 3000 K, since it represents a well-sampled isotherm within the training dataset as well as corresponding to an adiabat that lies reasonably close to the liquidus over the mantle pressure range. Since we are most interested in the behavior of the melt near the liquidus, this is the best strategy to ensure a physically sensible model while enabling easy evaluation of the parameters. Beyond following these general guidelines, the model is fairly insensitive to the chosen reference temperature, with all fitted parameters jointly adjusting their values to bestmatch the data for the chosen reference state. Uncertainties in the model parameters are obtained using standard error propagation techniques. For more details of the fitting procedure see App. D. The best-fit RTpress model (S11) is shown in Fig. 2, where each model isotherm is color-coded and matched to the corresponding data from Spera et al. (2011). As shown, the model captures the MD simulation data to high accuracy, with average residual uncertainties of only 0.6 GPa, 0.1Å 3 /atom, and 6 meV/atom in pressure, volume, and energy, with a corresponding coefficient of determination of R 2 ≈ 0.9998. This goodness of fit is comparable to the polynomial RT model fit by Spera et al. (2011), despite having five fewer free parameters, thus providing a more physically and statistically robust description of the thermodynamic properties. The model parameters and individual uncertainties are given in Tbl. 1. A similar procedure was carried out for the FPMD data of de Koker and , and the resulting fitted parameters are also given in Tbl. 1 (see App. D for details). Covariations between the fitted parameters are expressed by the parameter correlation matrix (Fig. D.2), which shows that significant covariations (magnitude ≥ 0.9) are limited to parameter-pairs within each of the three parts of the model (for example V 0 -K 0 , K 0 -K 0 , and b 3 -b 4 ). Correlations across model components are much smaller, all falling within the range of ±0.7. This pattern indicates that each component of the model  (70) is semi-independent of the others, lending further credibility to the model formulation.

Comparison with other Liquid EOS Models
As discussed in Sec. 3.4, the Rosenfeld-Tarazona Model has been previously extended to high pressure using a simple polynomial representation of the RT coefficients a(V ) and b(V ) (e.g. Saika-Voivod et al., 2000). A similar approach is to directly represent the Helmholtz free energy surface with a multi-dimensional Taylor expansion in volumetric strain and thermal deviation (developed in Lithgow-Bertelloni, 2005, 2011). As with the polynomial RT model, thermodynamic consistency is guaranteed, since it fully represents the free energy surface and all thermodynamic properties are obtained from appropriate derivatives.
The connection between these two methods is strong, since the thermal parameterization of the Taylor expansion approach is equivalent to the RT power-law dependence when truncated to first-order in temperature. The major difference is that higher order terms can be retained in the Taylor expansion method, though they tend to significantly amplify overfitting issues.
De  showed that for a range of silicate melts along the MgO-SiO 2 binary, a first-order expansion in temperature-equivalent to a generalized RT model-was sufficient to reasonably represent the MD simulation data (pure silica was the exception, which required second order thermal terms). Unlike in all of the classical MD studies employing the polynomial RT model, de Koker and  found that the fitted power-law exponent deviated significantly from its theoretical value of 3/5 (Rosenfeld and Tarazona, 1998). To test this possibility we also explored a model in which the power-law exponent was included as a free parameter, though we found that the goodness of fit to the data showed negligible improvement (R 2 increased by only 0.00002), with pressure and volume residuals decreasing by only 0.01 GPa and 0.1Å 3 and the energy residual actually worsening by 0.1 meV. When repeated for the deKoker2009 training data, we found that while the best-fit value increased to about m fit ≈ 1.2, it still suffered from similar insignificant improvement in the goodness of fit (with the R 2 value improving by only 0.0004). Given that the data were not found to strongly constrain the power law exponent, we favor fixing it to its theoretical value of m = 0.6 (as reported in Tbl 1).
Though these previous EOS methods have been reasonably successful, both have issues with overfitting and are generally restricted to molecular dynamics data only. The RT polynomial approach relies on abundant and extremely precise isochoric potential energy heating curves. Such datasets, in practice, are restricted to well-sampled MD simulations, where the conditions of the liquid can be precisely chosen, and the calculation is efficient enough to allow very long equilibration runs. Furthermore, both methods encourage the use of a large set of fitting parameters that are not directly physically meaningful, and are thus lead to significant overfitting. Both de Koker and Saika-Voivod et al. (2000) discuss how EOS models of silicate melts suffer from unphysical behavior when extrapolated outside the fitted data region. These melts appear to decompose into two immiscible liquids of equal composition outside the simulated P-T range. This behavior can easily arise from unphysical wiggles in the free energy surface introduced by excessively high order polynomials not adequately constrained by the data. Furthermore, it is also difficult to impose any constraints on the fitted parameters for these models to ensure sensible behavior because the physically relevant quantities are complex functions of multiple parameters. Instead, we must rely upon having enough accurate data coverage to constrain the model coefficients, which then indirectly determine the thermodynamic properties. This approach imposes necessarily large correlations between the fitted coefficients, introducing greater difficulties in the inter-comparison of different EOS studies and the evaluation of model uncertainties.
The RTpress model thus provides a highly valuable and new contribution to this family of liquid EOS models. By stitching together multiple modelcomponents already familiar to the mineral physics community, we are able to focus our attention on only physically meaningful model parameters. This eases model assessment and enables the use of priors (when necessary) to help constrain parameters to reasonable values when modeling limited or lower accuracy datasets. Furthermore, the number of parameters are considerably smaller, requiring 3-10 fewer parameters, depending on order of truncation required for competing models. The effect of this parameter reduction can be plainly seen in the modest uncertainties and correlation values, especially in the limited correlations displayed between parameters belonging to differing model-components, emphasizing how the data can provide independent constraints on each component of the RTpress model. Most importantly, this model sets itself apart from other liquid EOS forms in its applicability to a wide variety of data types, including not only molecular dynamics simulations but high P -T experimental studies as well.
In order to directly assess the quality of our model, we compare its predictions of the shock Hugoniot for MgSiO 3 melt to the experimental gas gun measurements reported in Mosenfelder et al. (2009), including collected results from numerous other authors (Akins, 2003;Akins et al., 2004;Luo et al., 2004b;Simakov and Trunin, 1973). The shock Hugoniot represents the collection of thermodynamic states that are possible behind a shock wave (arising from the conservation of mass, momentum, and energy). The final thermodynamic state depends sensitively on the initial state of the sample, and thus we show experimental data for both MgSiO 3 glass and crystalline enstatite starting samples, shown in black and gray, respectively, in Figure 3. We compare these experimental measurements to the theoretically predicted curves from our new EOS, showing results for both the S11 and dK09 parameter sets in solid and dashed lines, respectively (see App. E for more details). As clearly demonstrated in the figure, the region of pressuretemperature space relevant to the mantles of large rocky planets like the Earth (roughly 50-150 GPa and 3000-6000 K), are directly sampled by the glass Hugoniot. For the initially glassy samples, we predict both densities and temperatures that are in good agreement with experiments (though the FPMD-based model dK09 is less able to accurately capture the thermal properties of the liquid, under-predicting shock temperatures by ∼1000 K). The enstatite Hugoniot data lie at much higher pressures and lower temperatures, indicated by the gray region, which is less directly relevant to conditions during planet formation and evolution. For these data, both EOS where both models generally agree with one another except for a noticeable reduction in shock temperatures for dK09. The RTpress model provides a good match to densities and temperatures for glass starting material, which occupies the relevant region of P − T space for deep magma oceans. The Hugoniot data for enstatite starting material is also shown in gray, which is less well represented by either fitted model, but is less relevant to terrestrial magma oceans and these data may suffer from incomplete melting. models show much weaker agreement with experiments (with dK09 slightly outperforming S11). As pointed out by Thomas and Asimow (2013), some of these enstatite data likely suffer from incomplete phase transformation and may also depart from thermodynamic equilibrium. We therefore rely upon the glass Hugoniot data-where our model shows remarkable predictive accuracy-both for its geologically relevant P-T conditions and its freedom from difficulties with incomplete melting.

Adiabatic Compression of MgSiO 3 Melt
By directly modeling the reference Grüneisen compression profile, the RTpress EOS is nicely able to capture both the pressure-and temperaturedependence of adiabatic liquid compression. Using the standard thermodynamic derivative approach, we obtain the expression for the Grüneisen parameter of the RTpress model: where C V (V, T ) is the constant-volume heat capacity given by Eq. B.4. When evaluated along the reference adiabat, T = T 0S , the γ expression reverts to the reference γ-profile, as expected.
To evaluate the effectiveness of this temperature-dependent Grüneisen parameter formulation, we pay close attention to the choice of the reference gamma-model, which must be capable of reflecting the paradoxical compression behavior characteristic of liquids. In particular, it must allow for an initially increasing γ value with compression. This general behavior of liquids was discussed in detail in Wolf et al. (2015) and a semi-quantitative mechanism was proposed, relating to structural evolution in the liquid, which affects the representative vibrational frequency of the material. Upon compression, the liquid re-orders itself, adopting more compact and solid-like structures. As a consequence, this structural evolution is self-limiting: eventually, the liquid attains a near-close-packed structure where further rearrangement is no longer possible. Thus, all liquids must eventually give way to solid-like compression behavior, where gamma decreases with compression at high pressures.
The self-limiting character of the liquid-like Grüneisen compression is not typically well-represented in many common γ-models. The most popular power-law γ-model, γ(V ) = γ 0 (V /V 0 ) q , assumes that the fractional trend in the Grüneisen parameter is constant with compression, and is thus incapable of representing the required γ-trend turnover; the same is true of many other common forms. The finite strain γ-model, initially developed by Stixrude and Lithgow-Bertelloni (2005) as a more physically-motivated model for solids, however, guarantees the presence of a turnover when applied to liquids, where the Grüneisen parameter always returns to solid-like behavior after sufficient compression. This finite-strain model rests upon the variation of a characteristic vibrational frequency with compression, which was shown by Wolf et al. (2015) to naturally capture the entropic properties of MgO melt while automatically incorporating the physically required selflimiting property. To test out the affect of our selection of the finite-strain γ-model, we compare our preferred model to one using the standard powerlaw formulation, which yields a considerably worse fit to the internal energy curves with energy residuals that are ∼3.5 times larger (15 meV/atom). For these reasons, the finite strain model stands out as the best choice for accurately representing the γ-compression of liquids. The adiabats from the two RTpress models are in strong agreement shown in solid (S11) and dashed (dK09) lines color-coded by their potential temperature. For comparison, the mantle melting curve (liquidus) is represented by black dashed  and dotted (Andrault et al., 2011) lines. The lower panel compares how competing models differ from the RTpress model in terms of deviation percentage from the 2500 K adiabatic temperature gradientwith the shock-based model of Mosenfelder et al. (2009) as dotted and FPMD-based model of Stixrude et al. (2009) as crossed lines. The RTpress models nicely capture the average thermal gradient from the shock wave data (near-zero average deviation), while the Stixrude2009 model is systematically steeper everywhere by >30% (despite having the same training data as our dK09 model). The inset panel shows the compression evolution of the Grüneisen parameter, which controls the steepness of the adiabats, highlighting its liquid-characteristic temperature-dependence.

Pressure [GPa]
The temperature-dependent Grüneisen evolution of MgSiO 3 melt along adiabatic profiles is shown in the inset of Fig. 4. The compression behavior of γ is dominated by an initial very sharp rise (at all temperatures) until roughly ∼15 GPa, after which its value begins to plateau. This sudden change in slope primarily reflects the rapid rise in the bulk modulus of liquids as their structure becomes more solid-like with compression, as modeled by characteristically small values of K 0 and large values of K 0 (see Tbl. 1). Temperature acts as a second-order modification of the γ-trend, with high temperatures reducing the γ value by a relatively constant amount over much of the mantle pressure range. Roughly speaking, the Grüneisen parameter is seen to asymptotically approach a limiting value just above 1 as the liquids become increasingly solid-like, with temperature lowering its value by about 0.1 per 1000 K increase.

Crystallization of a Magma Ocean
Our EOS model for MgSiO 3 melt provides an approximate prediction for the thermal structure of an early Earth magma ocean. Given this reasonable simplified mantle chemistry, we can integrate the adiabatic pressure gradient (Eq. 2) to determine paths of constant entropy, which nicely approximate the thermal profiles of simple convecting systems. The main panel in Fig. 4 shows the adiabatic profiles of MgSiO 3 melt over a range of magma ocean temperatures (color-coded by their potential temperature, which is defined at 0 GPa) for both fitted RTpress models (S11 and dK09). Despite their parameter value differences, both models show strong agreement in terms of their predicted adiabatic profiles, closely matching one another in terms of both slope and curvature over most of the geologically relevant region.
We compare our theoretical simulation-based EOS to the prediction of Mosenfelder et al. (2009), which was based entirely on laboratory shock wave experiments. Visual comparison of the adiabatic temperature gradient is shown for the 2500 K adiabat in the lower panel of Fig. 4, in terms of percentage deviation from the RTpress model. It should be noted that the EOS formalism used by Mosenfelder et al. (2009) was not thermodynamically self-consistent, due to conflicts between its heat capacity law and Grüneisen γ-model, thus affecting its prediction of adiabatic profiles. Given these constraints, we expect that although the precise pressure-dependence of the adiabatic temperature gradient might suffer some inaccuracies, the average value should remain robust. As can be seen in the figure, though the EOS model of Mosenfelder et al. (2009) deviates noticeably from our model-being shallower at low pressures and steeper above ∼70 GPa-it has an average deviation near zero. This reflects the strong agreement, to within 3%, between the adiabatic profiles predicted by both RTpress models and those of Mosenfelder et al. (2009)-as well as a later minor correction to this model by Thomas and Asimow (2013)-with an average slope of about 11 K/GPa over the mantle pressure range.
In contrast, however, we find that our EOS departs from the adiabatic trends reported in Stixrude et al. (2009), based on the first-principles MD simulations of de Koker and   (Stixrude and Karki, 2005, originally reported in ). For that model, the adiabats are systematically 25-50% steeper than our prediction throughout the lower mantle, amounting to a temperature at the core-mantle boundary that is ∼600 K hotter than the results of both our RTpress EOS and the shock-wave EOS of Mosenfelder et al. (2009). Because we do not find this same thermal behavior when fitting the RTpress model to the data of de Koker and Stixrude (2009), we know that this systematic overestimate of the temperature gradient lies with the EOS form, and is not intrinsic to the FPMD data themselves. It should be noted that this mismatch does not imply strong disagreement in all areas between these EOS models. Rather, the disagreement appears to be mostly confined to sensitive thermal properties like the adiabatic slope and Hugoniot, whereas other properties like density as a function of pressure maintain reasonable agreement, as discussed by Spera et al. (2011).
The physical significance of these predictions are seen by comparing these adiabats to the melting curve for the mantle, indicating the crystallization pathway of an early Earth magma ocean. Fig. 4 shows two estimates of the mantle liquidus in black: the dashed curve is based on a combination of MD melt simulations and high-pressure melting experiments  and the dotted line is based on direct in situ experimental measurements of synthetic chondrite (Andrault et al., 2011). The relative slopes of the liquidus and adiabat determine where the mantle will begin to crystallize. If the adiabats are shallow compared to the melting curve, the magma ocean will crystallize from the bottom up, as assumed in previous magma ocean models (Solomatov and Stevenson, 1993;Abe, 1993). In contrast, if the liquid adiabats lie tangent to the melting curve at mid-mantle depths, then first crystallization can occur well above the base of the mantle, possibly allowing a center-outwards style of crystallization. Fig. 4 shows that MgSiO 3 melt adiabats are sufficiently steep to enable either style of crystallization-initiating at either ∼80 GPa or at the base of the mantledepending on the curvature of the liquidus, indicating that further work is required to settle this issue.

Conclusions
We present a new Equation of State for liquids at high P-T conditions, the RTpress model. This EOS provides a simple and physically motivated extension of the Rosenfeld-Tarazona model to high pressure applications (Rosenfeld and Tarazona, 1998), while avoiding issues with overfitting that affected previous efforts to develop liquid-appropriate EOSs. This results in an intuitive description, which separates the liquid's behavior into coupled compressive, adiabatic, and heat-capacity models that can each be constrained by either simulated or experimental data. We apply this model to the MgSiO 3 system, obtaining an accurate fit of the MD simulations of both Spera et al. (2011) and de Koker and Stixrude (2009)-given by models S11 and dK09, respectively-which compare favorably with the experimental shock-wave data and corresponding EOS of Mosenfelder et al. (2009). The adiabatic profiles of MgSiO 3 melt are determined from the fitted model and compared to two competing determinations of the mantle melting curve to determine the style of magma ocean crystallization. We find that the MgSiO 3 adiabats are steep enough to imply center-outwards crystallization assuming the high-curvature liquidus of Stixrude et al. (2009), but still results in bottom-up crystallization with the nearly-flat liquidus of Andrault et al. (2011). Further work must be done on the high-pressure melting properties of a realistic mantle chemistry in order to resolve which of these qualitatively distinct evolutionary paths best-describes the early evolution of rocky planets.

Acknowledgments
We would like to thank Mark Ghiorso, Frank Spera, and Lars Stixrude for useful conversations and helpful criticism throughout this project. We also thank the editor Kei Hirose and anonymous reviewers, whose comments were extremely valuable in improving the manuscript. ASW was partially supported by the Turner Postdoctoral Fellowship Grant at the University of Michigan. Funding from the ETH Zürich Postdoctoral Fellowship Program enabled DJB to conduct this research.

A. The Reference Adiabat and Grüneisen Parameter
The isentropic temperature path is described in terms of a weighted integral of the thermodynamic Grüneisen parameter (defined by Eq. 1). The expression for the thermal profile is given by integration: where γ 0S (V ) is the Grüneisen parameter variation along the reference adiabat. Note that throughout this manuscript, we use adiabat and isentrope interchangeably, reflecting their equivalence for reversible equilibrium thermodynamics. Calculating the Grüneisen parameter at any arbitrary temperature (Eq. 17) uses the thermal perturbation approach to modify this reference value (see Sec. 5.2). There are numerous analytical forms available to describe the reference γ-profile, but unlike the reference compression profile, the differences between these models significantly impact the overall fit and behavior. While the power-law γ-model is both simple and popular (see Poirier, 2000), we adopt the γ-model proposed by Stixrude and Lithgow-Bertelloni (2005). The derivation is based upon a Taylor expansion of the characteristic vibrational frequency for a material in terms of finite strain, f : where V and V 0 are the molar volume and its 0 GPa reference value. The resulting model for the Grüneisen parameter is: with, a 1 = 6γ 0 and a 2 = −12γ 0 + 36γ 2 0 − 18γ 0 where the Taylor expansion coefficients, a 1 and a 2 , are expressed in terms of the value and derivative of the Grüneisen parameter at zero GPa, γ 0 and γ 0 = V 0 (dγ/dV ) 0 . Note that this is parameterized in terms of the absolute Grüneisen derivative, γ 0 = γ 0 q 0 , since it avoids the mathematical singularity at γ 0 = 0, which has been observed for some liquids such as molten quartz (e.g., Asimow, 2012).
The corresponding isentropic temperature profile is derived from Eq. A.1 by inspection, yielding: This model has much more realistic liquid γ-compression properties compared to other two-parameter models, as discussed in detail in Sec. 5.2.
To derive the full temperature-dependent expression for the Grüneisen parameter (Eq. 1), we use its thermodynamic definition: where the simplification in terms of the isothermal entropy derivative is arrived at via the cyclic relation. This required entropy derivative is obtained from the derivative of the entropy expression (Eq. 12): where γ 0S = γ 0S (V ) is the reference Grüneisen function, defined along the reference adiabat.

B. Generalized RT model details
The Generalized Rosenfeld-Tarazona model (see Sec. 3.4) relies upon a number of details that are not fully presented in the main text for clarity, but are instead included here. The temperature derivatives of the thermal deviation (Eq. B.1) are given by: where the n th -order derivative, f (n) , is expressed by the recurrence relation above. The volume derivative of the thermal coefficients is given by: As in the standard RT model, the kinetic energy contribution is fixed to the theoretical high-temperature limit: where N is the number of atoms per formula unit and k B is Boltzmann's constant. This assumption should be generally reasonable to describe atomic motions above the melting point (see e.g. de Koker and Stixrude, 2009), however it does neglect any additional possible contribution from thermally excited electrons, which play a significant role at extremely high temperatures well above the liquidus and are separately handled and discussed in App. C. The total heat capacity is given by the sum of potential and kinetic contributions: where the volume dependence is given by the RT coefficients, b(V ), and the temperature dependence is given by the thermal deviation derivative, f T .

C. Electronic Entropy
In most studies of mineral and melt properties, we can safely assume that electrons are found within their ground-state. At very high temperatures, however, electrons can be promoted into higher energy states above the Fermi Energy Level. This effect can be calculated approximately using hightemperature smearing of electronic energy levels by assuming simple Fermi-Dirac smearing function (e.g., Kresse and Furthmüller, 1996). When higher energy states are available to the electrons, this adds an additional entropy term to the system, since an individual electron can be found in either the ground state or one of the excited states. To play an important role, the thermal energy must be large relative to the Fermi energy levels for the material being studied.
The first-principles MD simulations of De Koker and Stixrude (2009) explored the role of electronic entropy in liquids across the MgO-SiO 2 binary, finding a small but non-negligible electronic entropic contribution that kicks in at high temperature and grows with increasing temperature and decreasing pressure (large volumes). De Koker and Stixrude (2009) adopts a physically motivated empirical description of the electronic effects that incorporates a step-function-like electronic contribution to the heat capacity that switches on above some characteristic electronic temperature T el : where the thermo-electronic heat capacity coefficient, ζ, and the electronic temperature are both assumed to have a power law dependence on volume: For temperatures above the electronic temperature (T > T el ), integration yields non-zero contributions to the entropy, free energy, internal energy, and pressure : To assess the relevance to planet formation, we should compare the temperature-pressure dependence to plausible magma ocean conditions. Immediately following a giant impact, the magma ocean will adopt a strongly stably stratified thermal profile (e.g. Nakajima and Stevenson, 2015;Lock and Stewart, 2017)-with much hotter material overlying cooler denser materialinitiating a complex transient period of slowed cooling. Eventually, this stably stratified mantle should cool back toward a simple near-adiabatic profile, and thus for the sake of simplicity we consider the threshold melt adiabat just before crystallization, providing a reasonable estimate of the P T conditions. Fig. 4 shows that the threshold adiabat has a potential temperature of around T ≈ 2750 ± 100 K (solid gray curve), for both of the two leading estimations of the mantle liquidus. Mantle temperatures thus vary with depth, rising from 2750 K at 0 GPa to around 4500 K at the core-mantle boundary. Due to the strong volume-dependence of the electronic entropy contribution, reported in de Koker and , this adiabatic thermal profile maintains a very small contribution for electronic entropy of S el < 0.05 k B /atom through the mantle. Electronic entropy is thus safe to ignore for most magma-ocean applications, though it may play some role during the early transient extreme temperatures just after a giant impact.
To ensure that our model can handle the full range of applications relevant to planet formation, we fit this electronic correction (using least squares minimization) to the first-principles-derived electronic entropies reported in  (2009), and for completeness, we report the fitted parameter values in Tbl. C.2. It is important to note that even though much of the training data lies above the electronic temperature threshold, the gradual onset of the electronic contribution ensures that it is entirely negligible for all of the simulated data of Spera et al. (2011) and nearly all but the highest temperature data of de Koker and . This exercise enables us to include the electronic contributions whenever they are non-negligible, as is likely the case for the transient period just after a giant impact.

D. Model-fitting Details
The RTpress model is fit to data according to the χ 2 minimization procedure described in Sec. 4. To obtain a model that retains an accurate description of the training dataset over the full pressure range, we introduce an effective pressure uncertainty term,σ P , that accounts for residuals in both pressure and volume. To derive this effective pressure uncertainty, we recognize that there exists an opposing alternative to minimizing pressure misfit, by instead focusing on volume misfit and fitting volume as a function of pressure. Given that most compression equation of state forms are not mathematically suited to this swapping of independent variables, we can approximately transform pressure residuals into equivalent volume residuals using the bulk modulus: where K T is the isothermal bulk modulus evaluated at each datapoint location. Thus, if we were to adopt σ P (V ) as our pressure uncertainty, it would return the same approximate result as directly minimizing volume misfit (while remaining computationally fast). Since our goal is to minimize both volume and pressure misfit simultaneously, we introduce the effective pressure uncertainty as an average of the two endmember approaches: where using this definition provides a χ 2 fitness metric that is exactly halfway between the volume-focused and pressure-focused endmembers. Finally, to further increase the incentive to respect the properties at ambient pressure, we also extend the training dataset by including an additional datapoint obtained from the experiment-based volumetric model of Lange (1997). This additional datapoint has a volume of 12.80Å 3 /atom at 1 bar and 1673 K. During the fit, it is assigned an uncertainty 10 times smaller than the rest of the training data, ensuring it has sufficient statistical weight to properly influence the fit. Taken together, these practical steps allow the model to retain an accurate description of the training data in both the high-pressure and low-pressure regimes.
In principle, directly applying the χ 2 minimization procedure described in Sec. 4 is sufficient to obtain a best-fit (or maximum likelihood) model, however it can sometimes run into optimization challenges in practice. This is easily handled, however, by obtaining an initial fit to each component of the model in stages, and finally refitting to obtain the global fit and parameter uncertainties. The parameters are thus roughly determined in following fitting stages: • Fit isothermal compression curve using P V T data at reference temperature T 0 • Fit individual energy-temperature isochores with standard RT model, represent volume dependence of coefficients using logarithmic polynomial • Determine entropy along isothermal compression curve by comparing observed E − V trend with predicted F − E trend from fitted isothermal EOS. Integrate up in temperature using RT coefficients to obtain adiabatic compression path for temperature.
• Fit the γ-model to reference adiabatic temperature path. After this step-wise procedure, all three of the sub-model components have been nicely constrained by each piece of the data. The final step is to then refine all parameters using joint least-squares regression, producing the The correlated uncertainties of the model are then determined using the standard uncertainty estimation method using the curvature of the χ 2  surface in parameter space: where Σ is the parameter covariance matrix and H is the "Hessian" or curvature matrix in parameter space. The 1-σ uncertainties, reported in Tbl. 1, are then given by the square-root of the diagonal elements of the covariance matrix, σ i = √ Σ ii . Additionally, the correlation coefficients, which express the uncertainty trade-offs between the parameters, are determined by scaling the covariance matrix: where the correlation coefficient between parameter i and j takes on values −1 < f corr,ij < +1, expressing how strongly deviations from the best fit correlate for the two parameters. The correlation matrix for the RTpress model of MgSiO 3 is shown in Fig. D In the main text, we focus primarily on the results of the fit to the high resolution data of Spera et al. (2011), but we also explore the consequences of fitting the FPMD data of de Koker and . Though the best fit parameters for each model are fairly different (see Tbl. 1), the resulting equations of state remain in general agreement in terms of their high pressure properties. This is made clear by the near complete agreement between their predicted adiabatic profiles shown in Fig. 4. It should be noted that although the dK09 model has considerably higher uncertainties due to the smaller and noisier training dataset, it possess properties that are generally more accurate at near-surface pressure conditions when compared with models based on 1-bar experiments (Lange, 1997;Ghiorso, 2004;Ai and Lange, 2004). For completeness, we include the fitness diagrams associated with the dK09 model as well.

E. Shock Hugoniot
To calculate the theoretical shock Hugoniot for material undergoing melting, we follow a modified form of the standard procedure outlined by Mosenfelder et al. (2007). The Rankine-Hugoniot expressions impose conservation of energy, mass, and momentum on the shocked material, providing the link between changes in pressure ∆P H , volume ∆V H , and internal energy ∆E H for points lying along the shock Hugoniot: To calculate the model Hugoniot, we must determine the internal energy change from the initial solid state at ambient conditions up to the molten high density final state. To calculate this energy change, Mosenfelder et al. (2007Mosenfelder et al. ( , 2009) consider a metastable melting transition at ambient pressure and temperature, followed by heating and compression up to the final conditions. While technically correct, this approach relies on the liquid equation of state remaining accurate at the highly metastable conditions of room temperature and pressure. Instead, we follow an adjusted procedure where the material is heated up to its melting point and then fully melted at ambient where E liq (V, T ) is the internal energy of the liquid given by its equation of state, T fus is the ambient pressure melting point (or fusion temperature), and V fus is the volume of the liquid at the melting point. The transformation energy, ∆E tr , is the energy associated with heating and melting the initial starting material at ambient pressure: where the heating energy, ∆E heat , and the energy of fusion, ∆E f us are determined from experiments. Using this method, the liquid equation of state is only ever used within its stability range, improving accuracy and safely avoiding any potentially unphysical behaviors in the strongly metastable regime.
To determine the Hugoniot, we must first collect the experimental data on the conditions of melting. Thieblot et al. (1999) provide a detailed study of the heat capacity of enstatite all the way up to its congruent melting point. They find that MgSiO 3 has a fusion temperature of T f us = 1816 K. The heating energy, given as the energy required to heat enstatite up to the melting temperature from 300 K is 192.26 kJ/mol. Richet and Bottinga (1986) provides data on the fusion energy, which is reported as 73.2 kJ/mol for enstatite and 31.1 kJ/mol for enstatite glass (reduced by the vitrification energy). The final values of the transformation energies in atomic units are thus: ∆E tr (En) = 0.550 eV/atom ∆E tr (En-glass) = 0.463 eV/atom (E.4) Underlying this calculation is an approximation that the heating energy up to melting is similar for both enstatite crystal and enstatite glass. Given these data values, the Hugoniot is then calculated by solving for the temperature, at each liquid volume, that satisfies Eq. (E.1). This procedure is applied for both the RTpress model from this work as well as the EOS model from Stixrude and Karki (2005); de Koker and  to produce the Hugoniot comparison in Fig. 3. We compare our model Hugoniot curves with the available gas gun shock-wave data presented by Mosenfelder et al. (2009), which also includes collected experimental results from a variety of previous authors (Akins, 2003;Akins et al., 2004;Luo et al., 2004b;Simakov and Trunin, 1973). These are not the only available shockbased measurements, however, as Spaulding et al. (2012) and Bolis et al. (2016) both used the decaying laser-shock method to study MgSiO 3 melt at extreme pressures. We do not include these data in our model assessment since they are focused on ultra high pressures primarily relevant to super-Earths (generally ≥200 GPa and upwards of 800 GPa at temperatures as high as 20,000 K). Though our current focus is primarily on the terrestrial magma ocean regime, future extension of this equation of state to such extreme pressures and temperatures could benefit from a full consideration of these data.