Lower mantle geotherms, flux, and power from incorporating new experimental and theoretical constraints on heat transport properties in an inverse model

An inverse method is devised to probe Earth's thermal state without assuming its mineralogy. This constrains thermal conductivity (κ) in the lower mantle (LM) by combining seismologic models of bulk modulus (B) and pressure (P) vs. depth (z) with a new result, ∂ln(κ) $/$ ∂P ∼ $\mathrm{7.33}/B_T$, and available high temperature (T) data on κ for lengths exceeding millimeters. Considering large samples accounts for the recently revealed dependence of heat transport properties on length scale. Applying separation of variables to seismologic $\partial B/\partial P$ vs. depth isolates changes with T. The resulting LM dT $/$ dz depends on ${\partial }^{\mathrm{2}}B/\partial P^{\mathrm{2}}$ and $\partial B/\partial T,$ which vary little among dense phases. Because seismic $\partial B/\partial P$ is discontinuous and model dependent ∼ 200 km above the core, unlike the LM, our results are extrapolated through this tiny layer (D^′′). Flux and power are calculated from dT $/$ dz for cases of high (oxide) and low (silicate) κ. Geotherm calculations are independent of κ, and thus of LM mineralogy, but require specifying a reference temperature at some depth: a wide range is considered. Limitations on deep melting are used to ascertain which of our geotherm, flux, and power curves best represent Earth's interior. Except for an oxide composition with miniscule ${\partial }^{\mathrm{2}}B/\partial P^{\mathrm{2}}$, the LM heats the core, causing it to melt. Deep heating is attributed to cyclical stresses from > 1000 km daily and monthly fluctuations of the barycenter inside the LM.

1 Introduction and background

Heat moves when a temperature (T) difference exists, where the net flow is from hotter to colder regions. This phenomenon is important to Earth because it is dynamic. But as a consequence, the outcome of a laboratory experiment is greatly influenced by the time dependence of the applied heat (e.g., Tye, 1969), which has led to overlooking the length-scale dependence of heat transport and misunderstandings of experimental limitations and uncertainties as well as of microscopic mechanisms (Hofmeister, 2019, 2021). Thermal models in Earth science are particularly affected by these shortcomings, due to wide variations in relevant length scales, temperature, pressure, and material properties, such as transparency to thermal radiation. An improved understanding, based on a new theory and accurate data on mineral heat transport, is described next.

1.1 Recent findings on heat transport properties relevant to mantle studies

One incorrect presumption is that the physical properties representing heat flow (thermal conductivity, κ, or its close relative thermal diffusivity, D) are independent of the distance along the thermal gradient. This static view is inconsistent with Fourier's heat equation, as follows. Its simplest one-dimensional form is

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial T}{\partial t}}=D{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}},\end{array}$$

where t = time and z represents the direction of heat flow. Equation (1) holds for temperature changes being sufficiently small that the relevant properties vary negligibly. Dimension analysis of this effectively constant T condition provides

$$\begin{array}{}\text{(2)}& D\sim L^{\mathrm{2}}/\mathit{\zeta }\text{or}D\sim uL\text{where}u\sim L/\mathit{\zeta },\end{array}$$

where L is the distance over which heat travels, ζ is a time constant, and u is a characteristic speed. Thermal conductivity likewise depends directly on L because it is proportional to D,

$$\begin{array}{}\text{(3)}& \mathit{\kappa }=\mathit{\rho }{c}_{P}D\equiv CD,\end{array}$$

and because the multiplying parameters, density (ρ) and specific heat (c_P) at constant pressure (P), are independent of L.

The required length-scale dependence has been masked by experimental limitations, including ubiquitous use of similar sample lengths of > 1 to < 5 mm. Many experiments are steady state, so time and ζ are irrelevant. Other common techniques are periodic, where these oscillations about quasi-equilibrium involve another, different time constant (see, e.g., Tye, 1969; Zhao et al., 2016). The transient technique of laser-flash analysis (LFA), which avoids heat losses from physical contacts (see Vozár and Hohenauer, 2003) and monitors thermal evolution across L with time, has confirmed that D linearly depends on L below about 1 mm for electrical insulators, glasses (Hofmeister, 2019, chap. 7), semi-conductors, metals, and alloys (Hofmeister, 2021). Results (Fig. 1) are consistent with a linear response when L is small, as in Eq. (2).

Misunderstandings also stem from reliance on the historic kinetic theory of gas (KTG) to depict heat transfer in solids. However, heat and matter move together across long expanses in a gas, which is unlike a solid where these motions are decoupled. Furthermore, KTG assumes elastic collisions, through which temperature cannot change. Neither non sequitur is addressed by morphing molecular collisions in a gas into elastic scattering of pseudo-particles denoted as phonons in a solid. Because gas data are collected under negligible T gradients to avoid convection, assuming random fluctuations in all three directions is reasonable and provides formulae mostly compatible with gas data. Yet, the ratios of the transport properties are not correctly described, while the ubiquitous emission of thermal radiation from all states of matter remains unexplained. Accounting for inelasticity in molecular collisions addresses both shortcomings in KTG (Hofmeister, 2019, chap. 5).

Regarding condensed matter, Fourier assumed heat flows into, across, and out of the stationary solid, whereby part of the heat is stored in the elements along the path. The process is diffusion, which is underscored by Fick constructing his formulation after Fourier's.

Fourier defined flux as heat per area per time and realized that

$$\begin{array}{}\text{(4)}& \mathrm{\Im }=-\mathit{\kappa }(T,P){\left.{\displaystyle \frac{\partial T}{\partial z}}\right|}_P=\mathit{\kappa }(T,P)\times \left(-{\left.{\displaystyle \frac{\partial T}{\partial z}}\right|}_P\right).\end{array}$$

One dimension suffices for discussion since heat flows down the thermal gradient per the second law of thermodynamics. Equation (4) is fundamental: taking its spatial derivative, conserving energy, and simplifying using the definition of Eq. (3) leads to Eq. (1).

Figure 1 Dependence of thermal diffusivity at 298 K on thickness, for various samples, as labeled. Fits are least squares. (a) Insulators with simple compositions. Parallel and perpendicular orientations are shown for quartz. After Hofmeister (2019) – this is the revised and augmented version in Hofmeister (2021). (b) Nearly pure elements and two alloys. At low thickness, the trends are linear and converge. After Hofmeister (2021), with a Creative Commons license.

Experiments and theory show that light is the diffusing entity in solids (Hofmeister et al., 2014; Criss and Hofmeister, 2017), which is real and pure energy. Light, unlike a phonon, crosses interfaces. Attenuation of light across the sample provides the length-scale dependence of Fig. 1. These recent discoveries led to new formulae for the dependence of κ on P and T and the absorption spectra of a material, which were verified against reliable data on κ below 2 GPa (Fig. 2) and on D and κ from a few kelvins to well above ambient T (Hofmeister, 2019, 2021). LFA measurements of D at high T (Fig. 3a), combined with Eq. (3), show that κ above 1000 K is nearly constant for structures or chemical compositions more complex than Al₂O₃ (Fig. 3b). These advances are used in the present paper to evaluate of thermal conductivity in the lower mantle (LM) while accounting for ambiguities in the temperature and mineralogy for this immense region of the Earth.

1.2 Reliability of available information on lower mantle heat transport

Thermal models are based on transport properties. To achieve high pressures appropriate to the deep Earth, diamond anvil cell (DAC) experiments probe tiny samples. Accurately determining P near 1 atm in devices geared for extreme compression has not been achieved. Hence, results from DAC heat transport experiments are benchmarked against independent measurements of D or κ at ambient P (e.g., Hsieh et al., 2009). However, ambient data are collected from L > 1 mm, which are ∼ 100× larger than sample thicknesses used in DACs. Extrapolation from large to small L was not done and is non-linear (Fig. 1). Very high P studies are difficult, leading to additional problems, as discussed in detail by Hofmeister (2009, 2010b, 2019, 2021). To summarize, large thermal gradients preclude use of Eq. (1) while requiring knowledge of the T dependence of D (or κ) at P, which is the unknown sought. For tiny samples, heat flow is two-dimensional but one-dimensional equations are used. At high T, cooling occurs by ballistic radiation to the surroundings (e.g., to the detector used to ascertain T), which is not addressed in Fourier's description of heat diffusion (conduction). Thermal gradients changing direction during the experiments of McWilliams et al. (2015) and Konôpková et al. (2016) (see figure 13 in Hofmeister, 2021) were not addressed in their analysis. Thermoreflectance methods (e.g., Hsieh et al., 2009) assume the length scale over which heat diffuses, which dictates the results per Eqs. (1) to (4). Low P experiments using ∼ millimeter lengths and other techniques utilized in Fig. 2 lack these difficulties, as discussed in previous work and Sect. 3, and are utilized here.

Importantly, thermal gradients inside Earth are low. Even in the lithosphere, $\partial T/\partial z$ only reaches ∼ 20 K km⁻¹. Thermal transport properties vary little over a few degrees (Fig. 3). Because T varies less than 4 K over L > 5 km inside Earth, its thermal length scale is immense, and so heat transfer therein is always diffusive and isothermal properties are relevant. But in the laboratory, L is ∼ 10⁵ smaller, permitting ballistic (boundary-to-boundary) transport to augment diffusion, as recognized in minerals and rocks by Kanamori et al. (1968) and further documented by Pertermann and Hofmeister (2006), Branlund and Hofmeister (2007), and Merriman et al. (2018). Laser-flash experiments reduce and remove ballistic effects via sample coatings and via models (e.g., Blumm et al., 1997; Hahn et al., 1997). Our large and growing LFA database (e.g., Hofmeister, 2019) and the associated theoretical model are essential to ascertain heat transport at high mantle temperatures.

Seismic models provide velocities and density inside the Earth. Pressure is well constrained, since Earth's mass and moment of inertia provide independent boundary conditions (e.g., Anderson, 2007). Mineralogy is based on comparing laboratory data on minerals to radial models, such as the preliminary earth reference model (PREM) of Dziewonski and Anderson (1981) since radial changes depict average values. Comparison with laboratory studies in a forward-(fitting) approach is used but leads to equivocal results because temperature is not known independently. For the Earth, temperatures are changing, heat is moving, and seismic waves contribute energy to the rocks during their attenuation. Hence, conditions are not adiabatic, as previously assumed in forward-(fitting) models. In addition, minerals vary greatly in possible chemical compositions and structures.

Assessing the lower mantle is particularly uncertain because no rocks have been exhumed from below 670 km. Inclusions in diamonds only indicate P and T conditions when a single inclusion contains multiple phases because most, if not all inclusions, predate their diamond host (Nestola et al., 2017). The inference that lower mantle material is preserved in microdiamonds is based on separated inclusions of (Mg,Fe)O and enstatite (Stachel et al., 2000). The tetragonal garnet phase TAPP (now jeffbenite) once considered to form in the lower mantle is now known to be stable above 13 GPa, i.e., in the transition zone (Nestola et al., 2016).

Evaluating temperatures from seismic models via forward-fitting requires knowledge of the mineralogy (e.g., Cammarono et al., 2003). A thermal model is needed to account for Earth's heat being lost to space (i.e., non-adiabatic gradients). For the lower mantle, a wide range of κ values is possible due to the ambiguities in mineralogical models, even if experimental uncertainties were small. An alternative approach to fitting seismic velocities is needed to better understand this immense region of the thermally evolving Earth.

Figure 2 Measured pressure derivatives of thermal conductivity as a function of sample compressibility (the inverse of bulk modulus). All data were collected from millimeter-sized samples below 2 GPa. The box lists parameters for a linear fit with no intercept. Metals and Si (red crosses) are included in the fit. KBr exemplifies hydroscopic alkali halides which are soft (bulk moduli < 16 GPa). Details on the 24 heat transport studies of this figure are given in Hofmeister (2021), Table 3. Bulk modulus has been measured many times for the samples and is listed in the compilation of Bass (1995) and several others.

Figure 3 Heat transport properties of selected dense solids with a logarithmic y scale. (a) Thermal diffusivity on single crystals from laser-flash analysis, measured up to ∼ 1800 K. Fits using three parameters are shown. Diamond from Hofmeister et al. (2014). Oxides from Hofmeister (2012). Forsterite and olivine are 001 plates (Pertermann and Hofmeister, 2006). Perovskites from Hofmeister (2010a). (b) Thermal conductivity calculated from data on D, c_P, and density using Eq. (3). For the oxides, the dotted lines use constant, high-T values for c_P and ρ, whereas their other curves use temperature-dependent data. Orthorhombic perovskite is an estimate for a wide range of compositions; see Sect. 3.1.

1.3 Purpose and thesis

In view of limited knowledge of the LM, an analytical inverse approach is used here to decipher its thermal state from a seismic reference model with minimal assumptions. As in previous large-scale mineralogical or thermal studies (e.g., MacDonald, 1959; Anderson, 2007; Murikami et al., 2009; Criss and Hofmeister, 2016), average, radial temperatures are sought to describe Earth's structure, which is reasonably represented as spherical shells. Surface heat flux being remarkably similar for the continental and oceanic crusts, despite the great contrast in their heat-generating elements (e.g., Veiera and Hamza, 2018), points to the radial thermal gradients dominating Earth's thermal state and evolution. The low measured surface emissions of ∼ 60 mW m⁻², corresponding to ∼ 100 W km⁻³ of underlying rock, show that thermal evolution is now slow; i.e., conditions are quasi-steady state. Hence, angular (lateral) motions of heat are unimportant to describing Earth as a whole.

Our mathematical analysis is based on decades of mineral physics efforts which show that (1) pressure derivatives of diverse physical properties vary far less than ambient values do and that, as P climbs, all properties increase more weakly with P. (2) Physical properties at high T behave similarly, as illustrated by the Dulong–Petit limit representing heat capacity at high T. (3) Second-order P or T derivatives of physical properties are small, which means that cross derivatives are small, and so separation of variables reasonably describes many physical properties. Tabulated data on diverse properties and phases (e.g., Anderson and Isaak, 1995; Bass, 1995; Fei, 1995; Knittle, 1995) illustrate these points.

If variables are separable, a property of interest (Υ) follows the form

$$\begin{array}{}\text{(5)}& \mathrm{{\rm Y}}(T,P)={\mathrm{{\rm Y}}}_{\mathrm{ambient}}f\left(T\right)g\left(P\right),\end{array}$$

where f and g are independent, dimensionless functions. Equation (5) describes bulk and shear moduli from diverse elasticity experiments (e.g., Anderson and Isaak, 1995; Bass, 1995). This finding is important to heat transport, as bulk modulus is the prime descriptor of κ(P) per dimensional analysis (e.g., Dugdale and MacDonald, 1955).

For the lower mantle, variations in velocities from available seismic reference models differ negligibly except for the ∼ 200 km above the core (D^′′) where variations among studies are small, despite larger uncertainties for this region (see figures in Kennett et al., 1995). Utilizing PREM suffices (see Section 2.1 for further discussion). In the LM, excluding D^′′, velocity changes are slow and smooth, leading to interpretation of invariant chemical composition. Since T changes far more slowly with distance in the Earth than in experiments, an isothermal bulk modulus represents the mantle values. The present paper assumes that changes in mineralogy of the LM are secondary, i.e., that the main changes in its seismic radial profile are from P and T, which permits use of Eq. (5). General behavior of bulk moduli for dense materials from both compression and elasticity studies supports this contention.

It is most fortunate that the derivatives are simply described:

$$\begin{array}{}\text{(6)}& \partial B/\partial P=B^{\prime }+B^{\prime \prime }P\text{and}\partial B/\partial T=\text{constant},\end{array}$$

where for dense and hard materials compatible with the LM, the constant B^′ = $\partial B/\partial P$ is commonly near 4 and $B^{\prime \prime }={\partial }^{\mathrm{2}}B/\partial P^{\mathrm{2}}$ is negative and sufficiently small to require very high pressure for its resolution (e.g., Sinogeikin and Bass, 1999; Zha et al., 2000). Results center on $B^{\prime }=\mathrm{4}$ because this value corresponds to a harmonic interatomic potential (e.g., Hofmeister, 1993) and anharmonicity links to T, not P, changes (e.g., Wallace, 1972). Most measurements provide $\partial B/\partial T$ as a constant. Although second-order T derivatives exist, these are small (if even resolvable) for hard oxide and silicate minerals, as shown in compilations and more recent work (e.g., Aizawa et al., 2004).

1.4 Synopsis of our novel, analytical inverse approach, and organization of the report

Section 2 shows how to extract the LM temperature gradient ($\partial T/\partial z$) from pressure and depth derivatives of radial seismic reference models by using Eqs. (5) and (6), in an inverse approach. Fitting is not used, which describes familiar, forward modeling. The extraction uses generic values for $\partial B/\partial T$ and ${\partial }^{\mathrm{2}}B/\partial P^{\mathrm{2}}$, which describes dense phases, including the rock salt and perovskite type structures thought to occur in the deep mantle, to explore the possible range of $\partial T/\partial z$ and its depth dependence. Thus, our thermal model is independent of what phases with what compositions might exist in the LM. The method is new, so details are provided.

Section 3 sets upper and lower bounds on heat transport properties for the LM based on verifiably accurate methods. We derive a simple formula for $\partial \ln \left(\mathit{\kappa }\right)/\partial P$ from Fourier's heat equation, which confirms that the result of Hofmeister (2021) is an identity. The sole parameter in the identity (other than B_T and P) appears be a constant, as suggested in Fig. 2 and independently by previous work (e.g., Chopelas and Boehler, 1992). The resulting bounds on κ(T,z) for the LM, combined with $\partial T/\partial z$ derived from PREM (Sect. 2), provide flux and power across the LM, without assuming its mineralogy.

Section 4 constructs geotherms using a reference temperature at a shallow level that avoids melting of a peridotite composition anywhere in the LM. These geotherms are independent of κ. Then, we ascertain which of our geotherms, fluxes, and powers are compatible with additional constraints, such as phase equilibria and latent heat of melting.

Our inverse model, which is based on radial seismic changes and high T and P behavior common to dense phases, indicates that the LM has a heat source which is warming the outer core, while causing the inner core to melt. Possible heat sources are discussed in Sect. 5, along with implications of our results.

2 Extraction of geothermal gradients from radial seismological models in an inverse approach

2.1 Features of PREM

Seismic reference models represent Earth's average interior; i.e., they are radial. Aspherical images of the Earth's internal structure are represented as perturbations to a reference mode (e.g., Ritzwoller and Lavely, 1995).

Reference model results are displayed as the fairly smooth functions of velocities, density, and pressure as a function of depth (z) or radius (s), or similarly as plots of bulk and shear moduli, the quantities of which are also fairly smooth, being derived from ρ and the two velocities. Seismic discontinuities are present as kinks, most of which are small in these typical representations. In contrast, large jumps dominate plots of derivatives of variables vs. depth (Fig. 4). Hence, this paper makes use of the derivatives.

Figure 4 Derivatives from the seismologic model PREM (Dziewonski and Anderson, 1981) as a function of depth. Variables as tabulated by Anderson (2007), from which we calculated the pressure derivatives of the bulk (black) and shear moduli (grey) as well as the depth derivative of pressure (dots).

Taking derivatives accentuates differences, as this mathematical operation is the converse of integration, which averages and smooths. The pattern exhibited by velocity derivatives (not shown) is similar to moduli derivatives, whereas the density derivative (not shown) is relatively smooth, more like the pressure derivative, and so the depiction of Fig. 4 is inherent to PREM.

Smooth and continuous $\partial B/\partial P$ describes the lower mantle but only between depths of 871 and 2741 km (Fig. 4). Importantly, other reference models such as Ak135 differ negligibly from PREM velocities in this restricted region (see Fig. 12 in Kennett et al., 1995). Our approach (below) applies to continuous functions only, so PREM suffices to represent all reference models of this volumetrically immense region. However, D^′′, where seismic models differ, cannot be quantitatively analyzed. We can only extrapolate into this tiny shell. For brevity, “lower mantle” or “LM” refers to its central region from 871 to 2741 km only, unless specified otherwise. Extrapolation of results for the LM into the underlying D^′′ layer and up to 671 km, which defines the transition zone (TZ), is discussed in Sects. 4 and 5.

2.2 Separation of variables

The geothermal gradient is defined by

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial T}{\partial z}}={\displaystyle \frac{\partial T}{\partial P}}{\displaystyle \frac{\partial P}{\partial z}}={\displaystyle \frac{\partial P}{\partial z}}\left[{\displaystyle \frac{\partial T}{\partial B}}{\displaystyle \frac{\partial B}{\partial P}}\right]={\displaystyle \frac{\partial P}{\partial z}}{\displaystyle \frac{\partial B}{\partial P}}{\left({\displaystyle \frac{\partial B}{\partial T}}\right)}^{-\mathrm{1}}.\end{array}$$

PREM provides the input quantity $\partial P/\partial z$ (Fig. 4). PREM values for B and $\partial B/\partial P$ as a function of depth are affected by both compression and heating of the minerals inside the LM. Our goal is to utilize available data to distinguish the effects of P and T on B from PREM. This is possible through separation of variables and decades of data acquisition.

Equation (6) provides the input for $\partial B/\partial P$. Because an inverse approach is being used, we consider values compatible with many dense phases, oxides, and silicates. The changes (differences) below z = 871 km are of interest, so the value of $\partial B/\partial P$ at 871 km serves as a reference point. This approach links B^′ and B^′′, as shown graphically in Fig. 5. Hence, obtaining thermal gradients from PREM via Eq. (7) requires some estimates of $\partial B/\partial T$ and ${\partial }^{\mathrm{2}}B/\partial P^{\mathrm{2}}$ but not of $\partial B/\partial P$. This is understood by considering two end-member cases.

The case B^′′ = 0 for the lower mantle sets an upper limit since B^′′ is not positive. If this commonly used limit (e.g., Knittle, 1995) applies to the LM, then $\partial B/\partial P$ with depth solely results from compression and thus is invariant (horizontal dashed–dotted line in Fig. 5). Consequently, rising temperature causes the linear decrease in PREM $\partial B/\partial P$ with z immediately below 871 km under separation of variables. For the second case, we consider B^′′ = −0.015 GPa⁻¹, which reproduces the decrease in PREM $\partial B/\partial P$ with z below 871 km. With this match, changes in B of PREM solely result from compression; i.e., T is constant for z slightly below 871 km. Hence, B^′′ between 0 and −0.015 GPa⁻¹ depicts a LM that is both compressing and warming below 871 km, whereas B^′′ > −0.015 GPa⁻¹ depicts a cooling and compressing LM just below 871 km. This case is not shown because T in the outermost layers of the Earth increases with depth, and so the top of the LM should behave likewise.

Once z reaches 1300 km, PREM $\partial B/\partial P$ curves have become rather flat, but as z increases further, deeper than ∼ 2200 km, PREM $\partial B/\partial P$ curves take on a positive slope with z, which is linear just below 2741 km. The broad minimum near 2000 km suggests that a maximum temperature may exist in the LM, where its manifestation depends on mantle values of $\partial B/\partial T$ and ${\partial }^{\mathrm{2}}B/\partial P^{\mathrm{2}}$ (discussed below). The positive slope in PREM $\partial B/\partial P$ curves at great depths, assuming that Eq. (6) represents the LM, shows that its deepest extent is shedding its heat downwards while compressing, discussed further in Sect. 4.

The above findings are general. Here we emphasize that following the discovery of heat generation by radionuclides, it was recognized that Earth may be heating up rather than undergoing progressive cooling. This question was considered in some remarkable papers (e,g., MacDonald, 1959). Notably, Takeuechi et al. (1967) devoted an entire chapter of their book to this subject.

Importantly, temperature differences from the starting point at 871 km are germane. Consequently, the input for $\partial B/\partial P$ in Eq. (7) is the PREM curve less a line for the P response of B, i.e., Eq. (6). Figure 5 shows the two end-member cases, discussed above, and one intermediate case. These three examples show that B^′ is controlled by the choice of B^′′, in order to match the starting point at 871 km. The range of values of 3.8 to 4.3 is compatible with mineral data (e.g., Knittle, 1995), with 4 being the harmonic value (Hofmeister, 1993).

Figure 5 Construction of possible curves for the depth response of the P component of mantle $\partial B/\partial P$ at 871 km and deeper. The edges of the transition zone and inner core define the x axis. Discontinuities shallower than 871 km and in D^′′ are not addressed in our approach. If ${\partial }^{\mathrm{2}}B/\partial P^{\mathrm{2}}$ is small, near 0, then PREM indicates that warming first occurs as z increases but at greater depths temperatures decrease. This behavior describes the LM and outer core.

2.3 Constraints on $\partial B/\partial T$and ${\partial }^{\mathrm{2}}B/\partial P^{\mathrm{2}}$ from elasticity and volumetric measurements of diverse phases

Regarding the thermal input parameter, $\partial B/\partial T$ equals −0.023 GPa K⁻¹ from the average of four experiments on MgO (compiled by Bass, 1995) and is only slightly larger in magnitude, −0.029 GPa K⁻¹, for MgSiO₃ with the perovskite structure (Aizawa et al., 2004), now known as bridgmanite. A large range of values is possible, per Aizawa et al. (2004) and work cited therein. More importantly, uncertainties are high while similar values are observed for corundum, spinel, five compositions of olivine, orthopyroxene, zircon, five garnets, seven other oxides, and four metals, whereas framework silicates, diamond, and alkali halides have $\partial B/\partial T$ near −0.01 GPa K⁻¹ (Bass, 1995; Anderson and Isaak, 1995).

An input value of −0.026 GPa K⁻¹ is used, the average of which is within reported experimental uncertainty of measurements of LM candidate minerals and moreover describes dense silicates and oxides in general. Because $\partial T/\partial z$ is inversely proportional to $\partial B/\partial T$, the effect of varying this parameter on the results is easily ascertained. In contrast, a non-linear response is associated with B^′′, since a difference with PREM is involved and PREM curves are non-linear with depth (Figs. 4 and 5). Hence, possible B^′′ values are the focus.

From the compilation of Bass (1995), B^′′ is positive for silica glass and near 0 for MgAl₂O₄ spinel, which is also disordered. Otherwise, B^′′ for oxides and silicates ranges from −0.03 to −1.6 GPa⁻¹. The largest magnitude depicts orthopyroxene, which has unusually high B^′ as well. For bridgmanite, B^′′ was not resolved even at compression to 155 GPa (Dorfman et al., 2013), consistent with incompressibility of high-pressure, dense phases. Aluminum being present makes no difference (Zhu et al., 2020). However, ubiquitous use of the Birch–Murnaghan equation of state, which involves trade-offs between B and B^′ (see, e.g., Knittle, 1995), prevents resolution of B^′′ for these stiff structures. Polynomial fits are needed to ascertain B^′′.

Many studies exist of MgO, but some ambiguity exists because elasticity data are fit to a polynomial in pressure, whereas volumetric data are analyzed using an equation of state, which is a general formulation for the P and T dependence of volume (V). The second-order polynomial coefficient for B(P) is $B^{\prime \prime }/\mathrm{2}$. To make sure this convention (i.e., Eq. 6) was used, original studies were consulted. Elasticity measurements of MgO by Sinogeikin and Bass (1999) provide B^′′ = −0.04 ± 0.02 GPa⁻¹. The X-ray diffraction study of Yoneda (1990) is consistent with the range of 0 to −0.029 GPa⁻¹. Zha et al. (2000) reached the highest pressures and found that the null value reasonably represents elasticity and volumetric data combined. For the dense LM, B^′′ is small. Again, B^′ = 4 and negligible B^′′ describes a harmonic solid (e.g., Hofmeister, 1993). We consider 0 to −0.02 GPa⁻¹, mostly in steps of 0.0025 GPa⁻¹, to calculate $\partial T/\partial z$ from Eq. (7).

2.4 Thermal gradients from 871 to 2741 km

Thermal gradients are shown in Fig. 6, with an example of a fit to B^′′ = −0.005 GPa⁻¹. The sign convention used here is based on flow from a central source (s = 0) moving outwards. All curves are well represented by third-order polynomials (Table 1) with similarly high residuals. Fits can be scaled to address variations from an input value of −0.026 GPa K⁻¹. Constrains on $\partial T/\partial z$ are covered in Sect. 4.

Figure 6 Dependence of thermal gradients from PREM (black symbols) on depth and the choice of ${\partial }^{\mathrm{2}}B/\partial P^{\mathrm{2}}$ for LM material. Fits are in color, where the inset lists coefficients for B^′′ = −0.005 GPa⁻¹, which has a typical residual. For positive gradients, T increases inwards.

Table 1 Polynomial curve fits to the thermal gradient in kelvin per kilometer (K km⁻¹), with parameters for the geotherm and power curves.

The polynomial for $\partial T/\partial z$ is $T_{\mathrm{0}}+T_{\mathrm{1}}z+T_{\mathrm{2}}{z}^{\mathrm{2}}+T_{\mathrm{3}}{z}^{\mathrm{3}}$. n/a: not applicable because large $\left|B^{\prime \prime }\right|$ leads to unsupportable melting that occurs near 670 km. n.d.: not determined. ^a These values satisfy the criteria that the LM is not melted and that the outer core melts above 2750 K as measured for the Fe–S system (see text). ^b Range at 2741 km from estimates of low and high κ. Power at the core mantle boundary (CMB) is lower, by about 1 TW. Positive sign indicates heat flow from the core into the LM, discussed below. ^c From the broad maximum in power which is nearly the same for low and high κ. ^d Most likely to represent the LM; see text.

3 Lower mantle transport properties from theory and experiment

3.1 Dependence of thermal diffusivity on temperature

Three-parameter fits describe measurements of thermal diffusivity for diverse solids above room temperature up nearly to melting that are neither affected by physical contact losses nor by spurious radiative transfer gains:

$$\begin{array}{}\text{(8)}& D\left(T\right)=F{T}^{-G}+HT\text{or}D\left(T\right)=F^{\ast }{\left({\displaystyle \frac{\mathrm{298}}{T}}\right)}^G+HT.\end{array}$$

The fitting coefficient G is near unity and H is small (Hofmeister et al., 2014). When H = 0, the parameter F* =F(298)^G on the right-hand side equals D at 298 K. The general applicability of these formulae has been established by additional measurements now encompassing over 200 substances (Merriman et al., 2018; Hofmeister, 2019). Equation (8) represents sample thickness L > 1 mm, i.e., bulk samples, and high temperatures and thus is appropriate to the mantle.

Examples of D(T) for dense phases are shown in Fig. 3. Generally, H is quite small, ∼ 0.0002 mm² s⁻¹ K⁻¹, but is essential to represent high-temperature behavior (T > 1200 K) of structures involving unit cell formulae more complex than Al₂O₃. For simple materials such as MgO and alkali halides occupying the cubic B1 and B2 structures, H = 0 within uncertainty. However, Si with the diamond structure has non-negligible H when impurities are present, whereas graphite, which has a more complicated anisotropic structure and is generally impure, has a substantial H term (Hofmeister, 2019, chap. 7). LFA data on glasses show that large H is commonly associated with high Fe cation content (e.g., Sehlke et al., 2020). These findings point to absorption bands above ∼ 1000 cm⁻¹ and into the visible region being associated with the HT term. To provide H when LFA data on (Mg_xFe_1−x)O are not available, we consider corundum, rather than MgO, to represent an oxide-rich lower mantle. For a silicate LM, systematic behavior of orthorhombic and cubic perovskites with various chemical compositions (Hofmeister, 2010a) is considered. The average D from the three orientations of NdGaO₃ is used to compute the T dependence of D for orthorhombic perovskite, since this agrees with D near 298 K for unoriented MgSiO₃ from Osako and Ito (1991). A periodic technique was used, and their sample was polycrystalline. Contact losses are more important than ballistic gains, because the latter is reduced by physical scattering between grains. Hence, D(T) for NdGaO₃ in Fig. 3a represents a minimum for a complex silicate phase in the LM.

Differences among dense silicates at high temperature are not large: this is the basis of D = 1 mm² s⁻¹ being commonly used in geophysical models. Near independence of D from T for complex solids at high T considerably simplifies calculations (below).

3.2 Dependence of thermal conductivity on temperature

Figure 3b shows examples of κ(T) calculated from Eq. (3). For MgO and Al₂O₃, c_P and ρ as a function of T are well constrained even at high T (e.g., Ditmars et al., 1982; Fiquet et al., 1997; Chase, 1998). Heat capacity data on MgSiO₃ perovskite (Akaogi and Ito, 1993) are limited to near-ambient T due to back-conversion problems. Using a model to extrapolate c_P is unnecessary, due to uncertainties in D(T) for silicate perovskite. Since D varies more strongly with T compared to c_P and ρ, which moreover respond in opposite directions during heating, ambient values are used to estimate κ(T) for a silicate perovskite mantle. To ascertain the uncertainty in this approach, analogous computations were made for κ(T) of periclase and corundum. Accurate and estimated trends for MgO differ little (Fig. 3b), but the Al₂O₃ estimate is significantly steeper with T than using exact values in Eq. (3). Corundum D(T) is very flat and is better fit to a polynomial than to Eq. (8), so the discrepancy is likely due to extrapolation beyond the temperatures actually measured. Another factor is that the T variations of both D and κ with T are weak when the corresponding ambient values are low, as is indicated in Eq. (9) and evident in Fig. 3. Thermal conductivity of other silicates depends similarly on T as our estimate for perovskite (for examples, see Hofmeister et al., 2014, and references therein).

Therefore, we estimate high-T mantle thermal conductivity in terms of constant, limiting values. For an insulating silicate LM, κ is taken as 2.7 W m⁻¹ K⁻¹, whereas for a thermally conductive oxide LM, κ is taken as 7 W m⁻¹ K⁻¹. The generic value of D used in geophysical models corresponds to ∼ 3.5 W m⁻¹ K⁻¹.

3.3 Dependence of thermal conductivity on pressure

Many formulae have been proposed for P derivatives of transport properties, based on dimensional analysis. An exact thermodynamic relationship,

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }}}{\left.{\displaystyle \frac{\partial \mathit{\kappa }}{\partial P}}\right|}_T& =-{\displaystyle \frac{\mathrm{1}}{V^{\mathrm{1}/\mathrm{3}}}}{\displaystyle \frac{\partial V^{\mathrm{1}/\mathrm{3}}}{\partial P}}+\mathit{\alpha }{\displaystyle \frac{\partial {\mathit{\alpha }}^{-\mathrm{1}}}{\partial P}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{3}{B}_T}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}{\displaystyle \frac{\partial \mathit{\alpha }}{\partial P}}={\displaystyle \frac{\mathrm{1}}{B_T}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}+{\mathit{\delta }}_T\right)\end{array},\end{array}$$

was derived and confirmed using reliable available data on 20 different homogeneous solids at pressures up to 2 GPa (Hofmeister, 2021). Equation (9) excludes a typographic error in the earlier report. The physical properties, other than κ, are part of the equation of state (EOS). Thermal expansivity is defined by

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\mathit{\alpha }}_P=& {\displaystyle \frac{\mathrm{1}}{V}}{\left.{\displaystyle \frac{\mathrm{\Delta }V}{\mathrm{\Delta }T}}\right|}_P={\displaystyle \frac{\mathrm{1}}{V}}{\left.{\displaystyle \frac{\partial V}{\partial T}}\right|}_P\text{ or}\\ & {\mathit{\alpha }}_{L,P}={\displaystyle \frac{\mathrm{1}}{L}}{\left.{\displaystyle \frac{\partial L}{\partial T}}\right|}_P={\displaystyle \frac{{\mathit{\alpha }}_P}{\mathrm{3}}},\text{if isotropic}.\end{array}\end{array}$$

Compressibility ($=\mathrm{1}/B_T$, the bulk modulus) is defined by

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\mathit{\beta }}_T=& -{\displaystyle \frac{\mathrm{1}}{V}}{\left.{\displaystyle \frac{\partial V}{\partial P}}\right|}_T={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\left.{\displaystyle \frac{\partial \mathit{\rho }}{\partial P}}\right|}_T\text{ or}\\ & {\mathit{\beta }}_{L,T}={\displaystyle \frac{\mathrm{1}}{L}}{\left.{\displaystyle \frac{\partial L}{\partial P}}\right|}_T={\displaystyle \frac{{\mathit{\beta }}_T}{\mathrm{3}}},\text{if isotropic}.\end{array}\end{array}$$

The dimensionless Anderson–Grüneisen parameter (δ_T) describes the opposing effects of T and P on V,

$$\begin{array}{}\text{(12)}& {\mathit{\delta }}_T=-{\displaystyle \frac{B_T}{{\mathit{\alpha }}_P}}{\left.{\displaystyle \frac{\partial \mathit{\alpha }}{\partial P}}\right|}_T=-{\left.{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_P{B}_T}}{\displaystyle \frac{\partial B_T}{\partial T}}\right|}_P,\end{array}$$

while clearly showing that the efficiency of expansion and compression for any given solid is related. Hence, the far-right-hand side (RHS) of Eq. (9) describes how temperature components of the EOS, not just pressure components, regulate changes in heat conduction during compression.

3.3.1 Derivation of $\partial \ln \left(\mathit{\kappa }\right)/\partial P$ from Fourier's equation

Due to the importance of compression to mantle heat transfer, we explain why Eq. (9) is exact. The original derivation considered diffusion of thermal radiation. A simpler approach is covered here.

Since ℑ is independent of pressure in experiments, taking the P derivative of its definition (Eq. 4) suffices to relate $\partial \ln \left(\mathit{\kappa }\right)/\partial P$ to EOS parameters. The algebra is simple and not specified here. However, one must recognize that a negative temperature gradient ($\partial T/\partial z{|}_P$) is associated with positive signs for ℑ and κ in Eq. (4). But Eq. (10) defining thermal expansivity is a scalar quantity, being based on volume, which has no direction. A positive sign for α, typical of most materials, requires the thermal gradient and its inverse, $\partial z/\partial T{|}_P$, to be positive in an isotropic solid. In contrast, heat flow has a well-defined direction from some origin and so is a vector quantity. Maintaining consistent signs for both κ and α during algebraic manipulations after differentiation of Fourier's equation leads to Eq. (9).

3.3.2 Experimental validation

The hot wire/hot strip and Angstrom techniques accurately measure thermal transport in metallic samples at low P (< 2 GPa) and low T (< 1000 K) since metal–metal contact losses are low and ballistic radiative transfer gains are negligible (e.g., Andersson and Bäckström, 1986; Jacobsson and Sundqvist, 1988). Moreover, use of standard L near millimeters in piston–cylinder and multi-anvil apparatuses permits direct comparison of results on diverse materials. Reliable data for insulators at high P over millimeter length scales have been obtained near 298 K using the hot strip/hot wire or Angstrom methods (e.g., Andersson, 1985; Osako et al., 2004). Samples are single crystals, glasses, and disks of fine-grained soft powder (alkali halides) that were compacted prior to study. Unlike metals, systematic errors exist due to interface thermal resistance and ballistic transfer, but taking a logarithmic pressure derivative minimizes these problems. Low-pressure transport property measurements and EOS results of over 20 solids, mainly from elasticity data, are summarized in Table 3 of Hofmeister (2021). The response of thermal conductivity to compression (Fig. 2) points to δ = 7 at ambient conditions representing the average for silicates, oxides, metals, alloys, and alkali halides.

Further verification is provided by a well-studied material with a special, negative sign of thermal expansivity. A negative sign for pressure response for thermal conductivity is expected and was indeed observed for silica glass by Andersson and Dzhavadov (1992) and Katsura (1993). Consistency with Eqs. (9) to (12) is demonstrated, since fused silica also has uncommonly positive $\partial B/\partial T$ (Spinner, 1956).

3.3.3 Previous EOS evidence for nearly constant δ_T

Chopelas and Boehler (1992) constrained mantle values of δ from 5 to 6 for metals, oxides, and alkali halides by assessing thermal expansivity at high pressure and temperature. Anderson et al. (1992) argued for an ambient value δ₀ = 6.5 and a weak volume dependence. The focus of these studies, along with Helffrick (2017), who further modified the V dependence, is the effect of compression on α.

Larger δ₀ = 7 was obtained from Fig. 2, which compares measurements of ∂κ/∂T to 1/B for the same types of solids. Hofmeister (2021) calculated δ₀ from temperature derivatives of bulk moduli, which are more accurate than P (or V) derivatives of α for several reasons. (1) Elasticity measurements determine $\partial B/\partial T$ as a first derivative and so are more accurate than x-ray diffraction (XRD) studies which determine $\partial B/\partial T$ as well as $\partial \mathit{\alpha }/\partial P$ as second derivatives. (2) Linear dependence of B on T exists over a wide range of temperatures (e.g., Anderson and Isaak, 1995), which simplifies establishing this parameter and reduces its uncertainty. (3) Because B is large in magnitude whereas α is small, derivatives of B are easier to determine accurately. Average δ₀ = 7 obtained in Fig. 2 better agrees with the EOS study of Anderson et al. (1992) because they included elasticity data in their assessment. (4) Data on thermal expansivity at pressure from XRD methods give α as an average over the temperature ranges explored, which makes this approach to P derivatives of α very uncertain. High uncertainty in α, let alone its P derivative, is evident in the compilation of Fei (1995).

Pressure-independent δ is consistent with separation of variables describing the bulk modulus, i.e., if

$$\begin{array}{}\text{(13)}& \begin{array}{rl}B=& B_o\left(\mathrm{1}+aT+b{T}^{\mathrm{2}}\mathrm{\dots }\right)\\ & \left(\mathrm{1}+B^{\prime }P/B_o+\mathrm{0.5}{B}^{\prime \prime }{P}^{\mathrm{2}}/B_o\mathrm{\dots }\right).\end{array}\end{array}$$

Then from the RHS of Eq. (12) follows

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\left.{\displaystyle \frac{\partial {\mathit{\delta }}_T}{\partial P}}\right|}_T& ={\left.{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_P^{\mathrm{2}}{B}_T}}{\left.{\displaystyle \frac{\partial \mathit{\alpha }}{\partial P}}\right|}_T{\displaystyle \frac{\partial B_T}{\partial T}}\right|}_P+{\left.{\displaystyle \frac{B_{\mathrm{0}}^{\prime }}{{\mathit{\alpha }}_P{B}_T^{\mathrm{2}}}}{\displaystyle \frac{\partial B_T}{\partial T}}\right|}_P\\ & -{\left.{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_P{B}_T}}{\displaystyle \frac{\partial \left(B_{\mathrm{0}}^{\prime }\right)}{\partial T}}\right|}_P={\displaystyle \frac{{\mathit{\delta }}_T}{B_T}}{\mathit{\delta }}_T-{\displaystyle \frac{B_{\mathrm{0}}^{\prime }}{B_T}}{\mathit{\delta }}_T\\ & +\mathrm{0}={\displaystyle \frac{{\mathit{\delta }}_T}{B_T}}\left({\mathit{\delta }}_T-B_{\mathrm{0}}^{\prime }\right).\end{array}\end{array}$$

Equation (14) equals 0 when ${\mathit{\delta }}_{\mathrm{0}}=B^{\prime }$ and is small otherwise. Equation (14) suggests that δ is a constant on the order of 4, the value of which for B^′ constitutes the simple Murnaghan EOS.

Similarly exploring the middle term of Eq. (12) indicates that δ also weakly depends on temperature when thermal expansivity is described by separation of variables. Constant δ is thus a reasonable first-order approximation. Exact evaluation is difficult due to the generally small sizes of all derivatives, limitations of a polynomial representation, and assumptions underlying the forms for the EOS. Compressible alkali halides are very important for this endeavor but are hydroscopic, easily deformed, and transparent in the infrared as single crystals.

3.3.4 Lower and upper bounds for thermal conductivity in the LM

The pressure (or depth) dependence of κ in the LM is much weaker for perovskites than oxides (Fig. 7), assuming δ₀ = 7. Corundum and perovskite have similar bulk moduli, so their different curves in Fig. 7 relate to relative efficiency of heat transfer at high T only.

Figure 7 Thermal conductivity in the lower mantle for oxide and silicate perovskite compositions assuming temperature independence and using Eq. (9) with δ = 7. Solid curves are computed for B^′′ = 0. Symbols utilize larger magnitude B^′′ than likely exists, given the results from PREM, but even this makes little difference.

Pure MgO would have lower κ than Al₂O₃ with depth due to T increasing beyond 1500 K (Fig. 3b) but would have higher κ with depth due to P increasing. For this reason, corundum, which has the same κ as periclase at 1500 K, is used to represent an oxide lower mantle (Fig. 7).

4 Results

Calculations of geotherms, flux, and power from 871 to 2471 km are presented here, which are extrapolated through D^′′. Comparison is made with possible melting temperatures to eliminate cases unlikely to represent Earth's interior. The limiting case of B^′′ = 0 is unexpected but is useful for comparisons. Importantly, our geotherms do not utilize data on thermal conductivity and thus are essentially independent of mineralogy. In contrast, flux and power are independent of the reference temperature, but use κ and thus are affected by LM mineralogy. However, as κ varies little at high T (Fig. 3), only the proportion of complex silicates to simple oxides matters.

4.1 Calculation of temperatures across the lower mantle from PREM and a reference point

Geotherms are calculated by integrating the thermal gradients (Fig. 6) which were obtained from PREM, from 871 km downwards, using various B^′′ values (Table 1). Results (Fig. 8) consider three values for the 871 km reference temperature (T_ref) representing the top of the lower mantle. The minimum T_ref of 1500 K is based on temperatures for basalt extruding at the surface (e.g., Falloon et al., 2008), whereas the maximum T_ref of 2500 K is based on dry melting of peridotite at 670 km (Zhang and Herzberg, 1994). The intermediate (T_ref = 2000 K) corresponds to Takahashi's (1986) melting curve of peridotite at ∼ 410 km, which probably involved tiny amounts of moisture, based on sensitivity of melting to hydration.

All geotherms (Fig. 8) are flat near 871 km, due to PREM derivatives being linear for the shallowest LM (Fig. 6). Thus, the inverse method suggests that temperatures change little from 871 km up to the shallower depth of the 671 km seismic discontinuity. Chemical composition could be changing near 670 km, as descending slabs are resolved from earthquakes down to these depths but not below (summary figures are in Hofmeister, 2020, chap. 7, and in Hofmeister et al., 2022).

Regarding the base of the LM, the geotherms smoothly decrease with depth. Extrapolation of LM results from 2741 to 2891 km presumes similar bulk moduli derivatives for D^′′ and the LM. Although this seems unlikely, since reaction with core material is possible, discontinuities at these depths preclude robust analysis (Sect. 2), and modeled velocities in this region are less certain than in the LM (e.g., Kennett et al., 1995). But, as discussed in Sect. 2.3, as evidenced by compilations and recent work, derivatives of B vary little with structure, composition, and bond type.

Melting curves of Fe–S and Fe–Ni–S systems (Chudinovskikh and Boehler, 2007; Morard et al., 2011; Mori et al., 2017) set a minimum of ∼ 2750 K at 2891 km, whereas peridotite melting (Fiquet et al., 2010) sets a maximum near 4100 K (Fig. 8). Only for the combination of the highest T_ref and the limiting case of B^′′ = 0 is peridotite melting reached. Many, but not all, geotherms exceed the Fe–S eutectic. For example, if T_ref = 2000 K, then $\left|B^{\prime \prime }\right|$ must be smaller than 0.0075 GPa⁻¹. Table 1 lists the range of reference temperatures consistent with the above-mentioned phase equilibria for each B^′′ value considered. Only for the hottest possible T_ref = 2500 K can B^′′ reach −0.01 GPa⁻¹. Small second derivatives are consistent with difficulty in resolving these in experiments on dense materials, unless extreme pressures are reached (see Zha et al., 2000) and polynomial fits are used (Sect. 2.3).

Figure 8 Families of geotherms for possible values of second P derivatives of LM bulk modulus and three different starting temperatures at 871 km. All curves use the T derivative of B for silicate perovskite from Aizawa et al. (2004). For the dry peridotite solidus, results of Zhang and Herzberg (1994) are merged with high P data or Fiquet (2010). The thick black line is the eutectic melting curve of the Fe–S system (e.g., Morard et al., 2017).

Irrespective of the temperature values, the shape of the geotherms require a thermal maximum inside the lower mantle. For large magnitudes of B^′′, a maximum is indicated roughly near 670 km or slightly shallower. For the smallest $\left|B^{\prime \prime }\right|$, the maximum T in the lower mantle is reached at its interface with the core. In all cases consistent with phase equilibria, maximum LM temperatures are reached for z below 2200 km. Thus, from analyzing PREM curves, LM temperatures climb inwards for most cases considered. A heat source located in the LM is supported by flux and power calculations, below.

Altering our input value for $\partial B/\partial T$ from −0.026 GPa K⁻¹ will either expand or contract the splayed patterns of Fig. 8, which rest on thermal gradients of Fig. 6 obtained from Eq. (7). Comparison of such revised curves with phase equilibria would then change the ranges of B^′′ and T_ref that avoid melting in the LM while allowing melting in the core (Table 1, RHS), resulting in quite similar shapes. Thus, geotherms from PREM are robust, given the subsidiary information on melting relations, whereas the specific input values are interdependent. The shapes of the geotherms are compatible with the process of heat diffusion, i.e., thermal conduction, even though the calculations (depicted in Table 1 and Figs. 6 and 8) did not incorporate thermal conductivity values.

4.2 Calculation of flux across the lower mantle from PREM and κ

Flux, defined by Eq. (4), describes spherical geometry for radially changing T. Temperature values are not needed to compute the amount of heat being moved inside and across the LM, when κ only weakly depends on T. The gradients of Fig. 6 and Table 1 lead to families of flux vs. depth curves which depend on B^′′ values for each of low and high κ. The sign convention used here portrays heat from the center of a sphere moving outwards.

All cases (Fig. 9) provide a broad peak for ℑ across the lower mantle. For comparison, surface flux values are larger, averaging ∼ 60 mW m⁻² for either crust (Veiera and Hamza, 2018). The height of the peak in ℑ increases as B^′′ approaches its null limit. Only for this limiting case (B^′′ = 0) is flux within D^′′ large and positive, the behavior of which signifies that heat emitted from the core contributes to the flux in D^′′ and at large z in the LM. For the next larger value of B^′′, flux near D^′′ is positive but near zero. Thus, the vast majority of our calculations point to a lower mantle source whose heat is being shed to its adjacent layers. Section 4.3 provides further discussion.

Note that the curvature and a maximum in ℑ are inherent to PREM (Figs. 4 to 6; Sect. 2.3). The increase in thermal conductivity with P serves to flatten these curves at great depths near the outer core and, importantly, means that the efficiency of heat conduction increases with depth. Consequently, heat from a source deep in the Earth is conveyed more readily inwards than outwards.

How much heat is carried depends on the material, i.e., on the high-T value at ambient pressure. Our model (Fig. 7) rests on heat transport values that are established via measurements. The thermal response at modest temperatures (< 2000 K) of any given material is largely controlled by infrared fundamentals and near-IR overtones, whose frequencies overlap with the associated blackbody curve (Hofmeister, 2019, chap. 11). We have not accounted for electronic transitions of Fe²⁺ augmenting heat transfer significantly above ∼ 1000 K, as observed in various glasses (Sehlke et al., 2020, and references therein), since the chemical composition of the LM is not known. The purpose of this paper is to ascertain Earth's thermal state with minimal assumptions and input parameters. As accuracy is not possible given the scant definitive information on LM mineralogy, such as samples, salient features not affected by the details are pursued here.

Nonetheless, enhancements in thermal conductivity with depth would raise the flux near 2471 km and straighten the curves, resulting in melting in the deepest lower mantle. We suggest that κ cannot be significantly larger than that considered in Fig. 7 or the boundary layer D^′′ would be significantly larger and mostly molten, which contradicts observations of shear waves in this region (cf. behavior of velocities in the molten outer core from PREM, Fig. 5, to those in the solid layers).

4.3 Calculation of power across the lower mantle from PREM and κ

In spherical coordinates, power is provided by

$$\begin{array}{}\text{(15)}& \mathrm{\wp }=\mathrm{4}\mathit{\pi }{s}^{\mathrm{2}}F.\end{array}$$

For all parameters explored, ℘ has a broad peak in the lower mantle (Fig. 10). The depth where the maximum ℘ occurs (Table 1) points to the location of a heat source. This source is in the lower mantle, unless $\left|B^{\prime \prime }\right|$ is larger than considered, which suggests a location near or in the TZ.

Power is supplied to the LM from the core only for the case of high κ with quite small $\left|B^{\prime \prime }\right|$, both of which are unexpected. For this unlikely situation to exist, the core would also need to have a heat source with a sufficiently large associated outward flux to overtake the flux inwards from the LM source. Here, we pursue a simple explanation, consistent with realistic input parameters, that a source in the LM supplies heat to the core region. In detail, half of the four curves require B^′′ = 0, which is unexpected, and a third curve points to negligible power from the core. The fourth case of an oxide (high κ) with low B^′′ = −0.0025 GPa⁻¹ is probably incompatible parameters. This inference is consistent with large $\left|B^{\prime \prime }\right|$ being associated with compressible solids like salts (e.g., Bass, 1995). Periclase, which has been viewed as a lower mantle phase, is much more compressible than silicates with perovskite-type structures, considered to dominate the LM (see Sect. 2.3). If the mantle is mostly composed of dense silicates, as is the current view, the LM is an inefficient heat transmitter and stays hot where the heat is produced.

Heat conveyed from the lower mantle to the core could cause melting. The latent heat of ∼ 450 kJ kg⁻¹ for iron melting at high P and T (Aitta, 2006) is close to latent heats for basalt and many other materials. A constant rate of melting over geologic time requires 5.8 TW to make the outer core. This value is most compatible with the case of low κ (silicate) and B^′′ = −0.010 GPa⁻¹. If the source is winding down with time, as is likely, then B^′′ should be more positive. The case of low κ (silicate) and B^′′ = −0.0075 GPa⁻¹ is compatible with a wide range of temperatures for the top of the LM and a sulfide-rich core, which addresses the expected sulfur concentration from meteoritic models (see tables in Lodders and Fegley, 1999). These parameters compose our best estimate of the thermal gradient, geotherm, flux, and power, while suggesting that heating occurs near z = 1900 km (Table 1).

Figure 9 Calculated flux for high oxide and low silicate thermal conductivity, as well as various second P derivatives of the bulk modulus, as labeled. The same patterns are used in Figs. 6 and 8. Arrows are extrapolations across D^′′, showing that most cases indicate a heat source in the LM. When flux is positive, the temperature is increasing inwards.

Figure 10 Power calculated from the flux curves in Fig. 9 and Eq. (15). Positive power is associated with temperature increasing inwards. Arrows indicate the direction of flow across the core–mantle boundary region and also into the TZ.

5 Implications and conclusions

Earth's largest zone by mass or volume is the lower mantle, which has only been accessed remotely, via seismic data acquisitions and processing. PREM, Ak135, and other reference models provide nearly identical velocities (Kennett et al., 1995). The models yield smoothly varying properties and their derivatives over most of the LM, the changes of which are taken in this report to represent combined effects of P and T varying with depth. However, a smooth variation in chemical composition is not precluded in utilizing separation of variables (Sect. 2). Gradual changes in chemical composition could be hidden in our choices for B^′′ if changes in mineralogy lead to a linear response of $\partial B/\partial P$. Specifically, mineralogy depending linearly on density and thus on P would lead to an equation equivalent to Eq. (6). Indeed, systematic dependence of velocity (and thus of B) on density exists and has received much attention (see, e.g., the seminal paper of Shankland, 1972). But because $\partial B/\partial P$ varies little among dense phases, as underscored by commonly assuming harmonic values of 4 in analyzing data (see Sect. 2.3), geotherms calculated from PREM via the inverse approach developed here are largely independent of mineralogical details. Values considered for B^′′ depict all three different thermal situations that are possible below 670 km: LM temperatures can increase with depth below the TZ, or are constant, or may decrease. Only one case was considered for the last, due to strongly decreasing T being incompatible with a molten outer core (Fig. 8).

The thermal state and gradients inward from the top of the LM indicate that a power source exists in the lower mantle. Locations of the peak in ℘ or flux are little affected by the high-T values of thermal conductivity considered in the inverse model: instead, the location of the power source is largely controlled by B^′′. As discussed in Sect. 2.3, B^′′ is not well constrained. This ambiguity effectively lumps compositional changes with those solely due to compression. Contrastingly, the heights of the peaks in Figs. 9 and 10 are affected by values for κ and for $\partial B/\partial T$. Neither variable depends strongly on mineralogy. Note that the large contrast (in κ) explored considers vastly different simple oxide vs. complex silicate compositions.

Notably, thermal models show that the scale length for cooling over geologic time is ∼ 1000 km (Criss and Hofmeister, 2016), so heat generated in the LM cannot reach the surface to any appreciable degree. A similar scale length is obtained from Eq. (2): the Earth cools slowly because it is large and spherical with a refractory outside. Flat geotherms near z = 1000 km with a source near 2000 km are compatible with the thermally insulating nature of rocks and our earlier cooling models. Criss and Hofmeister (2016) used a constant κ that is independent of pressure and so assumed a much more thermally insulating mantle. The flat and low κ case underestimates ℑ in the lower mantle. The consequence is a narrower thermal maxima with little heating of the core in the calculations of Criss and Hofmeister (2016). The broad geotherms inferred from PREM without knowing thermal conductivity are thus consistent with the strong dependence of thermal conductivity on pressure that is demonstrated by experiments (Fig. 2) and, importantly, is inherent to Fourier's definition of ℑ (Sect. 3.3).

In short, results (Figs. 6, and 8 to 10) obtained by analyzing the seismologic representation of Earth's interior are consistent with studies of the equation of state, phase equilibria, and thermal transport properties. Table 1 lists the most likely thermal gradient for the parameters explored. Although other values are possible, the mostly likely gradient is unlikely to change significantly, as likelihood was deduced from melting temperatures for LM and core candidate materials.

Heat is produced in the LM, but how? Commonly considered sources are discussed next. The paper concludes with a proposal.

Heat produced in rocks by radioactive decay has been the focus of many studies. Yet, it is well known that U, Th, and K are concentrated in the continents, leaving very little for the mantle, if meteorites represent the bulk Earth. Flux from the oceans suggests mantle production of <100 W km⁻³. Such a tiny source can only heat the interior if deeply buried, but for this case excessively high T results, due to the time evolution of the heat-generating isotopes (Criss and Hofmeister, 2016).

This long-standing problem in geochemistry led to considering primordial heat as another source. This hypothesis is based on gravitational contraction producing heat and rests on Kelvin's discounted proposal for the generation of starlight. Changes in gravitational potential produce motions as per elementary physics textbooks. Conversion of gravitational potential to spin and orbits quantitatively accounts for the high kinetic energy for Earth and sister planets today, as well as high spin observed for young stars (see Hofmeister and Criss, 2012). Some accretionary heating is expected in the final stages, but this source is a small fraction of the gravitational potential, non-renewable, shallow, and winding down. Likewise, core formation is not a source of heat: rather, the planet would need to be already melted for a homogeneously accreted object to sort, since self-compression itself provides a stable density stratification.

Heat is a by-product of motions, when accompanied by deformation, non-elastic in particular, or friction. Motions are produced by forces which on large, planetary scales are gravitational in origin. On this basis, and because the previously explored sources of radiogenic and primordial heat are inadequate to describe the workings of Earth (e.g., the hypothesis of mantle convection: see Bercovici, 2007) as well as seismologic detection of a molten core, which is hot, our recent efforts have focused on forces and motions. Hofmeister et al. (2022) argue that the location of the barycenter, where the immense solar pull and orbital centrifugal forces balance, differing from that of the geocenter results in imbalanced stresses and forces that are cyclical with periods of both 1 d and 1 month, with plate tectonics being a consequence. Cyclic stresses promote failure (Schijve, 2009). Spin is important, as the force field is axially symmetric, which explains orientation of the mid-ocean ridge fracture system.

The barycenter is a point in space that the Earth spins through. Its location relative to the geocenter is defined by masses of the Earth and Moon plus the lunar distance, which varies over the month. The depth range of ∼ 1450 to 2050 km, shown in Figs. 8 to 10, includes the position of the power source indicated by B^′′ more negative than −0.005 GPa⁻¹. Over the day, this point in space moves through the LM and rarely lies in the equatorial plane, due to Earth's tilted spin axis. Thus, our results for Earth's thermal state are compatible with cyclical stresses heating the LM. This region is strong and plastic (or elastic) rather than brittle, like the lithosphere, but both are underlain by fluid layers. Liquids flow under any stress, the lack of rigidity of which adds stress to the overlying layers. The amount of heat generated is small, ∼ 1 TW, and is conducted away from the source in both directions. Investigating this proposal further is beyond the scope of the present report.