The pressure dependence of thermal expansivity affects mineral density at pressure and is an extrapolator for calculating self-compression adiabats of a self-gravitating body. I review different models for the pressure dependence of expansivity and how to decide which performs best. A finite strain model, proposed here, performs better when used to calculate adiabatic temperature lapses in both the solid silicate and liquid metal parts of a planet than either an ad-hoc exponential dependence on pressure or a commonly used mineral physics model. Choosing a particular thermal expansivity pressure dependence leads to significantly different temperatures in planetary interiors, and to inferred subsolidus properties related to homologous melting temperature. In particular, thermal expansivity in liquid metal in planetary cores at pressures comparable to Earth's core is significantly affected. The universality of the parameterization provides a simple way to model rocky planet interiors in our solar system and exoplanet interiors.

Equation of state, thermal expansion, adiabatic gradient, planetary accretion, finite strain

Planetary accretion is the process by which a planet grows from a nucleation site in the nebular dust and gas disk surrounding a young star into a self-gravitating body in orbit around the star. The nascent planet grows through stages governed by the dominant forces driving accretion: adhesive, electrostatic, and then gravitational (Armitage 2010). The growing planet matures from a planetesimal, to an embryo stage, and finally to a planet (Righter and O'Brien 2011). The nucleation site in the compositionally zoned disk controls whether the evolved planet is dominantly gaseous or rocky. After a rocky planetesimal reaches a state where it warms sufficiently, whether heated by short-lived radioactivity, by impact heating, or by adiabatic heating due to the internal pressure increase, it differentiates into metal—core and silicate—crust and mantle.

The details of the differentiation process rely on knowledge of the internal temperature structure of the growing planet. When bodies are small, thermal diffusion dominates and the disk temperature and short-lived radiogenic element abundance control the planetesimal's temperature (Šrámek et al. 2012). After planets grow sufficiently large to differentiate, solid-state convection in the silicate mantle and liquid state convection in the metallic cores govern the thermal structure (Breuer et al. 2010). These are essentially adiabatic temperature profiles set by the conditions at the convective boundary layers (the surface or the core-mantle boundary). Because the thermal expansivity (α) along with gravity (g) and heat capacity (CP) are involved in the calculation of the adiabatic gradient,
(1)
an accurate description of α's pressure dependence is needed to describe the temperature. The behavior of g with radius, in contrast, is simply parameterized (essentially two linear segments; see Fig. 1) and the pressure dependence of CP is small enough to be neglected if the mineralogy is not known (Appendix1).
For a given mass, a planet's size is governed by its density structure. In turn, the density is set by the proportions of the planet's constituent minerals and the equation of state (EoS) of those minerals. Because α's definition is
(2)
it represents the variation in volume (V) or density (ρ) with temperature

There is a difference between α's role in Equations 1 and 2. The thermodynamically astute reader will recognize a fallacy in this claim, and indeed there is: through the unity of thermodynamic relations, α is the same property in Equations 1 and 2. In Equation 1 however, α need not represent any real object. An example is the hard-sphere liquid (Hansen and McDonald 2013). Its free energy may be written explicitly (Lee 1995) and its thermal expansivity calculated from derivatives of the expression with respect to pressure and temperature. However, no experiment can measure α by heating a hard-sphere liquid and measuring its change in volume, which is the natural interpretation of Equation 2.

In Equation 1, α represents a pressure dependent bulk property of the material and can simply be a suitably chosen function of ρ or r that reproduces an adiabatic planetary density profile such as PREM (Dziewonski and Anderson 1981). An adiabat calculated that way might also be compared with a melting curve for metal or peridotite to determine melting conditions to assess whether a magma ocean might arise or a core might segregate in a growing planet, such as Labrosse et al. (2015) did. In convection modeling, α governs the buoyancy force arising from temperature variations in the bulk convecting fluid, liquid or viscous solid (Turcotte and Schubert 2002). Driscoll and Olson (2011), for example, used a pressure-dependent bulk α in their study of magnetic field strength around exoplanets, where the material was classified as iron, peridotite, perovskite, and post-perovskite.

In contrast, α is an intrinsic property of a mineral obtained through measurement of V vs. T and modeled with Equation 2 and then incorporated as part of an EoS. As an example, Stixrude and Lithgow-Bertelloni (2011) built a detailed mineralogical model of the mantle and calculated thermal expansion of the various assemblages met along P-T trajectories through it, leading to a detailed, and discontinuous description of the material.

If used to represent a bulk property, α might not ever represent a value for any particular mineral or mineral aggregate. Moreover, in the absence of knowledge of the constituent mineralogy of, say, an exoplanet, α's pressure dependence captures the mineralogical tendency to adopt denser forms at higher pressures in a general way. Thus, the need to parameterize self-compression and mineral behavior lead to different α model choices, which is the subject of this article.

Material equation of state (EoS)

To model the stages of planetary accretion of a rocky planet, a simple material parameterization is desirable, essentially due to one's ignorance of the identity of the specific materials and of their proportions. The two basic constituents are metal and silicate that I treat as single component phases in the thermodynamic sense. For computational simplicity I use a polythermal Murnaghan equation of state for each because it can be evaluated in closed form for ρ(P,T), the density at a particular pressure and temperature. Explicitly,
(3)
with ρ0 a density at P = 0 and reference temperature T0, K is the isothermal bulk modulus at P = 0, and T0 and K′ is its pressure derivative. Iα represents the integrated thermal expansion effect on density from the reference density, ρ0. Again, for simplicity, I assume that dα/dT is zero [a high-temperature—high-pressure approximation (Chopelas and Boehler 1989)], but that α is pressure dependent. Hence one can integrate Equation 2 to define
(4)

Pressure dependence of thermal expansion

The decrease of thermal expansivity with increasing pressure is well established observationally and theoretically (Chopelas and Boehler 1989; Anderson et al. 1992). One simple way to parameterize this is through an exponential decrease with increasing pressure (Tosi et al. 2013). Using the material bulk modulus as an internal pressure scale, one can write
(5)
with α′ the scaled rate of pressure decrease from the zero pressure value α0. If, say, α decreases to 50% of its ambient pressure value at the CMB [P = 135 GPa; (Stacey 1992)] then α′sil = log(2) × 135/Ksil, and for metal, α′met = log(2) × (360 – 135)/Kmet. Table 1 lists these parameters.
An alternative parameterization (Chopelas and Boehler 1989; Anderson et al. 1992) is to relate the pressure dependence to the volume change on compression V/Va. Chopelas and Boehler (1989) proposed
(6a)
with δ = 5.5 ± 0.5, whereas a generalized version of this is (Anderson et al. 1992; Wood 1993),
(6b)
with δ0 = 6.5 ± 0.5 and κ = 1.4. These forms lead to either a power law (Eq. 6a) or exponential dependence (Eq. 6b) on volume,
(7a)
or
(7b)
The equivalence of Equations 6a and 6b at small compressions may be seen by letting V/V0 = (1 – ε). Then from Equation 7a, (V/V0)δ ≈ (1 – δε). From Equation 7b, (V/V0)κ – 1 ≈ –κε and exp[(δ0/κ)(–κε)] ≈ 1 – δ0ε. Hence the two forms are identical for small compressions if δ ≈ δ0, and Equation 6b offers more control over extrapolation to higher compressions through κ. Table 1 contains the values used.
A final alternative for α's pressure dependence recognizes the similarity of the dependence on V/V0 to the finite strain parameter f = ½[(V/V0)–2/3 – 1] (Birch 1952). Thus one can also relate α to f (Driscoll and Olson 2011):
(8)
where ϕ is some positive, monotonically decreasing function of f. Lest this characterization be too vague, the particular choice used here is
(9)
In their planetary modeling, Driscoll and Olson (2011) used a simpler expression, ϕ(f) = (1 + 2f)–9/2.

Planetary P, T, and g profiles

To show the consequences of different choices for the pressure dependence of thermal expansivity, one needs to calculate consistent pressure (P), temperature (T), and gravitational acceleration (g) profiles. For a given planetary mass, I take the silicate and metal masses proportional to those in the Earth (Table 1). Either a differentiated profile may be calculated from the metal and silicate equations of state, or an undifferentiated profile may be calculated from a mechanical mixture of the components. A consistent P-T profile is obtained iteratively from initial conditions assuming separate adiabatic profiles in the mantle and in the core, or a single adiabatic profile if homogeneous. Iteration stops when the fractional change in the body's gravity and radius is <10–5.

One or two temperature fixed points are specified for each profile: the temperature at the surface and, if differentiated, the temperature at the CMB. Given the mass of the planet M, calculating the P, T, and g profile involves these steps:

  1. Set P(r) = 0, T(r) = constant (mantle and/or core).

  2. Calculate radius of the CMB and planet R with prevailing P(r), T(r) by integrating dM/dr = 4πr2ρ(P(r),T(r)).

  3. Calculate g(r) = G(Mr/r2), where Mr is the mass within radius r.

  4. Using the identity dP/dr = –g(r)ρ(P(r),T(r)), calculate a new “cold body” pressure profile P(r)=

    ∫rR
    ρ(P(r),T(r))g(r)dr.

  5. Calculate a “cold body” T(r) using the adiabatic gradient (Eq. 1) fixed at the conditions of the surface (and if differentiated, the CMB).

  6. Compare to previous R and g(R); if fractional change <10–5, profile is converged.

  7. Not yet converged; return to step 2 with new “warmer body” P(r), T(r), and g(r).

The algorithm typically converges within 5 to 10 iterations. With the values in Table 1, and with an adiabatic profile initiated at the surface at 1623 K [a characteristic shallow mantle temperature (Parsons and Sclater 1977; Stein and Stein 1992)] and continuous with a core adiabat at the CMB, the planetary radius, core radius, and gravity are within 0.1% of the Earth (Dziewonski and Anderson 1981). Figure 1 shows a comparison with calculated P and g profiles for the Earth.

The choice of a finite strain-based model for the pressure dependence of α is not immediately obvious. My assessment process involved a suite of plausible formulas for ϕ(r) (Fig. 2). The simplest formulas do not decrease fast enough through the mantle and core range of f to reproduce the tabulated decreases compiled from geophysical sources (Stacey 1992). I found through experimentation that a product of monotone decreasing functions, exemplified by Equation 9, fit the trends best for both metal and silicate. Relative to that, the mineral physics parameterizations asymptotically flatten quickly with increasing strain. The consequences of this behavior will become clear once the various models are used to compute adiabats.

Figure 3 shows T profiles due to adiabatic heating. In all cases an Earth-mass ME planet with a fraction of metal to silicate ~0.32 is used (Table 1). Temperature at the surface is 1623 K and at the CMB is 4000 K. Unlike Figure 1, temperature is not forced to be continuous at the CMB; rather, the CMB temperature is the foot of a new adiabat. I also show two peridotite solidus curves, one as parameterized by Wade and Wood (2005) and the other by Fiquet et al. (2010) (Table 1). The planetary surface and CMB radii are slightly different given the different α parameterizations.

The slopes of the adiabatic curves all approach zero at the center of the Earth, due to the adiabat's dependence on g(r), which is zero there (see Eq. 1). However, even though the temperatures at the CMB are identical, the temperatures at the center are quite different as are the slopes of the curves. For the same CMB temperature and approximately the same core radii, the temperatures at the center are 4388, 5616, and 7334 K. Clearly, the choice of the thermal expansivity's pressure dependence is important when phenomena relative to an adiabatic temperature gradient are involved.

The methods yield notable temperature differences at the mantle side of the CMB. The adiabats projected using the finite strain and mineral physics models yield much lower temperatures. One would conclude from the low temperatures there that a significant thermal boundary layer would develop, driving convection in the mantle by bottom heating. In contrast, the degree of basal heating with the exponential model would be smaller, with a correspondingly lower potential to drive convection.

Another difference between the adiabatic trajectories are their curvatures in the mantle. The mineral physics and finite strain adiabats are quasi-linear there. However, the exponential adiabat is subtly concave upward. Figure 4 displays the mantle portions of the three curves relative to the peridotite solidus to highlight this behavior and its consequences. If an ~500 K warmer foot for the adiabat were chosen, the exponential model for the adiabat would intersect the solidus at two radii. Two solidus crossings would suggest that zones of melt could form at both the base of the mantle and at the surface, leading to a basal magma ocean (Labrosse et al. 2007). The other models would yield melting at outer planetary radii, or a surface magma ocean.

As a way of choosing, which is the preferred parameterization, I recruit another thermodynamic expression for the adiabatic lapse (Stacey 1992):
(10)
with γ the thermodynamic Grüneisen parameter and KS the adiabatic bulk modulus. Because γ in the outer core is a virtually constant value, 1.52 (Alfè et al. 2002), use of VP in the core liquid, along with g(r) calculated from PREM, provides a test for which model best describes compression in the core, and, to a lesser extent, the mantle. The models (Fig. 5) are of Earth-mass planets with a surface adiabat initiated at 1623 K and a CMB adiabat initiated at 4000 K. The comparison with PREM shows that the finite strain model for α most closely reproduces PREM's adiabat in the core liquid. The situation in the mantle is not as easily compared due to the phase transitions in upper mantle minerals and the material being polymineralic. Restricting the comparison to the lower mantle, where the mineralogy changes little, the finite strain and mineral physics models perform equally well compared to PREM, with γ ≈ 1.5. The lower mantle range for the Grüneisen parameter is 1 ≤ γ ≤ 1.4 based on γ estimates and the adiabatic lapse (Brown and Shankland 1981; Jackson 1998; Katsura et al. 2010; Stixrude and Lithgow-Bertelloni 2011). The α models yielding a comparable lapse lie marginally beyond the high end of the range.

The mineral physics based model and the finite strain model perform equally well in the silicate mantle, but the performance is notably poorer in the core for the mineral physics based model. It is worth asking whether the poor performance in the core is due to the choice of particular values for the parameters used, or due to an inappropriate physical model. To answer this, I determined the parameters α0, δ0, and κ that fit PREM's adiabatic lapse in the core best. They are α0 = 4.98 × 10–5K–1, δ0 = 2.15, and κ = 1 × 10–4. While α0 is indistinguishable from the value in Table 1, the κ value shows that Equation 6a is a better model for an Earth-like core than is Equation 6b—a surprising result for an improved physical model (Anderson and Isaak 1993). Moreover, δ0 is significantly different than its range of 4–6 for silicate minerals (Anderson et al. 1992), and differs from K′, deviating from the rule of thumb that δa ≈ K′ for silicates (Anderson et al. 1992). The poor performance of the mineral physics based model may not be surprising if one reflects that liquids and solids differ in their internal structure and thus the interatomic forces that give rise to α, and K and K′. However, it underscores the advantage of the finite strain model: it captures the properties of both solids and liquids simply and uniformly.

The adiabatic profiles, while they yield Earth-like surface and CMB radii and gravity are not very accurate density models everywhere. Compared to PREM (Fig. 6), the density is overestimated in the shallow mantle by up to 30%. Core densities are within ±2% in the outer core and the density gradient is close to PREM, but there is no provision in the model for a solid inner core and hence the densities are underestimated. Upper mantle densities are not particularly well described due to the transition zone phase changes that affect both the temperature structure and the density (Katsura et al. 2010; Stixrude and Lithgow-Bertelloni 2011). In most of the planet, however, the density profile is within ±5% of PREM's.

The α models explored here focused on three aspects of the resulting adiabatic profiles: (1) their convexity; (2) their temperature lapse; and (3) their approximation to the known density profile of the Earth.

All of the models yield Earth-like dimensions, gravity, and maximum pressures for Earth-mass objects that have Earth-like metal/silicate ratios. Of the three parameterizations, however, the finite strain-based choice yields an adiabatic lapse most closely resembling Earth's in both the silicate and the metal parts of the planet (Fig. 5). This is established though comparison with PREM and the independently known behavior of the thermal Grüneisen parameter γ. The finite strain model matches the core's properties best, and performs as good as the mineral physics-based model in the mantle. A variant finite strain model used by Driscoll and Olson (2011) is not as successful, showing that some care in choosing ϕ(f) (Eq. 8) is warranted.

The exponential model, though intuitive and mathematically and computationally straightforward (Tosi et al. 2013), has an undesirable curvature in a T–r plot (Fig. 4). The character of the curvature could lead to false inferences about magma ocean development and to inferences of homologous melting temperature that control silicate rheology and seismic attenuation (Stacey 1992). The mineral physics model, despite its solid theoretical and observational underpinnings, leads to a temperature lapse that is too low in the core (Figs. 2 and 3).

None of the models accurately reproduce density throughout the mantle and core (Fig. 6), mainly because in their need for simplicity the EoS used neglects the solid-solid phase transitions that characterize the compression of the shallow mantle. Once into the lower mantle, however, they yield densities that are ±5% of PREM densities and thus do nothing outré given our knowledge of material behavior. Whether or not the profiles match PREM's density is unimportant when used for estimating the conditions of exoplanets, when only mass and radius is known (Howard et al. 2013). The simple metal+silicate model reproduces Earth's gross properties well (Fig. 1).

One could imagine further efforts to improve an α model by relaxing the high-temperature–high-pressure approximation and incorporating a nonzero temperature derivative, or, indeed, a Suzuki-type Debye model for thermal expansion (Suzuki 1975). Whether the added complexity is warranted to improve the performance for the silicate planetary component is not obvious. The virtue of the approach advocated here is that it is implemented in a simple way and can be incorporated into planetary accretion modeling without undue computational burden.

The adiabatic gradient's definition involves α, but the implications of a particular choice for α's pressure dependence on the gradient's behavior are not immediately obvious. Even mineral physics-based forms might not accurately represent bulk material behavior. Different forms lead to unexpected curvature in self-compression profiles and to significantly different adiabatic temperature lapses, potentially leading to unwarranted inferences for melting, freezing, and phenomena linked to homologous temperature.

The finite strain model for α's pressure dependence fits Earth's adiabatic lapse the best and appears equally suited to silicate solids and metallic liquids. Modelers of exoplanet compositions and internal structure could benefit from the uniformity and simplicity of the formulation. On account of the higher thermal expansivity in planetary cores that the finite strain model prescribes, the role of thermal buoyancy in numerical dynamo simulations may need to be reassessed.

1Deposit item AM-17-86007, Appendix. Deposit items are free to all readers and found on the MSA web site, via the specific issue's Table of Contents (go to http://www.minsocam.org/MSA/AmMin/TOC/2017/Aug2017_data/Aug2017_data.html).

I thank Jamie Connolly and Matthieu Laneuville for detailed comments on an early draft, Dave Yuen for enthusiastic encouragement and insights, and Marine Lasbleis for comments on and interest in the effort. The reviewers, Don Isaak and Scott King, provided valuable comments that improved the manuscript's style and logic. Work partly funded by JSPS KAKENHI Grant no. 15H05832.