Multidisciplinary Constraints on the Thermal‐Chemical Boundary Between Earth's Core and Mantle

Heat flux from the core to the mantle provides driving energy for mantle convection thus powering plate tectonics, and contributes a significant fraction of the geothermal heat budget. Indirect estimates of core‐mantle boundary heat flow are typically based on petrological evidence of mantle temperature, interpretations of temperatures indicated by seismic travel times, experimental measurements of mineral melting points, physical mantle convection models, or physical core convection models. However, previous estimates have not consistently integrated these lines of evidence. In this work, an interdisciplinary analysis is applied to co‐constrain core‐mantle boundary heat flow and test the thermal boundary layer (TBL) theory. The concurrence of TBL models, energy balance to support geomagnetism, seismology, and review of petrologic evidence for historic mantle temperatures supports QCMB ∼15 TW, with all except geomagnetism supporting as high as ∼20 TW. These values provide a tighter constraint on core heat flux relative to previous work. Our work describes the seismic properties consistent with a TBL, and supports a long‐lived basal mantle molten layer through much of Earth's history.

Q ad , may exceed Q CMB if the thermal conductivity is high enough, which would cause a thermally stratified layer to develop at the top of the core (e.g., Gubbins et al., 1982). The Q CMB must be large enough to sustain the present-day magnetic field, but a more stringent constraint may come from the requirement that thermal stratification not exceed several hundred kilometers (Christensen, 2018;Olson et al., 2017). Otherwise, the morphology of the magnetic field would deviate substantially from the observed field. Based on recent core energy budget models with the revised high thermal conductivity, Q ad is ∼15 TW. Using this Q ad , we calculate that Q CMB must be at least 10.5 TW to avoid an unreasonably thick thermal stratification at the top of the outer core (see Section 6 for details). From the perspective of the geodynamo, an upper limit on Q CMB is imposed by the requirement of a dipolar field, but this value is probably too large to be practically useful.
Numerical mantle convection models can provide estimates for Q CMB by combining the mantle temperature gradient with material properties of mantle silicates. A simpler approach relies on the thermal boundary layer (TBL) theory (Turner, 1979), which relates heat flow to the temperature drop between the CMB and the interior of the convecting mantle. While this approach is commonly used in parameterized convection models (e.g., P. Driscoll & Bercovici, 2014), it requires extensions to account for strongly temperature-dependent viscosity or chemical heterogeneity in the mantle. Additionally, the presence of lateral variations in seismic velocities may be too large to be purely thermal anomalies. Compositional heterogeneity may introduce density anomalies that alter mantle flow and cause deviations from the heat flow predicted in simple models. For example, the Large Low-Shear Velocity Provinces (LLSVPs; e.g., Garnero et al., 2016) which are suggested to be compositionally distinct from the ambient mantle, may influence the TBL by acting as insulators, locally decreasing the heat flux across the CMB (Olson et al., 2015). On the other hand, the TBL theory offers useful expressions for the thickness and residence time of boundary layers at the base of the mantle. Estimates of Q CMB from numerical models of mantle convection depend on the imposed temperature boundary conditions and choice of material properties. Recent models have achieved a large range of values of Q CMB between 5 and 15 TW (Nakagawa, 2015;Olson et al., 2015). We discuss the simplified, analytical TBL theory in Section 3, and use numerical methods applied to the lowermost mantle to test this theory in Section 4.
Seismology also offers a test of the TBL theory. A temperature increase within the TBL should affect the elastic properties of lower mantle materials, manifesting as a decreased seismic velocity just above the CMB. Given that the thermal anomaly resulting from the TBL would be maximum at the CMB and decrease with distance above the CMB, the velocity structure at the CMB will be most indicative of a TBL. The CMB structure will be best resolved by seismic methods with high sensitivity at that depth.
Measurable perturbations in the travel time, incoming direction, amplitude, and reflection behavior of body waves, as well as the frequency and splitting of Earth's normal modes indicate radial and lateral variations in elastic structure. There is evidence of velocity variations in the lower mantle on a range of scales: the LLSVPs mentioned above (Garnero et al., 2016), the D″ reflector (for a review see Cobden et al., 2015), and the Ultra Low Velocity Zones (ULVZs; Garnero & Helmberger, 1995). However, these are primarily local structures, and are not representative of the thermal state of the global CMB. Looking outside these regions, several models of the 1D global average velocity structure of the mantle show shallower seismic velocity gradients in the vicinity of the CMB (Dziewonski & Anderson, 1981;Kennett & Engdahl, 1991;Kennett et al., 1995). The gradient observed in the preliminary reference Earth model (PREM; Dziewonski & Anderson, 1981) has been interpreted as a 400-1000 K temperature increase close to the CMB, although lower mantle anisotropy casts uncertainty on the observed velocity gradient (Doornbos, 1983;Doornbos et al., 1986;Stacey & Loper, 1983). Lateral variation in the seismic velocity gradient in the lowermost mantle suggests geographically varying temperature gradients (Young & Lay, 1987). However, concurrent observation of a density increase with a velocity gradient decrease close to the CMB may be better interpreted as evidence of compositional stratification, given the small effect of temperature on density (Lay, 1989). We explore the possibility of observing the TBL seismically in Section 5.
As mentioned above, the outer core may also host a TBL. Convective processes, chemical exchange with the mantle, and light element transport in the outer core may lead to a stably stratified layer at the top of the outer core (B. A. Buffett & Seagle, 2010;. Observational evidence of such a layer comes from geomagnetic (B. Buffett, 2014) and geodetic (Braginsky, 1984) studies. Additionally, some seismic studies using the SmKS family of body waves have shown evidence of a layer of reduced seismic velocities in the uppermost 100-400 km of the outer core (Garnero & Helmberger, 1995;Garnero et al., 1993;Helffrich & Kaneshima, 2010;Kaneshima, 2018;Kaneshima & Helffrich, 2013;Kaneshima & Matsuzawa, 2015;Kohler & Tanimoto, 1992; The heat flux across the CMB is controlled by and influences many other Earth processes. Here, we leverage all of these approaches to obtain self-consistent estimates for the Q CMB and the thermal structure of the CMB boundary layer. We test the geochemical and geodynamic assumptions and validity of the current TBL model in both a compositionally homogeneous and heterogeneous mantle. We then make predictions for the observability of such a layer with various seismic methods. Finally, we calculate the implications of the TBL on convection in the outer core, and thus the longevity of the magnetic field through time.

Energy Balance and Estimates of CMB Temperature
Important constraints on the CMB heat flow (Q CMB ) can be inferred from energy balances for the core and mantle (see Figure 1), together with estimates of temperature in the vicinity of the CMB on either side of the boundary. In this section, we summarize constraints from petrology and mineral physics to provide updated anchor points for the temperatures in the core and mantle symbols used listed in Table 1.

Energy Budget for the Silicate Earth
The energy balance of the bulk silicate Earth (BSE), defined as the mantle plus the crust, may be written as where M is the mass of the BSE (4 × 10 24 kg), C P is the average heat capacity of the mantle (∼1,000 J/kgK) (Olson et al., 2015), and ̄∕ is the mass-averaged cooling rate of the mantle. H tot is the radiogenic heat production of the BSE. Q surf and Q CMB are the heat flow across the Earth's surface and CMB, respectively. Q surf is measured to be 46 ± 3 TW (Lay et al., 2008), while Q CMB is unknown. Q CMB is dictated by the temperature gradient across the TBL, ∇T, as well as the thermal conductivity of silicates at the CMB, k mantle , and may be expressed as, where S is an outward normal vector and q CMB is the heat flux per unit area. Quantification of Q CMB using the energy budget requires the knowledge of the mantle cooling rate and radiogenic heat production. The latter can be related to surficial heat flow by the Urey ratio, defined as radiogenic heat production over total heat output from the mantle, that is, Ur = (H tot − H cc )/(Q surf − H cc ), where H cc is the radiogenic heat produced in continental crust, ∼7.5 TW (Rudnick & Gao, 2003). Estimates of the Urey ratio span a wide range −0.08 ± 0.05, 0.3 ± 0.1, and 0.7 ± 0.1 based on cosmochemical, geochemical, and geodynamical arguments for the magnitude of H tot , respectively (Sramek et al., 2013). We explore the implications of these values in Section 3.2.
The other unknown in Equation 1 needed to estimate Q CMB is the cooling rate of the mass-averaged mantle temperature (̄∕ ) . The mass-averaged temperature is defined as Figure 1. Summary of the energy balance of the Earth. The energy budget of the mantle is discussed in Section 2.1 and the energy budget of the core is discussed in Section 6.1. A schematic of the core-mantle boundary (CMB) thermal boundary layer (TBL) is given in the blown up image. The dashed red curve is the theoretical temperature extrapolated from the geotherm to the CMB and the solid red curve represents the actual temperature profile, following the TBL model.
where T is the geotherm through the mantle and V is the volume of the BSE. We assume that the geotherm is approximated by an adiabat, anchored to the mantle potential temperature, T P , at the surface. Using this assumption in Equation 3 allows us to relate ̄ to T P . An example using the geotherm of Stixrude and Lithgow-Bertelloni (2011) and the PREM model (Dziewonski & Anderson, 1981) gives ̄= 1.26 . Other geotherm models (see Figure 2) give constants in the range of 1.25-1.30. One exception is the geotherm of Anderson (1982) because it includes a temperature jump at the base of the upper mantle. For subsequent calculations we adopt the geotherm of Stixrude and Lithgow-Bertelloni (2011) and let ̄∕ = 1.26 ∕ . Changes in the potential temperature can now be quantitatively related to the imbalance between heat loss due to conduction through the Earth's lithosphere, and heat gain by radioactive heating in the mantle and conduction into the base of the mantle.
Changes in the potential temperature and the cooling rate of the mantle can also be inferred from petrological observations. Primary mafic magmas are formed by partial melting of the mantle and without further modification by differentiation. The chemical composition of primary mafic magmas is a critical indicator of physical and chemical conditions (i.e., pressure, temperature, composition, fluid availability, and accompanying melt production) in the mantle and of melt evolution in the lithosphere. Primary magma MgO and FeO contents increase with T P , which provides a petrological record of the thermal state of the mantle from which they formed, the age of which can be measured independently with radiometric dating. Figure 3 summarizes the T P values inferred from petrological modeling of 1,167 Archean and Proterozoic basalts derived from the GEOROC database. We determine T P using PRIMELT (Herzberg & Asimow, 2015) based on the MgO content of the primary magma. The relation between MgO and T P depends on assumptions concerning the thermodynamic properties of decompression melting. Previous studies of T P obtained with MgO show that primary magmas with 9-11 wt% MgO required a T P in the range of 1600-1700 K at the present time (Herzberg et al., 2010;Ganne & Feng, 2017;Servali & Korenaga, 2018). Here, the present-day mantle potential temperature is assumed to be 1625 ± 50 K (Herzberg et al., 2010;Korenaga, 2008;star symbol in Figure 3), to produce primary basaltic magmas having 10%-13% MgO. Of the initial 1,167 data points drawn from the database, we filtered the data to keep MgO >9 wt% and SiO 2 >45 wt%, resulting in a final data set of 54 non-arc basalts (blue circles and triangles in Figure 3), which we combine with 33 basalts reported by Herzberg and Gazel (2009) (yellow circles in Figure 3). The T P value for our filtered sample converges to the present-day ambient value, suggesting that these samples record the secular evolution of ambient mantle temperature. Most of the newly added points (Figure 3, blue circles) are located within the previously constrained region by Herzberg and Gazel (2009) (Figure 3, yellow circles), whereas several newly added samples show lower T P values following a different trend (blue triangles). These samples are flood basalts from the Mid-Continent Rift in North America (Paces & Bell, 1989), North Atlantic Craton (Hall et al., 1985), and Tanzania Craton (Manya & Maboko, 2003). They might represent thermal heterogeneities in the ambient mantle (e.g., the sublithospheric mantle; Ganne & Feng, 2017) but not the broader cooling trend. It is noticeable that the calculation can distinguish lavas that formed from ambient mantle from those that formed in mantle plumes: first, the data is filtered by MgO and SiO 2 content according to Herzberg et al. (2010) and Ganne and Feng (2017); second, successful primary magma solutions obtained by the program PRIMELT are internally consistent with primary magma compositions of fertile peridotite composition KR-4003 (Herzberg et al., 2010), which can represent the ambient mantle.
Figure 3 also shows predictions for T P versus time based on thermal evolution models of the mantle with a Urey ratio of 0.2 or 0.3, representing the uncertainty in the absolute abundance of long-lived radioactive elements (Korenaga, 2008), which are consistent with the petrological modeling. On the other hand, some T P values  Korenaga (2017), are consistent with the models of Ganne and Feng (2017) or G. F. Davies (2009), representing a lower bound in mantle cooling rate. Nevertheless, the cooling rate of the last 2 billion years is constrained to dT P / dt = 80-110 K/Ga based on both the broader petrological record (Herzberg and Gazel (2009); and our filtered data) and the mantle thermal evolution model (Korenaga, 2008).

Estimates of Temperature at the CMB
The temperature throughout the lower mantle and core is often approximated by assuming that the regions away from the boundaries are chemically uniform and thermally well mixed. This is a reasonable assumption in a vigorously convecting material, where heat is mainly transported adiabatically by fluid motion, and conductive and radiative heat transfer are negligible. In the mantle, the geotherm is an adiabatic temperature gradient with anchor points (Figure 2) based on the pressure and temperature of known mineralogical phase transitions. Meanwhile, the adiabatic temperature of the outermost core is anchored on the melting behavior of iron alloys at the conditions of the inner core boundary (ICB). In this subsection, we summarize the constraints on T CMB from the mantle and core sides of the boundary layer.

Mantle-Side Temperature Constraints
Uncertainties around composition, mineralogical phase changes, and material properties lead to a range of temperature profiles from Earth's surface to the CMB, summarized in Figure 2. All models assume an adiabatic temperature distribution through the mantle, with various deviations near the boundaries. An adiabatic temperature gradient can be described as where T is temperature, z is depth, g is gravitational acceleration, α is the coefficient of thermal expansion, and C P is the isobaric heat capacity of the aggregate material (e.g., Turcotte & Schubert, 2002). Gravitational acceleration can be estimated from the density profile of the Earth. At the high temperatures of the mantle, isobaric heat capacity does not change significantly with increasing pressure and temperature (Turcotte & Schubert, 2002). Thus, reliable constraints on the thermoelastic properties (e.g., α) of major mantle minerals, such as olivine, wadsleyite, ringwoodite, and bridgmanite are the basis of mantle adiabatic temperature estimation.
Using the equations of state of major mantle minerals (olivine, wadsleyite, ringwoodite, and bridgmanite), Katsura et al. (2010) re-evaluated the temperature at the 410 km discontinuity, and extrapolated it to both the deeper and shallower regions. Their adiabat yields a temperature of 2755 K at a depth of 2,800 km. Previously, Brown and Shankland (1981) proposed a slightly lower (by 70-110 K) temperature profile using a similar thermodynamic model to that of Katsura et al. (2010), but with the temperature at the 660 km discontinuity 50 K lower than the more recent study. Hence, the adiabatic temperature of 2449 K at 2,885 km depth derived by Brown and Shankland (1981) may serve as a lower bound. The temperature profile proposed by Stacey and Davis (2008) is in agreement with Katsura et al. (2010), but the temperature gradient in the lower mantle was assumed to be a constant 0.4 K/km, whereas the gradient shown in Katsura et al. (2010) was largely decreasing with depth based on the updated equations of state of mantle minerals. Similarly, with new experimental constraints, Ono (2008) estimated temperatures at depths of 660 and 2,700 km of 1850 and 2600 K, respectively, by comparing the depth of seismic discontinuities with pressures of post-spinel and post-perovskite transitions in pure MgSiO 3 . However, the effect of Fe concentration on phase transitions, and thus anchoring temperatures, was not considered. As an upper bound, Katsura et al. (2010) used an adiabatic temperature of 2937 K at 2,886 km depth, which Estimates of the geotherm through the mantle from previous studies (Anderson, 1982;Brown & Shankland, 1981;Katsura et al., 2010;Ono, 2008;Stacey & Davis, 2008;Stixrude & Lithgow-Bertelloni, 2011). All models assume an adiabatic temperature distribution through the mantle and some models include a TBL at the Earth's surface. A value for T CMB is obtained by adding a lower thermal boundary layer (TBL) to the adiabatic temperature in the lowermost mantle (TBL not shown). The study of Stixrude and Lithgow-Bertelloni (2011) gave an adiabatic temperature of 2692 K at a depth of 2,891 km. A very similar value was obtained by Ono (2008). We use the Stixrude and Lithgow-Bertelloni (2011) adiabat throughout the paper for calculating the TBL model at the core-mantle boundary.
was obtained from the model proposed by Anderson (1982) that included a temperature jump of 150 K at the base of the upper mantle.
In contrast to constraining anchor temperatures using the phase transitions of individual silicate minerals, Stixrude and Lithgow-Bertelloni (2011) focused on the phase equilibria of multiphase assemblages in the mantle. An algorithm for global minimization of the Gibbs free energy was used to self-consistently compute the isentropes of a six-component pyrolite model of the MORB source region (Workman & Hart, 2005), which can serve as a useful reference state from measuring radial and lateral temperature variations in the mantle. By comparing with experimental phase equilibria data, the self-consistently computed 1600 K isentrope was used to illustrate the mantle geotherm, yielding an adiabatic temperature 2692 K at the CMB depth of 2,891 km (Stixrude & Lithgow-Bertelloni, 2011). This is the mantle geotherm model that will be used here.

Core-Side Temperature Constraints
The temperature on the core side of the TBL can be estimated by anchoring the temperature to physical properties of the core itself. Commonly, the temperature at the CMB is estimated by determining the melting temperature of core material at the ICB and then integrating over the adiabatic profile through the outer core. This is reasonable given the well-stirred nature and very low viscosity of outer core material. The upper-bound of T CMB comes from assuming the core is made of pure iron, which has an estimated melting temperature at the ICB of 6230 ± 500 K (Anzellini et al., 2013). This results in a T CMB of 4700 ± 500 K, which is near the liquidus of peridotite and chondritic materials at CMB pressures Fiquet et al., 2010). Since large volumes of melt are not observed in the lower mantle, the CMB temperature cannot significantly exceed the solidus temperature of the mantle, and thus assuming the core to be pure iron leads to too high a CMB temperature.
C. Davies et al. (2015) developed three core compositional models consisting of Fe-O-Si alloys based on the density jump at the ICB. For simplicity, we chose to use the three representative ICB temperatures that result from these models. When extrapolated along an adiabat, the model inner core temperatures give T CMB of: 3900 K (∼17 wt% O and 2 wt% Si), 4100 K (∼13 wt % O and 8 wt%Si), and 4300 K (∼8 wt% O and 10 wt% Si). All of these values for T CMB fall acceptably close to the solidi of CI chondrite and peridotite-type material at CMB pressures (∼4180 K), but are several hundred degrees higher than the solidus of pyrolite-type material at the same condition (3570 ± 200 K; Nomura et al., 2014). Importantly, they are consistent with previous models of this type (e.g., Olson et al., 2015), which we use for comparison to our model in Section 3.

Section Summary
In this section, we summarized constraints and reasonable ranges on the temperatures at the top of the core and the bottom of the mantle from petrology and mineralogy. We considered several mantle adiabats to determine the . Mantle potential temperature (T P ) derived from the composition of primary magmas. 1,167 samples were compiled from the continental igneous rock record in the GEOROC database. Fifty-four samples (blue circles and triangles) have MgO >9 wt%, and SiO 2 >45 wt%. Thirty-three samples reported by Herzberg and Gazel (2009) are also included (yellow circles). Primary magma compositions and mantle T P values were modeled using the PRIMELT (Herzberg & Asimow, 2015). The red and orange curves represent the evolution of T P with present-day Urey ratio of 0.2 and 0.3, respectively, based on a secular thermal model back projected from a present-day mantle temperature of 1625 K and assuming whole mantle convection, fully described in Korenaga (2008) and updated by Korenaga (2017). The present T P is assumed at 1625 K (star symbol) with uncertainty of ±50 K Herzberg et al. (2010). G. F. Davies (2009) (black dashed curve) and Ganne and Feng (2017) (purple dashed curve) show significant differences in the overall cooling rate resulting from different parameterization (e.g., Urey ratios), which propagate to the present-day T P at 1575-1625 K. The present day is indicated as 0 Ga. temperatures profile throughout the mantle. Using these temperatures above and below the CMB and within the mantle, we can estimate the temperature change across the CMB and infer the associated Q CMB .

Models of Mantle Convection
Three-dimensional numerical simulations reproduce many hypothesized aspects of mantle convection (e.g., Heister et al., 2017a;Zhong et al., 2000) and can even replicate structures of similar size and shape to the seismically observed LLSVPs (McNamara & Zhong, 2005). Numerical models have included depth-and temperature-dependent material properties, phase transitions, radiogenic heating, compositional heterogeneity, and boundary conditions that reproduce surface plate motions (e.g., Nakagawa & Tackley, 2008;Olson, 2013), allowing us to explore how the CMB heat flow varies for a wide range of thermal conditions and physical properties. However, the computational demands of high-resolution, three-dimensional models would necessarily limit our search of parameter space. Instead, we use a simple theoretical boundary layer model to estimate the CMB heat flow for a range of thermal conditions. A single three-dimensional mantle convection model can be used to anchor or calibrate the theoretical model, allowing us to extrapolate to other conditions. This approach is only valid if we can demonstrate the validity of the boundary layer theory when realistic complications, such temperature-dependent viscosity or compositional heterogeneity are included. In this section, we show how the boundary layer model can be calibrated using a recent mantle convection model by Olson et al. (2015). In the next section, we demonstrate that the boundary-layer theory can be adapted to more realistic complexity.
The model of Olson et al. (2015) adopts a (superadiabatic) temperature change of 2500 K across the mantle. This value represents the sum of temperature changes across the boundary layers at the top and bottom of the mantle, and lies near lower limit of estimates given in the previous section. The goal is to extrapolate these results to other conditions. The average, present-day CMB heat flux from the simulation of Olson et al. (2015) is 84.8 mW m −2 , with regional and temporal variations predicted. The model simulates internal convection with subduction at the surface and dense piles sitting at the core-mantle boundary (formed from an initial dense layer placed in the model to emulate the LLSVPs). Higher-than-average CMB heat fluxes are often reported in regions below subduction and away from areas where the dense material is piled. Lower CMB heat fluxes are typically found below the dense, LLSVP-like material. Here, we develop a model to extrapolate CMB heat flux to a broader range of conditions.

Thermal Boundary Layer Model
We adopt the model of Olson et al. (2015) as a reference solution and use the simple TBL model to extrapolate the predictions. The TBL model is based on the assumption that convection in the vicinity of one boundary is independent of processes near the other boundary. This assumption leads to a simple relationship between the Nusselt number, Nu, and the Rayleigh number, Ra. Here, Nu is the ratio between the convective heat flux and the heat flux conducted in a static fluid. The Rayleigh number describes the strength of convection and is defined by where ρ is the fluid density, α is the coefficient of thermal expansion, κ = k/ρC P is the thermal diffusivity, k, is the thermal conductivity, η is the fluid viscosity (in Pa·s), and ΔT is the (superadiabatic) temperature change across the depth, d, of the fluid. On dimensional grounds, the relationship between Nu and Ra takes the form where q is the convective heat flux and c is a constant. Howard (1964) estimated Nu for the mantle, Nu m , by arguing that the TBL grows to the point of becoming locally unstable, as measured by the critical value for the Rayleigh number (denoted Ra c ). The result is = (2 ) −1∕3 , where Ra c ≈ 700 gives the constant c = 0.089 proposed by Kraichnan (1962). A slightly larger constant, c = 0.124, is recovered from experiments (Niemela et al., 2000).

Rearranging Equation 6 for the heat flux (heat flow per unit area) gives
where we introduce the kinematic viscosity, ν = η/ρ, and define the constant A as Recall that ΔT is the sum of the temperature change across the top and bottom boundary layers. It is more common to express the heat flux at each boundary in terms of the temperature drop across each of the two boundary layers. We derive the required expression by considering a much simpler configuration. Convection in a simple fluid layer has identical temperature changes across the top and bottom boundaries (i.e., ΔT TBL = ΔT/2). The resulting relationship between heat flow and ΔT TBL suggests that the heat flux across the CMB can be written as where ΔT CMB is the temperature drop across the TBL (ΔT TBL ) at the base of the mantle.
We compare the predictions of the TBL model with the results of the thermochemical model of Olson et al. (2015). An estimate of ΔT TBL at the top of the mantle is given by the mantle potential temperature of T P = 1600 K (Stixrude & Lithgow-Bertelloni, 2011) minus the surface temperature. Taking ΔT TBL = 1325 K implies ΔT CMB = 1175 K when ΔT = 2500 K. Using the experimental estimate for c, together with the material properties from Olson et al. (2015) (see their Table 1 or our Table 2), gives q CMB = 58.7 mW·m −2 , which is lower than the average value of 84.7 mW·m −2 from Olson et al. (2015). There are several reasons why the prediction of TBL model might depart from the numerical model. The influences of spherical geometry and temperature-dependent viscosity could have significant roles. Here we adopt a nominal lower-mantle viscosity of 10 22 Pa·s in the TBL model. A reduction in the average viscosity in the hot TBL (≈3 × 10 21 Pa·s) would be enough to reconcile the heat fluxes. Other possible factors include deviations in the average radial temperature in the interior of the mantle from an adiabat, particularly in the presence of compositional heterogeneity. A subadiabatic internal temperature (Bunge et al., 2001) would yield a larger value for ΔT CMB , and might also account for the higher value of q CMB in the numerical model.
Despite the simplicity of the TBL model, we find reasonably good agreement with the numerical model of Olson et al. (2015). We absorb the deficiencies of the TBL model into a change in the constant A and use this calibrated model to explore the dynamics of the mantle around the conditions adopted in Olson et al. (2015). Changing the value of A in Equation 8 by a factor of 2 brings the TBL model into agreement with the numerical model of Olson et al. (2015) when we adopt a constant viscosity. Separate assessments of the complications due to temperature-dependent viscosity and chemical heterogeneity are given in subsequent sections.
A prediction for the time-averaged temperature profile across the TBL is given by a solution of the thermal conduction problem where erf() is the error function, z is the height above the CMB, and l is the time-averaged thickness of the TBL. The average thickness is recovered from the CMB heat flux using q CMB = kΔT CMB /l and the time required for the TBL to grow by diffusion to the average thickness is τ = l 2 /πκ. Representative values are l = 55.5 km and τ = 43.3 Ma for the parameters adopted in Olson et al. (2015): ΔT CMB = 1175 K, k = 4 Wm −1 K −1 , and q CMB = 84.7 mW·m −2 . About 90% of the temperature drop across the TBL occurs over a distance of 1.3l from the CMB or about 70.6 km for our representative estimates.
We model the TBL using Equation 9. In Section 5, we predict the seismic anomaly that would results from of these temperature profiles and assess the observability of such a layer by various seismic methods. We then calculate a range of q CMB values, and compare these with the observed mantle cooling rate. We estimate q CMB k 4 W/mK ρ 5,570 kg/m 3 C P 1,000 J/kgK α 2 × 10 −5 K −1 η 1 × 10 22 Pa ⋅ s Table 2 Mantle Properties From Olson et al. (2015) using three representative values for ΔT CMB derived from the difference between the mantle adiabat of (Stixrude & Lithgow-Bertelloni, 2011) and low, moderate, and high estimates for the CMB core temperature based on different estimates of the melting curve of Fe-O-Si alloys: T CMB = 3900, 4100, and 4300 K (Section 2.2.2). We use reference material properties from Olson et al. (2015) (see Table 2). For the low, moderate, and high cases, respectively, we obtain q CMB = 89.4, 109.8, and 131.2 mW m −2 , which correspond to globally integrated heat flows, Q, of 13.6, 16.6, and 20.0 TW ( Figure 4a). These values differ from the present-day estimate of Olson et al. (2015) only because we have adopted different values for T CMB . An error function solution across the TBL implies TBL thicknesses of 70 ± 9, 66 ± 9, and 63 ± 9 km for the low, moderate, and high core temperature TBL models ( Figure 4b). Uncertainties in thickness are fully propagated from the range of likely values of k, C P , and η. We test k between 3.5 and 5.5 (X. Tang et al., 2014), C P between 1,000 J/(kg K) (Olson et al., 2015) and 1,300 J/(kg K) (Akaogi & Ito, 1993), and η between 1 × 10 22 Pa s (Olson et al., 2015) and 5 × 10 22 Pa s (Cizkova et al., 2012). These thicknesses are 1.3l which is, as discussed above, representative of the distance above the CMB at which the error function solution for temperature accounts for 90% of the temperature drop across the TBL (Turcotte & Schubert, 2002). The time to grow to the average thickness for each case is 40.3, 36.4, and 33.3 Ma, respectively.
The largest source of uncertainty in the simple TBL model is due to the change in temperature across the TBL, ΔT CMB . By using a range of likely core temperatures, as described above, we are able to account for possible variations. The uncertainty in the viscosity of the TBL is the next largest contributor to uncertainty in the thickness and stability of a TBL. Even though the viscosity appears as ν 1/3 in Equation 9 for the heat flux, the allowable range of viscosity for lower mantle materials is large. Changing the average viscosity of the TBL from 10 21 Pa⋅s to 10 24 Pa⋅s in our model yields heat fluxes between 238 mW·m −2 for a low viscosity and 23 mW·m −2 for a high viscosity for T CMB = 4100 K (implying ΔT CMB = 1408 K). Other examples can be computed from the ν 1/3 dependence. Comparable changes in the maximum TBL thickness and stability time scales are expected with the high and low viscosities. Knowledge of viscosity at the CMB is clearly important to estimate the CMB heat flow. The value of thermal conductivity used here is supported by experimental data (Geballe et al., 2020). Uncertainties in the thermal conductivity of the lower mantle are not expected to change our results significantly.

Comparison of Model Prediction With Observations
We use the observations of the Earth's secular cooling in Section 2.1 to evaluate our predictions of Q CMB from TBL theory. We specify the total radiogenic heat production (H tot ) by considering three representative values of Ur of 0.7, 0.3, and 0.08 based on geodynamical, geochemical, and cosmochemical arguments from Sramek et al. (2013), respectively. Each of our three estimates of Q CMB is based on a CMB temperature from the coreside constraints. We use Q CMB and the radiogenic heat production implied by the value of Ur to compute an instantaneous mass-averaged mantle secular cooling rate, which is further compared against the mass-averaged mantle secular cooling rate calculated based on the average mantle cooling rate in the last 2 billion years ( Figure 3). The instantaneous cooling rate of the mantle likely differs from its time average cooling rate in the last 2 billion years. However, we choose to ignore this difference for the sake of simplicity. The resulting mass-averaged mantle cooling rates are shown in Figure 5 as a function of the CMB temperatures. For example, a Urey ratio of 0.7 yields a temperature increase (i.e., Earth is heating over the geological time) for all CMB temperatures, although the uncertainties permit cooling for the lowest temperature of 3900 K. A Urey ratio of 0.08 is consistent with petrological constraints at the higher CMB temperatures considered here, whereas a intermediate Urey ratio of 0.3 favors relatively low CMB temperatures. Better constraints on T CMB would allow a better assessment of the Urey ratio at our predicted Q CMB .

Section Summary
In this section, we introduced a theoretical TBL model which, combined with our various estimates of the CMB temperature jump from Section 2 allows us to predict Q CMB and the thickness of the corresponding TBL. In subsequent sections, we test the validity and predictions of this model.

Testing the TBL Model With a Numerical Convection Model
Our simple theoretical TBL model was developed with assumptions that limit its applicability to the real Earth: for example, the assumption of uniform material properties. In the real mantle, both viscosity and thermal conductivity are expected to vary strongly with temperature. Second, there is an assumption of chemical homogeneity radially across the Earth. Radial heterogeneities in the Earth's composition will affect regional heat flow across the TBL (Olson et al., 2015), the expected regional temperature profiles, and the overall dynamics. For this reason, it is reasonable to question whether the standard TBL model (described in Section 3.1) is sufficient to predict changes in CMB heat flow when the adiabatic temperatures in the mantle and core are varied within allowable ranges. Laboratory experiments (Richter et al., 1983) suggest that the influence of temperature-dependent viscosity can be incorporated into the TBL model by setting the value of viscosity equal to the arithmetic average value at the top and bottom of the TBL. We test this approximation using a numerical convection model that allows large variations in viscosity over a small radial distance. We also use the numerical model to explore the influence of chemical heterogeneity at the base of the mantle. We demonstrate that the simple TBL model can be extended to include real-world circumstances, justifying our use of the model to extrapolate to other conditions. Our numerical model is based on the open source finite element code ASPECT (Heister et al., 2017a). We restrict the model to 2D, which permits large viscosity variations over a small radial distance and sufficiently high numerical resolution. One of the challenges in using strongly temperature-dependent viscosity is the development of a cold stagnant lid at the top of the convection model. In such a case, the style of convection is controlled by conduction through the stagnant lid (Moresi & Solomatov, 1995) and is unlike convection in Earth's mantle, where subduction of cold lithosphere cools the interior and maintains a larger temperature change across the lower boundary layer. For this reason, we focus on convection in the lower half of the mantle and represent cooling due to subduction with a uniform heat sink throughout the lower mantle. This setup avoids the cold stagnant lid in the calculation and allows us to investigate the development of thermal instabilities at the lower boundary in isolation from processes at the top boundary. In fact, the independence of the top and bottom boundary layers is fundamental to the TBL model. With this approach, we are able to test the influences of temperature-dependent viscosity in a setup that is conceptually consistent with the standard TBL model. Our numerical approach is unique relative to previous larger-scale models that incorporated a TBL at the CMB (Heister et al., 2017a;McNamara & Zhong, 2005;Olson et al., 2015) in that the much finer numerical resolution allows a more realistic representation of the 100-km thick TBL that was resolved over a much larger depth range in previous models. Consequently, our resulting vertical gradients of effective viscosity and temperature above the CMB are notably  Figure 3 constrain the secular cooling rate to a narrow range (red shaded bar). Note. We convert the time variation of the potential to that of the mass-averaged mantle temperature. Different Urey ratios (Ur) are based on geodynamical (red line), geochemical (blue line), and cosmochemical (green line) evidence taken from Sramek et al. (2013). CMB temperature is related to Q CMB through the thermal boundary layer model. Error bars of the cooling rate are calculated by propagating the uncertainties of Urey ratios (Sramek et al., 2013), the energy terms including the oceanic and continental heat fluxes (Jaupart et al., 2015), and the heat production of the continental crust (Rudnick & Gao, 2003). higher than these earlier studies, making the 2D models more appropriate for studying the dynamic behavior of the thin TBL. In addition, the TBL theory assesses the steady-state behaviors of convection, which is easier to achieve in 2D than 3D models.

Governing Equations for Mantle Convection
The mantle convection models solve the three conservation equations of mass, momentum, and energy, assuming the lower mantle to be an incompressible and Newtonian fluid: where u is velocity, P is dynamic pressure, η is dynamic viscosity, ρ m is density of the ambient mantle, α is the coefficient of thermal expansion, g is gravitational acceleration, T is temperature, κ is thermal diffusivity, ΔT = T − T 0 is the temperature anomaly relative to a constant reference temperature T 0 , and H is a heat sink.

Model Setup
We aim to resolve the nature of convection in the lower half of the mantle, thus all the simulations are performed in a 2D Cartesian box with an aspect ratio of 2:1, or 2,890 km × 1,445 km in the horizontal and vertical dimensions, respectively. Our models use an adaptively refined finite element mesh that is updated after every five time steps, and achieves a 1.4 km vertical resolution at the finest level of refinement. This provides adequate numerical resolution for properly capturing large dynamic gradients within fine-scale mantle features like the thin TBL (Moresi & Solomatov, 1995).
Boundary conditions on temperature include a constant CMB temperature on the lower boundary, and insulating conditions (i.e., zero heat flux) on the top boundary and side walls. A free-slip boundary condition is assumed at the top and bottom boundaries. The Boussinesq approximation means that the mantle adiabat collapses to a constant background temperature, such that T − T 0 represents the deviation from a mantle adiabat. Buoyancy forces in the calculation stem from (a) temperature perturbations relative to the constant reference value T 0 , and (b) chemical density anomalies associated with an initial compositional structure. The thermal component of buoyancy in our model mainly originates in the TBL.
The applied heat sink, H, continually cools the lower mantle and draws a heat flow across the CMB when this boundary is maintained at a fixed temperature. In reality, we do not expect subduction to behave as a uniform heat sink, although vigorous mixing should tend to homogenize the cooling. Our expectations are similar to the assumptions usually made to justify an adiabatic temperature profile in the convecting part of the mantle. The time-averaged heat flow across the lower boundary is set by the heat sink. In other words, the average CMB heat flow is an input parameter for the model. An important quantity recovered from the solution is the temperature drop needed to sustain this heat flow. We use the relationship between temperature drop and heat flow to assess the TBL model.

Models With Temperature-Dependent Viscosity
The influence of temperature-dependent viscosity is modeled by letting the local value of viscosity vary according to where η r (z) is the depth-dependent background viscosity, T and T 0 represent the local and reference temperature of the mantle, respectively, and β describes the magnitude of the temperature dependence.
Representative values for β are often inferred from a scaling based on the homologous temperature (e.g., Stacey & Davis, 2008). This relationship relates viscosity to the melting temperature by assuming the activation enthalpy of self-diffusivity, H diff , is only a function of melting temperature (Yamazaki & Karato, 2001). Despite the semi-empirical nature of the homologous temperature scaling, it has been shown to work very well for many mantle minerals . The exponent of temperature dependence (β) is related to H diff by β = H diff /RT, where T is the melting temperature at zero pressure and R is the ideal gas constant. For the case of diffusion creep in MgSiO 3 bridgmanite under dry conditions we expect β ≈ 20.6 (Ammann et al., 2010). The results for bridgmanite may be relevant for our simulations because the lower mantle is likely to be dry and dominated by diffusion creep (Bercovici & Karato, 2003;S. Karato, 1981).
Geodynamic models rarely achieve realistic values for the temperature dependence of viscosity, so we gradually increase β from 0 to 20 to simulate progressively stronger temperature dependence. Figure 6a shows a representative numerical calculation for the case of a strongly temperature-dependent viscosity (β = 20). Warm material rises from the TBL to form plumes, whereas cold material descends toward the TBL. The coldest temperatures are usually found in isolated regions immediately above the TBL.
A profile of the horizontally-and time-averaged temperature and viscosity is shown in Figures 6b and 6c as a function of distance outward from the CMB for the case β = 20. The thickness of the TBL is determined by the distance where large fluctuations are detected in the temperature profile (blue star). These fluctuations represent the influence of thermal anomalies due to warms plumes and cold downwellings, which are treated as part of the convective interior.
The representative simulation in Figure 6 has a horizontally averaged viscosity that varies by a factor of 17.5 across the TBL. Lateral variations in viscosity are also present in the calculation, particularly in the region above the TBL. Temperature heterogeneity above the TBL causes the ratio of highest to lowest viscosity to exceed a factor of 32. Although the value of the viscosity variation is smaller than previously reported (Heister et al., 2017a), our models actually resolve larger vertical gradients of viscosity and temperature, as is what the TBL theory evaluates in effect. Despite these large viscosity variations, we do not observe any evidence for a separate cell of convection within the TBL (e.g., Ke & Solomatov, 2004). It appears that a much larger viscosity variation is needed to promote an isolated layer convection within the TBL.

Models With Chemical Heterogeneity
Here we investigate the influence of potential chemical heterogeneities associated with the LLSVPs on the CMB heat flow. In the models, we introduce an initial 40-km thick chemically distinct layer at the base of the mantle with an excess compositional density. We compare models with excess densities of 0%, 5%, and 10% relative to the ambient mantle. All these models have moderate temperature dependent viscosity (β = 3). Introducing chemical heterogeneity at the base of the mantle alters the spatial distribution of heat flow. Figure 7 illustrates how the dense material is organized into piles by convection. Thicker layers of dense material accumulate below upwellings, which draw material from the surrounding boundary region. Similar structures are observed in previous calculations (e.g., McNamara & Zhong, 2005). Distinct vertical profiles of temperature develop in regions with thin and thick layers of dense material (see Figure 7b). Regions with thick TBL coincide with the location of dense piles, whereas thin TBL occur where dense material is replaced by ambient mantle. Even though the pile morphology is transient due to ongoing convection, the temperature above these thermochemical boundary layers is relatively uniform due to rapid convective mixing. We define a temperature drop across the TBL using the horizontally averaged temperature profile. The corresponding heat flow is affected by both the thickness and horizontal distribution of the dense material. We seek to understand how these influences can be represented in a TBL model.

Testing the TBL Theory
The TBL theory predicts that the heat flux, q CMB , depends on Δ 4∕3 . We test this dependence by plotting q CMB versus ΔT CMB calculated from our models to estimate the best-fitting exponent. A common extension of the TBL theory to account for temperature-dependent viscosity assumes that the viscosity in Equation 9 can be replaced with , the arithmetic average of the viscosity at the top and bottom of the TBL (Richter et al., 1983). Numerical solutions with both constant and temperature-dependent viscosity can be compared if we plot ̄1 ∕3 versus ΔT CMB . In both of these calculations we expect a slope of 4/3. A minor modification of the fit to include the influence of κ, as defined in Equation 7, allows us to connect the y-intercept of the log-log fit with the constant A in Equation 8. Figure 8 shows a log-log fit of ̄1 ∕3 −2∕3 to ΔT CMB . The slopes are obtained by linear least squares regression. Our numerical models reveal fits with slopes close to the predicted value of 4/3. The fit for the case of constant viscosity is consistent with theory at the 2σ level. A slightly steeper slope is obtained for the case of strongly temperature-dependent viscosity, although we have overlapping uncertainties in the slopes for β = 0 and 20. A narrower uncertainty for the strongly temperature-dependent case means that the theoretical prediction lies  slightly outside the 2σ uncertainty. One possibility is that the uncertainties at β = 20 are underestimated due to fortuitously small misfits to the linear trend. The reported uncertainty is entirely based on deviations from the linear trend. Accounting for uncertainty in identifying the top of the boundary layer and recovering ΔT from the numerical models would likely increase the error bars on the slope.
We also observe a small offset in the computed heat flow for the β = 0 and β = 20 cases. We cannot alter the definition of the average effective viscosity without altering the slope of the trend. This implies that the discrepancy is likely attributed to a small difference in the numerical constant A (see Equation 8). For example, a small change in the assumed value for the critical Rayleigh number Ra c could account for the offset. However, the consistency in the recovered slopes for the cases of constant and strongly temperature-dependent viscosity suggests that the use of a depth-averaged viscosity across the TBL is an effective strategy for predicting the influence of strongly temperature-dependent viscosity in the TBL model. Even though the numerical model does not reach the values of ΔT CMB expected in the Earth, the scaling in Figure 8a supports our use of the TBL model.
Compositional heterogeneity at the base of the mantle is another complexity in a realistic TBL model. Our numerical solutions with a dense chemical layer predict lower heat flow than the standard model at the same ΔT CMB . Figure 8b reveals an offset in q CMB when layers with 5% and 10% excess density are present at the base of the mantle. The slope of the log-log relationship between q CMB and ΔT CMB at 10% excess density is remarkably similar to that obtained from calculations without a dense layer. This suggests that the TBL model can reproduce variations in q CMB with ΔT CMB . The offset in the predicted heat flow suggests a connection between layer properties and the model constant A. Quantifying this relationship would allow us to explore the influence of different layer thicknesses on the CMB heat flow.
The chemical layer with an intermediate excess density (5%) exhibits a steeper slope of q CMB to ΔT CMB (1.65 ± 0.05) than either of the cases with 0% and 10% excess density. This distinct behavior may be due to the influence of entrainment by the overlying convection. The layer with intermediate excess density experiences more entrainment than the layer with 10% excess density. In addition the rate of entrainment increases at larger ΔT because the vigor of convection is increased. When the duration of each simulation (in model time) is comparable we can expect more total entrainment when ΔT is large. Thus, part of the change in q CMB at 5% excess density may be due to a gradual reduction in the residual thickness of the dense layer. Full entrainment of the dense layer at large ΔT (or at large simulation time) should eventually approximate the conditions with no dense layer.

Section Summary
In this section we numerically simulate a TBL using constraints derived from petrology and mineralogy (Section 2) and TBL theory (Section 3). We test extensions of the simple theoretical TBL model to allow for temperature dependent viscosity and compositional variations. We show that TBL theory can be extended to these two specific complications, although other factors, such as phase transitions and compressibility are not taken into account.

Seismological Manifestations of a TBL
Seismic observations of the lower mantle serve as a test of the TBL theory. We take the TBL model derived above, predict the impact on several seismic observables, and discuss the likelihood of such a layer being detectable by various seismic methods.
To determine the impact of a TBL on the elastic properties of the lowermost mantle, we calculate the mantle temperature profile and use BurnMan (Cottaar et al., 2014) to calculate the resultant seismic velocity anomaly. We solve the equation of state using the third-order method of (Stixrude & Lithgow-Bertelloni, 2005), average the effects of multiphase compositions, and model the likely reductions in seismic velocity and density (Figure 9). We construct temperature profiles for the low, moderate, and high values of ΔT CMB (Section 3.1). For these three ΔT CMB values, the TBL model would result in velocity reductions between 3%-4% and 1.4%-1.7% for dVs and dVp, respectively, and a density reduction of 1.4%-1.7%, with the strongest effect for the highest ΔT CMB . These values are calculated for a pyrolitic mantle, where the composition is fit to PREM in Cottaar et al. (2014) (molar fractions: 0.55 Mg bridgmanite, 0.04 Fe bridgmanite, and 0.41 Mg-Fe periclase). Since the exact composition of the lower mantle is not known, we also test the velocity reduction of our TBL model on other mantle These models result in at most a 0.3% difference in dVs, and much smaller differences for dVp and dρ. While compositional heterogeneity would affect the absolute seismic velocity, its effect on the velocity perturbation of the TBL is likely too small to be detectable in seismic data, thus compositional variation is neglected in our seismic models.
The TBL is likely to result in a very thin layer of reduced velocity with the largest dV at the CMB. Therefore, waves diffracted along the core-mantle boundary will be most sensitive to velocity anomalies in the TBL. The vertical sensitivity of diffracted waves depends on the seismic frequency, with higher frequency waves having sensitivity concentrated closer to the CMB. Diffracted waves have been used extensively to study the seismic structure at the CMB (e.g., Alexander & Phinney, 1966;Wysession & Okal, 1989, 2013Wysession et al., 1992). We model the influence of a TBL on diffracted waves using axiSEM (Nissen-Meyer et al., 2014), a rotationally symmetric waveform simulator based on the Spectral Element Method. We impose the velocity and density anomalies predicted by BurnMan on top of the 1D velocity and density structure of ak135 (Kennett et al., 1995) and simulate the seismic wavefield at a maximum frequency of 0.5 Hz. Several differences are visible for the TBL case compared with the 1D case (Figures 9d-9f): distorted waveforms, delayed arrivals, and visibility at greater source receiver distances. Travel time delays would be up to 1.5 s for the moderate ΔT CMB model. Considering the precision of travel time measurements, it is possible that with careful analysis the impact of the TBL would be observable in high frequency Pdiff and Sdiff data.
A velocity reduction from a TBL would also influence the angle of take-off of diffracted waves from the CMB. This would manifest at the surface as a change in slowness of the wave relative to the reference case, where slowness = sin 0 , i is the incidence angle and v 0 is the local seismic velocity. We determine the slowness anomaly that would be caused by a TBL by computing synthetic seismograms at closely spaced stations, akin to a seismic array, and use beam forming (D. Davies et al., 1971) to measure the slowness of Pdiff and Sdiff waves observed at synthetic seismic arrays. We determine that our standard TBL model would result in slowness increases of 0.05 and 0.1 s/° for P and S-waves, respectively. We confirm the influence of a TBL using a 1D infinite frequency ray-tracing approach. We again find that our standard TBL model would result in slowness increases of 0.05 and 0.15 s/° for P and S-waves, respectively. These slowness anomalies are very close to the detection threshold for real data, thus would likely be difficult to observe. Moreover, a lower ΔT CMB would cause a weaker velocity anomaly in the TBL, which would be even less likely to be observed as slowness anomalies in real data. However, while both slowness and travel time anomalies could also be produced by a broader velocity reduction in the lower mantle such as the LLSVPs, the sharpening of waveforms and observation at greater distances are both more diagnostic features of a TBL/thin low velocity layer.
On an even finer spatial scale, the amplitudes of body wave phases reflecting off the CMB are influenced by the reflection and transmission coefficients, thus are sensitive to the velocity and density contrast between the mantle and core (Persh & Vidale, 2004;Rost & Revenaugh, 2004). We compute the overall transmission and reflection coefficients and calculate the amplitude ratio and travel time of core-reflected phases (PcP, ScS, and PcS) on interaction with the CMB with and without a TBL, for a range of incidence angles (Figure 10). Amplitudes are very similar for models with and without the TBL, thus impedance contrast is barely sensitive to the TBL. However, travel time anomalies approach 2 s for a ScP waves affected by a TBL layer, which may be observable in real data. Therefore, while it would be possible to distinguish a TBL from an adiabatic model with reflected phases like ScP and PcS, it will depend strongly on data quality and observation accuracy.
Free oscillations of the whole planet, known as normal modes, have been used to investigate the broad-scale structure of the mantle (e.g., PREM, Dziewonski & Anderson, 1981). However, their limited vertical sensitivity means that they are likely not capable of resolving a TBL (P. J. Koelemeijer et al., 2012). In contrast, the splitting of normal modes trapped at the solid-liquid interface of the CMB, that is, Stoneley modes (Stoneley, 1924), have been used previously to sample the CMB and infer the velocity and density structure (P. Koelemeijer et al., 2013Koelemeijer et al., , 2017. Stoneley modes may be sensitive to a TBL. We predict the effect of our TBL model on the center frequencies of several normal modes with strong sensitivity to the CMB, and also Stoneley modes, whose sensitivity is only at the CMB. We use the normal mode summation code Mineos (Masters et al., 2011) to calculate the perturbation to normal modes caused by PREM alone and PREM with our moderate TBL case, which we compare with observed mode frequencies (P. Koelemeijer et al., 2013). The observations are not well matched by PREM, and are better matched by our TBL model (Figure 10e), with the exception of several of the higher frequency modes (2S16, 2S25, and 3S26). The differences between the reference case and our TBL model are greater than the measurement error. While not conclusive on Earth structure, this demonstrates that a TBL would perturb modes enough to be detectable, and suggests that both Stoneley modes and other CMB-sensitive modes require a low velocity layer close to the CMB.
We compare predictions for 1D reference models with 1D reference models plus a TBL. Since reference models are constructed from data, they may include the influence of a TBL and thus this comparison may seem redundant. Figure 10. Variation of (a and b) amplitudes and (c and d) travel times of core reflected phases (PcP and ScP) as a function of incidence angle for the TBL model (red) and the reference model (black). (e) Center frequency anomaly observations (observation-prediction) for a selection of normal modes most sensitive to the core-mantle boundary, with predictions from a preliminary reference Earth model (PREM; black) and a PREM Earth with a TBL with 1.3l = 63 km (red). Error bars about zero show the measurement error (P. Koelemeijer et al., 2013). However, we use ak135 and PREM (Dziewonski & Anderson, 1981;Kennett et al., 1995), neither of which are constructed using diffracted waves, and PREM uses only a limited catalog of normal modes. Given what our simulations above show above, the reference models are not likely to be very sensitive to the TBL and so our method is appropriate.

Section Summary
We use our model of the temperature across the TBL (Section 3) to predict the likely impact on seismic velocities. We then test the sensitivity of various seismic observational tools to the TBL. Our analysis suggests that the perturbations to the elastic parameters of the lowermost mantle caused by realistic estimates of a TBL would be on the edge of our detection thresholds for reflected and diffracted waves, but may be clearer with certain normal modes. With present data, seismology only provides a loose upper bound on Q CMB .

Constraints From Power Requirements for the Dynamo
The structure and stability of the TBL in the mantle at the CMB is intrinsically linked to the energy budget of the core. Heat flow across the CMB drives convective motion of the outer core's liquid iron alloy, which is responsible for the long-term existence and observed variations of the geomagnetic field. Higher heat flow (i.e., high Q CMB ) results in more energetic outer core convection, providing more power for the generation of the geodynamo, but also implies rapid cooling rates and very hot temperatures for the planet going back in geological time. Any model of the TBL at the CMB must be consistent with the energy budget of the core and the observation of the geomagnetic field, but also lead to inferred past temperatures for the mantle that are reasonable given the properties of silicate materials. Geomagnetic evidence suggests that the temperature gradient across the outer core is close to adiabatic (i.e., Q CMB ≃ Q ad ); where the temperature gradient exceeds the adiabat there will be convection, whereas there may be a thermally stratified layer where the gradient is shallower than the adiabat (e.g., Braginsky, 1993;Lister & Buffett, 1998). In this section we assess the constraints on CMB heat flow from the magnetic field.
Presently, the geodynamo is driven by a combination of thermal and compositional convection induced by secular cooling and the solidification of the inner core. The energy balance of the core can be described by where Q S denotes secular cooling of original accretional heat, Q L is latent heat resulting from the solidification of core material at the ICB, Q G is the release of gravitational potential energy as the inner core expels light elements and the light material buoyantly rises, and Q R is radiogenic heat production.
The adiabatic heat flow across the CMB, Q ad , can be determined from the adiabatic temperature gradient at the top of the core ( ) : Here, r is the radius of the core and k is the thermal conductivity of outer core material at the CMB. In the case of the core, we define the thermal gradient as the derivative of the adiabatic temperature profile: where γ is the Grüneisen parameter and ρ is the density. Taking ρ CMB = 9,906 kg/m 3 from PREM (Dziewonski & Anderson, 1981;Labrosse, 2015), γ = 1.5 (Vočadlo et al., 2003), T CMB = 4100 K, and k = 100 W/mK, we find Q ad to be 15 TW.
Notably, there is no term in the energy budget of the core (Equation 15) for the generation of the magnetic field. Ohmic dissipation, Φ, is the conversion of electrical to magnetic energy that occurs entirely within the core, so the global energy budget is not affected. However, dissipation is a non-reversible process, so it contributes to the entropy budget of the core. The entropy budget is defined as where E k is the entropy associated with the adiabatic conduction of heat out of the core, E Φ is the entropy of ohmic dissipation, and the remaining subscripts correspond to those in Equation 15. The energy and entropy budgets of the core are related by the core cooling rate, , to yield the following relationship: We choose to neglect the radiogenic heating terms in Equation 19. This is justified by recent metal-silicate partitioning experiments showing that the radiogenic heat production within Earth's core is likely negligible, particularly at present day because much of the radioactive material has decayed away (Blanchard et al., 2017;Bouhifd et al., 2013;Chidester et al., 2017;Faure et al., 2020;Malavergne et al., 2007). Importantly, Equation 19 is valid while there is an inner core. Prior to inner core solidification, the terms related to the growth of the inner core (Q L , E L and Q G , E G ) disappear, and ohmic dissipation is maintained entirely by secular cooling and opposed by conduction along the adiabat.
To determine whether our mantle TBL model is consistent with the generation of the magnetic field at present day, we balance the energy budget using Equation 19 and our estimates of T CMB from Section 2 and the corresponding values for Q CMB and Q ad from Section 3.1. We calculate the entropy and energy production within the core using the fourth-order parameterization in density following Labrosse (2015), with the assumption that the outer core is well-stirred and that the radial heat flux into the mantle is laterally homogeneous. We use the same material properties for the outer core material as Labrosse (2015), except we take the core thermal conductivity to be independent of depth within the core and we test a range of thermal conductivity values at the CMB.
Experiments to measure the thermal conductivity of core materials at high pressures and temperatures are incredibly difficult, primarily due to the small sample size, low melting temperature of available experimental tools (e.g., conductive leads), and uncertainties in converting the measured electrical conductivity to thermal conductivity via the Wiedemann-Franz law. The range in accepted experimental values of k of pure, solid Fe are 33 W/ mK (Konôpková et al., 2016) and 226 + 71/−31 W/mK (Ohta et al., 2016). The effect of light elements is yet uncertain, but one study has shown that the addition of Si could decrease the thermal conductivity of Fe alloys significantly (Hsieh et al., 2020). On the other hand, ab initio calculations on the transport properties of Fe alloys suggest that the Wiedemann-Franz law does not hold under the pressures and temperatures of core conditions and k = 80-200 W/mK (de Koker et al., 2012;Gomi & Hirose, 2014;Gomi et al., 2013Gomi et al., , 2016Pozzo et al., 2012Pozzo et al., , 2013Xu et al., 2018). Given the large uncertainty in this value, and particularly with the difficulty in interpreting experimental observations, we chose to use the lower range of ab initio values 77-130 W/mK.
The results of our energy budget calculation are presented by plotting the entropy of ohmic dissipation as a function of Q CMB (Figure 11). The bold lines in Figure 11 were calculated at our test CMB temperatures with a thermal conductivity value of k = 100 W/mK, similar to that of Pozzo et al. (2012). We also give an upper bound on E Φ using k = 77 W/mK (Xu et al., 2018) and a lower bound using k = 130 W/mK (Nimmo, 2015b). A dynamo in the outer core requires that E Φ is at least greater than zero (above the gray shaded region), but most dynamo models suggest that the strength of the current magnetic field requires E Φ of ∼100 MW/K at minimum (B. A. Buffett, 2002;B. A. Buffett & Christensen, 2007;Christensen, 2010;Christensen & Tilgner, 2004;Gubbins et al., 2003;Labrosse, 2003;Nimmo, 2015a;Roberts et al., 2003;Stelzer & Jackson, 2013). Our TBL model (gold box in Figure 11) is consistent with the ohmic dissipation required to explain the presence of the magnetic Figure 11. Core entropy production as a function of core-mantle boundary (CMB) heat flow after Labrosse (2015). Solid bold curves are calculated at each of our estimated CMB temperatures with k = 100 W/mK, while lower and upper bounds (dashed curves) are given for k = 130 and 77 W/mK, respectively. Q ad (gray box) is the adiabatic heat flow that covers the range of feasible CMB temperatures (3900-4300 K) and thermal conductivities (77-130 W/mK). The range of CMB heat flows of our TBL model is shown by the gold region (with specific values of 3900, 4100, and 4300 K marked as dotted vertical lines). field for any of the T CMB values. Indeed, 100 MW/K of ohmic dissipation can be supported at the present day for any heat flow above ∼7.5 TW. This is largely due to the energy sources derived from the growth of the inner core.
The Q CMB values estimated from our TBL model are high, suggesting that the core is superadiabatic for most conditions explored (i.e., Q CMB ≥ Q ad ). Such high values are sufficient to keep the core well mixed. However, it is also possible for the core to be slightly subadiabatic, which permits a stably stratified layer to develop on the core side of the CMB (Gubbins et al., 1982, for example).

Constraints on Core Thermal State From Paleomagnetic Observations and Geodynamo Simulations
The thermal state of Earth's core is controlled by the ratio of the heat flux across the CMB and the heat conducted along the core adiabat, which as introduced in Section 3.1, is called the Nusselt number. For the outer core Nu o = Q CMB /Q ad . If Q CMB > Q ad , the outer core will convect throughout and its fluid will be well-mixed, but if Q CMB < Q ad then a portion of the outer core will be stably stratified. Convection will occur where the actual thermal gradient is steeper than the core adiabat.
A non-convecting layer in the core will have implications for the strength and structure of the magnetic field that can be generated by the convecting portion of the outer core below. The paleomagnetic record places constraints on Earth's magnetic field strength and structure through geologic time. The morphology and variations of the geomagnetic field help estimate the thickness of a stratified layer, and therefore the Q CMB for a current adiabatic heat flux Q ad (15 TW, calculated in Section 6.1 with k = 100 W/mK). The fraction of the core radius which is thermally stratified, d, is a function of the Nusselt number d ≈ 0.3(1 − Nu o ; Lister & Buffett, 1998). A discussion of further geomagnetic constraints on core stratification is found in Section 7.3.
Numerical dynamo simulations with stratification below the CMB show that stratified layers thicker than ∼400 km result in magnetic field structure dissimilar to Earth's, putting a lower limit on Nu o of 0.7 and a lower bound on Q CMB of 10.5 TW, assuming Q ad = 15, indicated by the red line and square in Figure 12 (Christensen, 2018;Olson et al., 2017). Simulations with stratified layers thicker than ∼1,000 km do not maintain a dynamo indicated by the black line in Figure 12 (Christensen, 2018;Nakagawa, 2015;Olson et al., 2017).
Observations of the paleomagnetic field morphology and reversal rate also help constrain how superadiabatic we believe the core to be. Numerical dynamo simulations are controlled by a set of non-dimensional parameters defined by the ratios between various forces which set the balances between terms in the magnetohydrodynamic equations the simulations solve (e.g., Rayleigh number: = buoyancy dif f usion , Ekman number: = viscous dif fusion rotation , Prandtl number: = viscous dif fusion thermal dif fusion , and Magnetic Prandtl number: = viscous dif f usion magnetic dif f usion ). In rapidly rotating dynamo simulations the convective heat transfer efficiency (Nusselt number, Nu, an output) scales with the strength of convection (Ra, an input) as ≈ ( ∕ ) 6∕5 (Christensen, 2010;King et al., 2010). Dynamo simulations with a range of Ra values display a transition in structure and behavior from dipolar-nonreversing to multipolar-reversing fields. Vigorous convection that is not dominated by rapid rotation (i.e., Ra and E are both too high) produces magnetic fields dissimilar to Earth's. Dynamo solutions are attainable with highly superadiabatic conditions (Nu ≫ 1), but their fields do not look like Earth's. It is difficult to place a precise upper limit on how superadiabatic Earth is by scaling from geodynamo simulations that use smaller than Earth-like Coriolis forces (i.e., too large Ekman numbers). Earth likely resides near the transition as it is dipolar in structure and reversing (Christensen, 2010;King et al., 2010;Olson & Christensen, 2006). Where this transition occurs in Nusselt and Rayleigh number space suggests Nu > 1 and Q CMB > Q ad for Earth (Kutzner & Christensen, 2002;Olson & Christensen, 2006). However, the sharpness of the transition to multipolar magnetic fields implies the thermal profile of Earth's core is not highly superadiabatic, that is, Q CMB ≈ Q ad . Figure 12. If the outer core is subadiabatic (Q CMB < Q ad , i.e., everything below the blue line), a stable thermally stratified layer will form at the coremantle boundary. The thickness of this layer is a function of the Nusselt number (Lister & Buffett, 1998). Constraints on the thickness of the stable layer from geomagnetic observations and geodynamo simulations may then be used to estimate Q CMB for an assumed Q ad . Estimates of Q CMB from the mantle models and of Q ad from conductivity models and measurements likely place Earth's core within the dashed rectangle.

Section Summary
In this section, we use our theoretical TBL model (Section 3), paleomagnetic observations of the strength and behavior of Earth's magnetic field, mineralogically and experimentally derived values of core thermal conductivity, and predictions from geodynamo simulations to constrain Q CMB .

Discussion
Here we investigate the connections between the different disciplines introduced above and the implications of constraints on the present day Q CMB .

Long-Term Thermal Evolution of Earth's Core
As discussed in Section 6.1, the CMB heat flows predicted by our TBL model in Section 3.1 (13.6-20.0 TW) provide sufficient energy to drive Earth's geodynamo at present regardless of the thermal conductivity of Earth's core material (Figure 11). Most of the energy presently available comes from the solidification and growth of the inner core in the form of latent heat (Q L ) and gravitational potential energy (Q G ; Nimmo, 2015a). This mechano-chemical source of energy allows the core heat flow to be subadiabatic and still provide sufficient energy for the geodynamo, because it is not relying on thermal convection alone. However, recent thermal evolution models (e.g., Labrosse, 2015;Nimmo, 2015a) suggest that the inner core began solidifying as recently as 500 Ma, while the recorded paleomagnetic field is much older, raising questions about the magnitude of heat fluxed across the CMB, and thus the viability of our TBL model, in the past. We assess paleomagnetic evidence for the longevity and intensity of the magnetic field to understand the core's ongoing contribution to the Earth's heat budget.
Geodynamo studies have been conducted to predict the geomagnetic signal of inner core nucleation and thermal evolution. It is an important target because independent confirmation of inner core nucleation would help constrain the history of Earth's dynamo and thermal evolution, and more fundamentally our knowledge of core conductivity and thus the heat flux out of the core. P. E. Driscoll (2016) found a transition from a multipolar to dipolar dynamo regime followed by another transition to a weak-field dynamo as Earth cools before inner core nucleation, and a strong-field dipole following inner core nucleation. Landeau et al. (2016) predict little signal of nucleation manifested in surface field intensity with a dipole-dominant strong-field dynamo likely both before and after inner core nucleation, with a slightly larger octopolar field contribution following inner core nucleation. The present day inner-core:outer-core aspect ratio is 0.35. Lhuillier et al. (2019) ran a suite of chemically driven dynamo simulations with a range of inner core sizes and found a transition from a "small inner-core" regime to a "large inner-core" regime at an aspect ratio of 0.2 where the field was weaker and polarity reversals more frequent.
Paleomagnetic studies have probed the rock record for these signals. Bono et al. (2019) infer an inner-core age of 565 Ma from low paleointensity values recorded by magnetic inclusions within single silicate crystals. In a review of all existing paleointensity data, Biggin et al. (2015) propose inner core nucleation occurred at ∼1.2 Ga based on an observed increase in paleointensity amplitude and scatter at 1.3 Ga due largely to measurements from Thomas (1993). Kodama et al. (2019) resampled the study area that contributed the high intensity values at 1.3 Ga, and did not reproduce the anomalously high paleointensity of Thomas (1993). There is not currently evidence of a clear, smoking-gun paleomagnetic signal of inner-core nucleation, and there are still many large time gaps in the database (Smirnov, 2017).
The dynamo has been active for most of Earth's history regardless of the driving mechanism. Prior to inner core solidification, the heat flux out of the core is described by the equation So, prior to inner core nucleation, the heat flow out of the core and into the mantle must have been superadiabatic to support a thermally driven dynamo. Here, we use our core energy budget from Section 6.1 to infer the changes in Q CMB and T CMB through geologic time.
Geodynamic calculations that simulate tectonic plate motions suggest that Q CMB has not varied significantly over the recent geologic past (Olson et al., 2015). Holding Q CMB fixed at our intermediate value (16.6 TW) and k at a moderate 100 W/mK and integrating the energy sources back through Earth's history leads to an inner core age of ∼570 Ma (solid black curve, Figure 13), consistent with previous models and paleomagnetic observations (Bono et al., 2019;Labrosse, 2015). In this case, the adiabatic heat flow, Q ad , is 15 TW at present and surpasses Q CMB at 1.36 Ga. Without the energy sources available from inner core growth, the entropy associated with ohmic dissipation, E Φ , and thus the magnetic field, drops below zero as the heat flow becomes subadiabatic. The age of the inner core is controlled by Q CMB , where higher heat flow results in a higher rate of cooling for the planet and a younger inner core. If instead we use our upper-bound CMB temperature (4300 K) and associated heat flow (20.0 TW), the inner core is slightly younger (473 Ma), but we can sustain a dynamo as far back as ∼2.5 Ga (gray curve, Figure 13). An inner core older than 1 Ga would require CMB temperatures and heat flows well below those estimated here.
Alternatively, if we force a marginal dynamo (E Φ = 0 MW/K) in our intermediate case and allow the heat flow to vary prior to inner core solidification (orange curve in Figure 13), we find that Q CMB first drops to match the adiabatic heat flow defined by Equation 16, then steadily rises as the temperature rises going back in time. In this scenario, the inferred initial CMB temperature 4.5 Ga is slightly higher than the steady heat flow case, and Q CMB exceeds ∼22 TW. On the other hand, dynamo simulations suggest that the entropy required for a magnetic field of Figure 13. The thermal evolution of the core based on heat flows from our thermal boundary layer model after Labrosse (2015) showing (a) Q CMB , (b) E Φ , (c) inner core radius, and (d) T CMB . Q CMB is considered to be constant while the inner core exists. Either Q CMB is held constant for the entire model (gray and black solid and dashed curves) or E Φ is held constant prior to inner core solidification and Q CMB is allowed to vary naturally (orange and purple curves). There is no magnetic field when E Φ drops below zero.
present-day intensity is at 100 MW/K at a minimum (B. A. Buffett, 2003;Roberts et al., 2003). If we hold E Φ at 100 MW/K prior to the solidification of the inner core (purple curve, Figure 13), Q CMB increases at the transition to inner core solidification by several terrawatts to match Q ad . The maximum heat flow in this scenario exceeds 29 TW at 4.5 Ga. Note that these sharp changes in Q CMB at the onset of inner core growth are not realistic, since the heat flow out of the core is dependent on the mantle receiving the heat, not the presence of the inner core. This structure in the models simply highlights the need for Q CMB to be at least equal to Q ad for the geodynamo.
The thermal evolution of the core, and thus the CMB, is strongly dependent on the thermal conductivity of outer core material, k, because the thermal conductivity is directly related to adiabatic cooling (see e.g., Equation 16). As discussed in Section 6.1, this value is not well-constrained. All of the thermal evolution models presented above were run holding k constant at 100 W/mK throughout the depth of the outer core. If instead we lower the thermal conductivity in our moderate test condition to our lower-bound estimate of k = 77 W/mK (Xu et al., 2018), we observe E Φ > 0 through all of Earth's history (black dashed curve, Figure 13b), while the heat flow, age of the inner core, and T CMB remain unchanged from the higher thermal conductivity model (solid black curve, Figure 13). This highlights the importance of narrowing the range of possible thermal conductivities and melting temperatures of core material in future studies. Constraining these material properties is crucial to developing a firm understanding of the Earth's thermal evolution.
The major limitation of all of our thermal evolution models, and many recent models in the literature that could explain the long-term behavior of the magnetic field (e.g., Labrosse, 2015), is that the required CMB temperatures early in Earth's history are much higher than the solidus of peridotite material at lower mantle conditions (Fiquet et al., 2010). This would result in a very thick and long-lived molten layer within Earth's lowermost mantle, which is inconsistent with the petrology and geodynamics estimates (see Section 3 for details). Here, we discuss several alternative scenarios that could explain the magnetic field without requiring unacceptably high early temperatures. First, it is possible that there are additional sources of energy within the core that would balance the energy budget while providing a geodynamo, allowing the core to cool much more slowly such that the inferred initial temperatures were much more reasonable. Work to identify possible energy sources within the core that may alleviate this problem is ongoing. As mentioned in Section 6.1, the possibility that radioactive decay was an important energy source in the core has been refuted (Blanchard et al., 2017;Bouhifd et al., 2013;Chidester et al., 2017;Faure et al., 2020;Malavergne et al., 2007). More recently, it has been proposed that strongly lithophilic elements, such as Mg and Si, dissolved into the core at the very high temperatures of core formation and subsequently precipitated out as oxides or silicates as the core cooled (Badro et al., 2016;Hirose et al., 2017;Mittal et al., 2020;O'Rourke & Stevenson, 2016;O'Rourke et al., 2017). Since precipitation is mostly likely to occur at the top of the outer core, this would provide a strong source of gravitational potential energy as the much denser supernatant metallic liquid will gravitationally sink. Whether these lithophile elements will partition into the core in high enough concentrations to saturate and precipitate out is an ongoing area of research (Badro et al., 2018;Du et al., 2017). It has also been proposed that mechanical stirring arising from tidal coupling in the Earth-Moon-Sun system may drive a dynamo, subverting the need for a thermally driven dynamo prior to inner core nucleation (Andrault et al., 2016;Reddy et al., 2018). An extreme alternative would bypass the geodynamic constraints on lower melting entirely: if the lower mantle was well above the liquidus temperature, it is plausible that a magnetic field was generated by conductive silicate liquids (Blanc et al., 2020;Stixrude et al., 2020;Ziegler & Stegman, 2013), which would support the very hot CMB early in Earth's history.

Seismic Confirmation of a Thermal Boundary Layer
Our tests predict the likely seismic properties of the manifestation of a TBL that is consistent with other geological and geophysical observations. Here, we interpret existing observations of CMB structure in the context of these results and make inferences about the TBL.
The most comprehensive study of Pdiff waveform variations are from Euler and Wysession (2017), who measured regional variation in slowness anomalies and amplitude decay constant for Pdiff as a function of wave frequency, which is inversely proportional to the thickness of the layer that waves are sensitive to at the CMB. Some regions that are distant from the LLSVPs, such as under northern Asia, displayed increasing slowness with increasing sensitivity close to the CMB on the order of 0.1 s/°, consistent with predictions from our TBL model. Moreover, in the same region the authors recorded a decrease in the amplitude decay constant, that is, less decay with distance, consistent with our waveform simulation showing greater amplitudes as a function of distance relative to the reference. While the study of Euler and Wysession (2017) lacks the spatial coverage and frequency range to be directly compared against our predictions, the observations of the patterns of core diffracted waves are encouragingly similar to the manifestation predicted for our TBL model. Euler and Wysession (2017) also investigated the variation of slowness with period for several models of lowermost mantle structure. We compare the predictions from our synthetics with these models (Figure 14). While separated by vertical offsets, the pattern of our TBL model best matches model RJK2705 of Valenzuela and Wysession (1998) derived from Pdiff observations which shows a reduction in P-wave velocity toward the CMB. Our TBL model also has a similar slowness pattern to that of the two-layer mantle convection model of Solomatov and Moresi (2002). This figure demonstrates that the effects of the TBL could be observable using slowness perturbations of Pdiff waves, although present studies are not yet sufficient.
Observations of the center frequencies of Stoneley modes give some resolution on the structure of the region immediately above the CMB (P. Koelemeijer et al., 2013). While predictions for our TBL model do not fully match observations (Figure 10e), they show better agreement than PREM, which does not have a region of strongly reduced velocity close to the CMB. This suggests that a reduced velocity region in the lowermost mantle is necessary to fitting the Stoneley mode observations.
By, combining multiple disciplines, we significantly improve compared with previous estimates of Q CMB . Our model and preferred Q CMB has important implications for the Earth, most notably the degree of partial melt in the lowermost mantle. Combining recent estimates of lower mantle solidus temperatures with estimates of core temperature, we may expect partial melting anywhere from 0 to 15 km, and 8-26 km above the CMB depending on core temperature, both of which are incompatible with seismic observations.
An important constraint on the maximum temperature of the TBL comes from the combination of petrology with seismology. As discussed above (Section 7.1), many thermal evolution models predict temperatures greater than the experimentally measured solidi at CMB pressures for likely mantle compositions: ∼4180 K for CI chondrite and peridotite-type material Fiquet et al., 2010), and 3570 ± 200 K for pyrolite material (Nomura et al., 2014); we may expect partial melting anywhere from 0 to 15 km, and 8-26 km above the CMB depending on core temperature. Seismic observations preclude a global melt layer at the present day. While ULVZs are often interpreted as partial melt, these are intermittent features of the CMB (McNamara et al., 2010;Rost et al., 2010), implying that the lower mantle does not contain a global melt layer at the CMB. Whereas, in other areas of the mantle, the core-mantle transition has been resolved to be less than 1 km thick (Vidale & Benz, 1992), again ruling out a global melt layer. Studies of reflected phases, which have the best vertical resolution of all methods used to study ULVZ, resolve the minimum resolvable thickness of ULVZs to be 7 km using ScP reflections . If the temperature of the TBL exceeded the solidus of lower mantle mineralogy, it would have to do so for less than 7 km. This implies that a pyrolite lower mantle would be hard to reconcile with seismologic data given that for all three of our TBL models (low, moderate, and high core temperatures) we would expect a global melt layer that is thick enough to be seismically visible (Figure 4b), which is not observed. A pyrolite lower mantle would only be consistent with seismologic data if the solidus temperature falls toward the very top of the possible range from Nomura et al. (2014), and the core temperature falls at the bottom of the range. In contrast, a CI chondrite or peridotite lower mantle would likely only generate a global melt layer for the hottest core temperature case Fiquet et al., 2010). We thus place constraints on the CMB temperature, depending on mantle composition. Figure 14. Pdiff slowness as a function of dominant period for filtered synthetic data. Measurements from our synthetics are shown by solid lines with circles to indicate the period at which the slowness was measured. These models are the 1D model ak135 (Kennett et al., 1995) (black), and ak135 plus our thermal boundary layer (TBL) model (red). Predictions from additional models of are shown by dashed lines. These models are: preliminary reference Earth model (PREM; green), PREM with an embedded ULVZ (purple), PREM with a positive gradient in the lowermost 200 km (gray), a second cell of convection in the lowermost 200 km (blue), and RJK2705 of (Valenzuela & Wysession, 1998) which has a strong decrease in velocity toward the coremantle boundary in the lowermost 150 km, most similar to our TBL model. These models are copied from Euler and Wysession (2017) who use a different period range than that used here.
A further effect of increased temperature in the lower mantle would be on seismic anelasticity. However, there is disagreement over the anelastic effects of temperature anomalies on seismic velocities (Trampert et al., 2001). Though it is argued that increased temperature could cause reduced velocities due to anelastic contributions (S. I. Karato, 1993), there are uncertainties on the quantitative effects due to limited mineralogical data (J. P. Brodholt et al., 2007;Matas & Bukowinski, 2007). This prevents directly interpreting velocity anomalies close to the CMB in terms of the anelastic effects of temperature.

Constraints on Outer Core Stratification
There is observational evidence for a stably stratified layer below the CMB from seismic and geomagnetic data (e.g., Braginsky, 1993;Young & Lay, 1987), which would modulate heat flow out of the core, thus is of importance here. Core stratification could develop in several ways. A compositionally stratified layer at the top of the core could result from the barodiffusion  or immiscibility of light alloying elements (e.g., O, S, and Si; Arveson et al., 2019), or from the diffusion of light elements from the mantle into the core (B. A. Buffett, 2010;B. A. Buffett & Seagle, 2010). It is also possible that the heat flux out of Earth's core is subadiabatic (see Section 6.1). A subadiabatic heat flux will result in a thermally stratified layer at the top of the core because heat conducted along the adiabat toward the CMB exceeds the heat removed by the mantle (i.e., Q CMB ≤ Q ad ). Warm, buoyant fluid accumulates at the top of the core unless strong chemically driven convection can mix this layer back into the interior of the core. Thermal stratification is strongly dependent on the thermal conductivity of the core, which is not well constrained. Some theoretical studies suggest the thermal conductivity of core material is quite high (up to 150 W/mK), while other models and some experimental results are lower (35-90 W/mK;de Koker et al., 2012;Konôpková et al., 2016;Ohta et al., 2016;Pozzo et al., 2012Pozzo et al., , 2013Xu et al., 2018). Thermal and compositional stratification are, of course, not mutually exclusive, and light elements may pool within a thermally stratified layer. These scenarios produce testable predictions about the properties and observability of such a layer.
The presence, strength, and thickness of a stratified layer below the CMB can be estimated by its ability to propagate waves, such as Alfvén, inertial, MAC (magnetic, Archimedes and Coriolis), and seismic waves. Observable evidence of such a layer can be sought in differential travel times of SmKS seismic waves, variations in the lengthof-day, and geomagnetic secular variation. The strength of stratification (i.e., the density gradient relative to the adiabatic gradient) is expressed in terms of the buoyancy frequency, which is the frequency that a parcel of fluid will oscillate in a density stratification without the effects of rotation or magnetic field. Compositional stratification is thought to be stiffer than thermal stratification, that is, it has a higher N. For example, a thermodynamic model of chemical interactions with the mantle predicts a compositionally stratified layer of light elements at the top of the core 60-70 km thick with N ≈ 55 Ω, where Ω = 0.729 × 10 −4 s −1 is the angular velocity of Earth's rotation (B. A. Buffett & Seagle, 2010). Such a layer would produce a 2% increase in seismic velocity, while seismic data are interpreted to show a decrease in seismic velocity of 2% (e.g., Helffrich & Kaneshima, 2010).
Alternatively, Lister and Buffett (1998) model the competition between vigorous chemical convection and subadiabatic thermal conditions and predict a thermally stratified layer of ∼100 km thickness with a N less than 10 −4 s −1 = 1.4 Ω. Such a layer may be geomagnetically observable: B. Buffett (2014) showed geomagnetic dipole field fluctuations and zonal core flow at the top of the core can be described with a model of MAC waves propagating in a stratified layer 140 km thick with a buoyancy frequency ∼1 Ω. This layer would have a temperature gradient of 0.12 K/km yielding a density gradient of 0.012 kg/km (B. Buffett, 2014), so it may not be seismically observable due to small temperature differences.
Dynamo simulations predict that the presence and strength of a stratified layer at the top of the outer core modulates how lateral heterogeneities in Q CMB impact outer core convection. With strong stratification (i.e., large N), Christensen (2018) found circulation at the very top of the stratified layer due to thermal heterogeneity, but the rest of the core convection is unaffected by lateral heterogeneity in Q CMB . On the other hand, Olson et al. (2017) modeled weaker stratification and found that convection in the lower portion of the outer core is affected by boundary heterogeneities. Similarly, the "locked" dynamos of Sreenivasan and Gubbins (2008) have convection and magnetic field morphology locked by the local CMB heat flux pattern (q CMB ), and produce fields which are consistent with high latitude magnetic flux patch concentrations that seem quasi-localized over North America and Russia in a 10 ka field model CALS10k.1b (Korte et al., 2011). The rapidly rotating convection simulations of J. E. Mound and Davies (2017) and J. Mound et al. (2019) found localized stagnant regions below LLSVPs rather than a global stratified layer.
It has been argued that the high-latitude magnetic flux patches observed in Holocene field models require radial core flow near the CMB to form, which would be suppressed by a thick and strong global stratification (e.g., Amit, 2014;Gubbins et al., 2007;Huguet et al., 2018;Lesur et al., 2015;Whaler, 1980). These studies use the secular variation of the geomagnetic field to invert for fluid velocity just below the core mantle boundary by assuming there is negligible ohmic diffusion (i.e., the frozen flux approximation). However, Barrois et al. (2017) argue that diffusion is not negligible and some CMB magnetic flux patches are due to diffufsion rather than advection.
A global stratified layer is consistent with observations of MAC wave propagation (e.g., B. Buffett, 2014;B. Buffett et al., 2016). MAC waves (both zonal waves with periods of approximately 30-115 years and equatorially trapped waves with periods ∼5 years) have been observed in geomagnetic secular variation, dipole field fluctuations, and variations in length-of-day. The observations are best fit by a stratified layer with a thickness in the range 20-150 km. This stratified layer model makes several strong predictions: a Nu o = 0.88-0.98, which place a bound on Q CMB = 13.3-14.8 TW assuming Q ad = 15 TW. This range of Nu o is shown in Figure 12 with a yellow area and square. It also predicts N = 0.84-5 Ω, which indicates the stratification is likely thermal or weakly chemical (B. Buffett, 2014;B. Buffett & Matsui, 2019;B. Buffett et al., 2016). Again, the observations of MAC waves do not distinguish thermal and compositional stratification, but the thickness and buoyancy frequency predicted are consistent with thermal stratification (B. Buffett et al., 2016) which indicates Q CMB > Q ad . This is consistent with the lower values of our TBL model (13.6-15 TW).
While many seismic studies report the existence of a layer of reduced seismic velocity in the upper outer core, at present there is no clear way to reconcile the observations with geodynamic and mineralogical constraints. Geodynamo studies matching geomagnetic observations limit thickness of a stratified layer to less than 400 km thick ( Figure 12). Meanwhile, seismic studies that require a layer of reduced velocity in the upper outer core disagree over the thickness ranging from between 90 km (Tanaka, 2007) to 450 km thick (Kaneshima, 2018), and other studies arguing for no velocity reduction (Alexandrakis & Eaton, 2010;Irving et al., 2018). Moreover, the models of the development of a stratified outer core layer from a concentration of light elements which are consistent with mineral physics, cosmochemistry, and geodynamics (Badro et al., 2014;B. A. Buffett & Seagle, 2010;Umemoto & Hirose, 2015) predict a layer of increased seismic velocity and reduced density. The SmKS waves have little sensitivity to density, thus density cannot be used as a diagnostic parameter. It has been argued that nonideal mixing of a light element could result in a seismic velocity reduction (Helffrich, 2012), although subsequent studies have not supported this finding (Badro et al., 2014;Umemoto & Hirose, 2015), although it has been suggested that elemental exchange could generate a low density layer with a seismic velocity reduction, although these scenarios may be unrealistic (J. Brodholt & Badro, 2017). As such, the structures required to fit seismology are in conflict with those which fit other measures. Uncertainties in the measurement errors on SmKS data (Souriau et al., 2003), the contaminating effects of mantle structure (Garnero & Helmberger, 1995;Tanaka, 2014), and the disagreement between seismic studies of the same layer do cast some doubt on the robustness of the seismic models. Our conclusions about the TBL still stand regardless of the cause of the outer core stratification, thus we do not conclude either way about the existence and properties of this layer.

Section Summary
In this section we integrated constraints from seismic and geomagnetic observations, in the context of Q CMB . We used geomagnetic observations and thermodynamic theory to place self-consistent constrains on Q CMB at both the present day and back in time. We seismically test the structures in the lowermost mantle and uppermost outer core. Finally, we examined the possibility of a stratified layer at the top of the outer core and the predicted magnetic and seismic signals such a layer would produce.

Summary and Conclusions
Here we present constraints on the heat flux across the core mantle boundary from multiple disciplines. We assess petrological and mineral physics constraints on the temperature of the CMB from the mantle side and the core side. We then use a simple analytical description of the TBL model to describe the temperature profile across the lower mantle to connect estimates of the temperature at the CMB with mantle adiabats. We test the suitability of our simple analytical model with a numerical simulation to describe CMB heat flux under more complex conditions, such as compositional heterogeneity and temperature dependent viscosity. Using these temperature profiles, we determine that the TBL is likely below our current seismic detection thresholds. Considering the core, we then investigate what CMB heat flux is required to explain the existence and persistence of the magnetic field and the possible existence of a stratified layer atop the outer core. Overall, our simple TBL model reconciles disparate constraints on the global heat flow across the core mantle boundary at the present day from petrology, mineral physics, geodynamics, seismology, and geomagnetism ( Figure 15) resulting in a most likely range for Q CMB of 13-15 TW. Though we apply several simplifying assumptions, our multidisciplinary approach significantly narrows the range for possible core heat flux at the present day. Moreover, we can predict that Q CMB was likely higher in the past to account for the observed magnetic field. Figure 15. Summary of the constraints on present day Q CMB from the various fields showing the ranges of constraints discussed in the text. Dashed region indicates range of agreement between constraints.