Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle

The Earth's magnetic field is generated by a dynamo in the outer core and is crucial for shielding our planet from harmful radiation. Despite the established importance of the core-mantle boundary heat flux as driver for the dynamo, open questions remain about how heat flux heterogeneities affect the magnetic field. Here, we explore the distribution of core-mantle boundary heat flux on Earth and its changes over time using compressible global 3-D mantle convection models in the geodynamic modeling software ASPECT. We discuss the use of the consistent boundary flux method as a tool to more accurately compute boundary heat fluxes in finite element simulations and the workflow to provide the computed heat flux patterns as boundary conditions in geodynamo simulations. Our models use a plate reconstruction throughout the last 1 billion years -- encompassing the complete supercontinent cycle -- to determine the location and sinking speed of subducted plates. The results show how mantle upwellings and downwellings create localized heat flux anomalies at the core-mantle boundary that can vary drastically over Earth's history and depend on the properties and evolution of the lowermost mantle as well as the surface subduction zone configuration. The distribution of hot and cold structures at the core-mantle boundary changes throughout the supercontinent cycle in terms of location, shape and number, indicating that these structures fluctuate and might have looked very differently in Earth's past. Our results have implications for understanding the Earth's thermal evolution and the stability of its magnetic field over geological timescales. They provide insights into the potential effects of the mantle on the magnetic field and pave the way for further exploring questions about the nucleation of the inner core and the past state of the lowermost mantle.


INTRODUCTION
The heat flux out of the Earth's outer core into the overlying mantle is one of the drivers for the geodynamo responsible for generating Earth's magnetic field.Past modeling studies have demonstrated that both the magnitude and spatial variability of this heat flux have a strong impact on the convection patterns in the outer core and the resulting magnetic field.Heterogeneous mantle forcing can organize flow near the top of the core (Mound et al. 2019), introducing non-zonal structure into the time-averaged magnetic field and giving it a morphology and secular variation matching Earth's modern field (Mound & Davies 2023).Geodynamo models also suggest that the amplitude and pattern of the heat flux heterogeneity across the core-mantle boundary (CMB), particularly near the equator, affect the average timing between polarity reversals of the magnetic field (Glatzmaier et al. 1999;Olson et al. 2010;Olson & Amit 2014), the deviations of the time-averaged field from a geocentric axial dipole (Bloxham 2000b;Heimpel & Evans 2013), the polarity transition paths during reversals Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 3 and excursions (Kutzner & Christensen 2004), the field strength (Takahashi et al. 2008), and the location of intense geomagnetic flux patches (Olson & Christensen 2002;Gubbins et al. 2007;Sahoo & Sreenivasan 2020).For example, simulations show that locations of minimum field intensity are clearly correlated with above average heat flux patterns at the CMB (Korte et al. 2022).Furthermore, core-mantle interactions are a major influence on the secular variation of the magnetic field computed in geodynamo simulations (Bloxham 2000a).
In addition, paleomagnetic data reveal variations in the magnetic field on much longer time scales (> 1 Myr) than characteristic for circulation in the outer core (centuries), instead matching typical time scales of mantle convection.For example, there is evidence that the frequency of magnetic field reversals varies in a periodic manner, with Phanerozoic superchrons (> 10 Myr periods with few to no magnetic reversals) occurring roughly every 200 Myr (Biggin et al. 2012), often being preceded by extended periods of hyper-reversal (∼10 reversals/Myr; Biggin et al. 2012;Meert et al. 2016).Olson & Amit (2015) find an inverse correlation between the reversal frequency and the activity of Large Igneous Provinces for the past 160 Myrs of Earth's history, suggesting a link between plume activity and the geodynamo.The changes in reversal frequency over time are suggested to correlate with changes in magnetic field strength, with lower field values occurring during periods of high reversal frequency, representing a highly unstable state of the magnetic field, and stronger values occurring during extended periods of non-reversal, representing a more stable state (Cox 1968).This apparent coupling between Earth's mantle and core is also evidenced by the correlation of the magnetic field calculated from geodynamo models using the lowermost mantle seismic velocity patterns as boundary condition with the present-day geomagnetic field (Gubbins et al. 2007).
These observations and modeling results suggest that the influence of the mantle is reflected in the paleomagnetic record.However, understanding this relationship remains difficult for several reasons.
Although we have reasonable estimates of the present-day heat flux at the core-mantle boundary and seismic tomography models offer insights into the distribution of hot and cold material along this boundary, heat flux patterns have likely undergone substantial changes throughout Earth's history.
Unfortunately, we lack a direct record of these past variations.Therefore, the magnitude of magnetic field variations caused by mantle convection is currently not well constrained, posing challenges for interpreting the paleomagnetic record.In addition, most dynamo simulations of Earth's outer core that have taken into account a heterogeneous heat flux across the CMB have either used simplified patterns based on low degree and order (usually ≤2) spherical harmonics (i.e., Bloxham 2000bBloxham , 2002;;Olson & Christensen 2002;Kutzner & Christensen 2004;Aubert et al. 2007;Takahashi et al. 2008;Sreenivasan 2009;Olson et al. 2010;Heimpel & Evans 2013;Olson & Amit 2014;Hori et al. 2014; have used seismic tomography (e.g., Masters et al. 1996) to infer the heat flux (i.e., Glatzmaier et al. 1999;Bloxham 2000a;Olson & Christensen 2002;Christensen & Olson 2003;Kutzner & Christensen 2004;Aubert et al. 2007;Gubbins et al. 2007;Aubert et al. 2008;Amit & Choblet 2009;Olson et al. 2010;Gubbins et al. 2011;Sreenivasan & Gubbins 2011;Olson & Amit 2014;Mound & Davies 2017;Christensen 2018;Terra-Nova et al. 2019;Sahoo & Sreenivasan 2020;Mound & Davies 2023), (for a review, see Amit et al. 2015a).Both types of patterns only represent the large-scale variations in heat flux, but do not capture smaller-scale variations or strong lateral gradients.While imposed heat flux heterogeneities proportional to seismic wave-speed anomalies better represent Earth, they still neglect non-thermal sources likely to contribute to the tomographic pattern.Efforts have been made to address this challenge (Amit & Choblet 2009, 2012;Amit et al. 2015b;Choblet et al. 2023), and, as an alternative to seismic tomography, geodynamo models have employed CMB heat flux patterns from a the present-day state of a mantle convection model to capture the relevant physical processes (Olson et al. 2015).Nevertheless, even for the present-day, some of the complexities of lower mantle dynamics are not captured in the heat flux patterns used in most geodynamo simulations.
To our knowledge there has only been one study that has incorporated CMB heat flux patterns corresponding to different times in Earth's history: Olson et al. (2013) impose a lower mantle history based on time-dependent convection going back to 330 Ma as boundary condition to their geodynamo.They find that these models more readily explain the slow variations in reversal frequency in the Phanerozoic Geomagnetic Polarity Time Scale than models with a heterogeneity pattern that does not change.However, they include only the largest scale components of thermal core-mantle interaction by truncating the CMB heat flux pattern at spherical harmonic degree 4 and therefore ignore smaller scale thermal perturbations such as plume formation.In addition, they rely on the heat flux of one specific mantle model (Zhang & Zhong 2011) and therefore can not take into account the uncertainty of the CMB heat flux related to uncertainties in lowermost mantle material properties and chemical composition.
A deeper understanding of the changes in heat flux patterns at the core-mantle boundary over time and their influence on the geodynamo is essential to determine which observed magnetic field changes can be attributed to mantle convection and which require alternative mechanisms, such as the nucleation of the Earth's inner core.To constrain these variations, we can use models of mantle convection as a link between available surface observations and the CMB heat flux pattern.Subducted slabs sinking down towards the base of the mantle and forming cold areas above the CMB cause a large heat flux out of the core.In addition, they push hotter material together in the regions between downwellings, creating hot regions with a lower CMB heat flux.We therefore expect characteristic changes in CMB heat flux patterns throughout the supercontinent cycle.Reconstructions of the motion Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 5 of the tectonic plates at the Earth's surface and their subduction history throughout this cycle can be useful tools to constrain these patterns and their temporal variations.
1.1 Temporal and spatial variations in core-mantle boundary heat flux Global 3D mantle convection models have been used for several decades to provide insights into spatial and temporal variations of the CMB heat flux.These simulations offer a view of how the temperature distribution in the Earth's mantle evolves over time, and when Earth's plate motion history is imposed as boundary condition, they can successfully reproduce the locations of hot regions below Africa and the Pacific inferred from seismic tomography for the present-day Earth (the Large Low Shear Velocity Provinces, or LLSVPs) (McNamara & Zhong 2005;Bull et al. 2009;Zhang et al. 2010;Davies et al. 2012;Bower et al. 2013;Bull et al. 2014;Flament et al. 2017;Cao et al. 2021b;Flament et al. 2022).
This result illustrates how subduction history controls the lowermost mantle structure and suggests that the present-day dominant degree-2 lower mantle structure is not a stable feature throughout Earth's history (Zhang et al. 2010;Zhong & Rudolph 2015;Müller et al. 2022) with both the number and location of hot structures likely having changed over time (Zhong & Rudolph 2015;Flament et al. 2022).While this topic is still debated (Zhong & Liu 2016), with some studies advocating for the long-term stability of the present-day degree-2 configuration (Bull et al. 2014;Cao et al. 2021b), there is agreement that the lowermost mantle structure strongly depends on plate motion history.
Prior studies have defined criteria that mantle convection models used to provide heat flux patterns for geodynamo simulations should fulfill (Olson 2016): (1) a realistic equation of state, considering phase transitions, mantle heterogeneities and initial and boundary conditions for temperature; (2) sufficiently complex rheologies considering pressure-, temperature-and strain rate dependent viscosity and possibly compositional variability; and (3) upper surface velocity boundary conditions derived from plate reconstructions over sufficiently long time spans.However, mantle convection models computing CMB heat flux patterns have either used simple equations of state in incompressible modelsneglecting important contributions to mass and energy transport and violating point (1), or they have studied CMB heat flux patterns arising in convection models without prescribed plate motions, so that plate-like behavior can emerge self-consistently, but generally does not resemble plate motions in the Earth's past-violating point (3).Material properties controlling the transport of heat out of the core that have been simplified in models with prescribed plate motions also include lowermost mantle viscosity and thermal conductivity.In particular the presence of a weak post-perovskite phase at the base of the mantle has been shown to have a strong impact on both pattern and magnitude of CMB heat flux (i.e., Nakagawa & Tackley 2011).The thermal conductivity strongly varies with both pressure and temperature in the Earth's mantle, significantly affecting heat flux patterns and the differences between surface and CMB heat flux (Tosi et al. 2013).Specifically, both mantle compressibility and the depth-dependence of the thermal conductivity act to decrease the volume of subducted slabs in the lowermost mantle: The lowered conductivity near the Earth's surface leads to slower cooling and thinner subducting plates, and the increasing density with depth causes slabs to decrease in volume as they sink.The slab volume is therefore overestimated in models that do not take into account these effects.
Neither of these complexities has been included in global 3D convection models with prescribed plate motions (e.g., McNamara & Zhong 2005;Zhang & Zhong 2011;Bower et al. 2013;Bull et al. 2014;Zhong & Rudolph 2015;Hassan et al. 2015;Flament et al. 2017;Cao et al. 2021b,a;Flament et al. 2022;Müller et al. 2022;MacLeod et al. 2023;Frasson et al. 2023), which instead assume a constant diffusivity.Finally, plate reconstructions that cover the whole supercontinent cycle have only become available recently (Merdith et al. 2021).Even though these models can suffer from large uncertainties related to ambiguous paleomagnetic data, especially in the Precambrian, and uncertain subduction zone locations before 150 Ma, they are opening up the opportunity to study the associated changes in CMB heat flux patterns over sufficiently long time spans.
We here present global 3D mantle convection simulations that fulfil the criteria given above.These compressible, multi-phase, thermo-chemical convection models apply a plate reconstruction throughout the last 1 billion years (Merdith et al. 2021) as surface boundary condition to constrain the spatial and temporal heat flux variations at the core-mantle boundary and how they are affected by the physical properties of the mantle.We discuss how our models can be used to prescribe this heat flux as a boundary condition for geodynamo simulations, and make this workflow freely available, together with our model outputs.Our work provides a tool for the geodynamo community to estimate the variability of the core-mantle boundary heat flux and to use these patterns in geodynamo simulations in the future.

METHODS
We set up our mantle convection models using the geodynamic modeling software ASPECT (Kronbichler et al. 2012;Heister et al. 2017;Gassmöller et al. 2018;Clevenger & Heister 2021).Specifically, we solve the equations for compressible mantle convection, using an equation of state that is based on mineral physics data (see Section 2.1) and an Earth-like rheology (see Section 2.2).Since the focus of our study is deformation in the Earth's mantle, we consider only viscous stresses, and we assume that the viscosity is isotropic and that we can neglect terms including the bulk viscosity (Schubert et al. 2001).This leads to the following equations for conservation of mass, momentum and energy: where u is the velocity, ε is the deviatoric strain rate, p the pressure and T the temperature.
Additionally, η is the viscosity, ρ is the density, g is the gravity vector, C p is the specific heat capacity of the material, k is the thermal conductivity, Q is the intrinsic specific heat production, and α is the thermal expansion coefficient.
Our model geometry is a 3D spherical shell encompassing the whole mantle, and we prescribe the velocity at the surface based on a plate reconstruction of the last billion years (Merdith et al. 2021, see Section 2.3).During the model runtime we compute the heat flux at the core-mantle boundary using an accurate consistent boundary flux method described in Section 2.6.To constrain the uncertain- ties in spatial heat flux variations, we compute several simulations with different material properties (Table 2), as described below.

Equation of State
We treat the mantle as a mechanical mixture of mid-ocean ridge basalt (MORB) and harzburgite, using the composition from Xu et al. (2008).To compute the density, thermal expansivity and specific heat of the different lithologies, we use the thermodynamic modeling software Perple X (Connolly 2009) together with a thermodynamic database (Stixrude & Lithgow-Bertelloni 2011).This approach automatically includes the effect of both compositional variations and mineral phase transitions on buoyancy, heat transport and volume changes/compressibility as described in previous geodynamic modeling studies (Nakagawa et al. 2009).The resulting density for basalt and harzburgite and the density difference between the two is shown in Figure 1.To compute the material properties of the mechanical mixture, we arithmetically average the material properties (based on a composition's volume fraction for density and thermal expansivity, and based on mass fractions for the specific heat).
We use a constant thermal conductivity of 4.7 W/m/K in our most simple setup.Since the thermal  2.
conductivity is expected to be much larger near the core-mantle boundary-increasing the CMB heat flux-we additionally investigate the effect of a pressure-and temperature-dependent formulation (see Table 2).Specifically, we use the thermal conductivity model from Tosi et al. (2013) above 660 km depth and the model from Stackhouse et al. (2015) below (Figure 2b).

Rheology
Our mantle rheology is depth-and temperature-dependent using the preferred viscosity profile of Steinberger & Calderwood (2006) and an Arrhenius law that separates radial and lateral viscosity variations: Here, η 0 (r) is the viscosity profile describing the depth-dependence, and H(r) is the depth-dependent activation enthalpy (as given in Steinberger & Calderwood 2006) defining the dependence on temperature.T ref is the reference temperature profile (the initial mantle adiabat), ∆T = T − T ref is the deviation from this adiabat, n is the stress exponent, and R = 8.314 J/(mol K) is the gas constant.
Since the focus of our study is the lowermost mantle, we do not take into account the strain rate dependence of the viscosity, but the stress exponent is assumed to be n = 3.5 in the upper mantle and n = 1 in the lower mantle to account for the change in dominant creep mechanism from dislocation to diffusion creep.We limit the viscosity to be between 10 20 and 5 × 10 23 Pa s.In particular the lower limit is chosen to make sure that convective processes do not occur on a smaller length scale than can be resolved by our numerical resolution.
We also include two model setups that test additional rheologic complexities (see Table 2).The first one includes a viscosity reduction within the post-perovskite phase.Since the amount of weak-ening remains uncertain with estimates ranging from 1 to 4 orders of magnitude (Hunt et al. 2009;Ammann et al. 2010;Goryaeva et al. 2016), we here choose a viscosity reduction of a factor of 100.
The other setup includes a composition-dependence of viscosity, increasing it by a factor of 10 in the basaltic material compared to the average mantle.This is motivated by the good fit to present-day mantle structure inferred from seismic tomography achieved in the models by Flament et al. (2022) using this value.All viscosity profiles are shown in Figure 2c.
Since our models do not include a visco-plastic rheology or any other mechanism that would weaken plate boundaries, the prescribed plate motions at the surface (Section 2.3) lead to an unrealistically large amount of friction at plate boundaries.This is in particular the case in subduction zones, where the material is cold and viscous, and causes shear heating to be overpredicted.We therefore limit shear heating in Equation (3) in our models to prevent unrealistically high temperatures associated with the forced surface plate motion.Specifically, we compute if material would deform plastically as estimated by a Drucker-Prager yield criterion with a cohesion of C = 10 MPa and a friction angle of ϕ = 0.085.We then limit the stress being used to compute the shear heating to not be higher than this yield strength σ yield .
Q shear = min(2η ε, σ yield ) : ε with (5) This ensures that average mantle temperatures evolve as expected for the Earth's interior in our reference (No. 1; Thermochemical) model, slowly cooling down over time.

Boundary conditions
To model the changing patterns of subduction throughout the supercontinent cycle, we use a reconstruction of the last 1 billion years of plate motion history (Merdith et al. 2021) and prescribe it as velocity boundary condition at the model surface.As the subducted slabs sink downwards into the lowermost mantle, they form regions of low temperature causing a high heat flux out of the core.
In addition, they push hotter and/or chemically distinct material together into plume clusters or thermochemical piles, which feature low heat flux out of the core.This mechanism provides a coupling between the plate motions at the Earth's surface and the heat flux at the core-mantle boundary.Since the modeled subduction zones are prescribed in the same locations they are thought to have been on Earth, we also expect the pattern of core-mantle boundary heat flux to be representative of Earth's history.We realize that there are significant uncertainties associated with plate reconstructions, especially going further back in time than the oldest ocean floor preserved at present-day.Consequently, some Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 11 subduction zones might be missing or be in the wrong location in our models.We nevertheless think that the use of this reconstruction is justified since the objective of our study is to predict characteristic changes in heat flux patterns throughout the supercontinent cycle rather than to constrain the exact heat flux pattern at the Earth's core-mantle boundary throughout the last billion years.
The bottom boundary of our model is closed, but allows for free slip in the direction tangential to the boundary.The temperature is fixed to 273 K at the top and 3700 K at the bottom boundary.
This core-mantle boundary temperature is well within the range of recent estimates.It can not be higher than the pyrolite solidus, which experimental data constrain to 3570±200 K (Nomura et al.

Initial conditions
Since the thermal and chemical state of the Earth a billion years ago is unknown, we here make the simplest assumption and start with an "empty" mantle (without any plumes or slabs) that is a mechanical mixture of 18% MORB and 82% harzburgite.Specifically, we assume that the initial temperature profile is adiabatic with a potential temperature of 1613 K with additional boundary layers at the top and bottom following a half-space cooling model.The top thermal boundary layer has an age of 70 Myr, resulting in an appropriate temperature profile for oceanic plates.The structure of the bottom boundary layer is different between our models setups (see Table 2).In the model without chemical heterogeneities (No. 2; Thermal model), the bottom boundary layer is assumed to have an age of 50 Myr since the low viscosity allows for the frequent rise of plumes, keeping the layer relatively thin.In all other models, we consider the effect of a purely basaltic layer at the base of the mantle with an initial thickness of 200 km.Since this layer is dense, insulating the mantle and delaying plume formation, the thermal boundary layer is expected to be thicker compared to the purely Thermal model.
We therefore set the initial temperature to be 800 K above the mantle adiabat in the lowermost 150 km of the mantle (approximately 3300 K), and compute the temperature above that depth based on a halfspace cooling model with an age of 150 Myr.The resulting temperature profile is shown in Figure 2a.
Note that since the thickness of the boundary layers is defined by age, it is different in the models with a pressure-and temperature-dependent thermal conductivity compared to the ones with a constant conductivity.

Numerical Methods
ASPECT uses an adaptive finite-element mesh to discretize the model domain, which results in a cell size in our models between 45 and 250 km (depending on the refinement level and distance from the core-mantle boundary).Specifically, we refine the mesh in regions where the gradients in temperature, composition and viscosity are high, and we additionally enforce that the mesh is always refined to the highest resolution in the thermal boundary layers (below 2500 km depth and above 80 km depth).
Since we use second-order finite elements, we achieve a resolution in terms of distance between the quadrature points of 39 km in horizontal and 22 km in vertical direction at the base of the mantle, where the core-mantle boundary heat flux is computed.
To minimize numerical diffusion, we use ∼90 million particles (using the implementation of Gassmöller et al. 2018) to track the evolution of the chemical composition.Specifically, our model contains two distinct compositions, harzburgite and mid-ocean ridge basalt.We assign each particle a composition that represents the fraction of basalt according to our initial conditions and interpolate particle properties to the finite-element mesh using a quadratic least-squares approximation.
Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 13 reported in the publications (gray lines).Note that the right panel only shows the case of Ra = 10 5 and the anelastic liquid approximation for different dissipation numbers to improve the visibility of the figure and that we used the case "RefVT" as reference value.

Computing heat flux with the consistent boundary flux method
The main purpose of our models is to compute the spatially and temporally variable heat flux density out of the liquid outer core and into the overlying mantle.However, accurately computing derived quantities like heat flux at the boundary of a numerical domain is challenging, because many numerical methods emphasize accuracy as an integrated quantity over the volume of a cell, and lose accuracy towards the faces and edges of a cell.In addition, the accuracy of a solution derivative is always reduced compared to the solution quantity itself.To circumvent these limitations, instead of computing the heat flux directly from the gradient of the temperature solution using Fourier's law, we utilize a consistent boundary flux (CBF) method as described for the heat equation in Gresho et al. (1987).CBF has been benchmarked in the geodynamics community as a very accurate method to compute another derived quantity-dynamic topography-and has also been suggested as a promising technique to compute heat flux before (Zhong et al. 1993;Weir 2019).
We have benchmarked our CBF implementation using the incompressible models described in Blankenbach et al. ( 1989) and the compressible models of King et al. (2010), and illustrate the accuracy im-provement of the CBF method in Figure 4.A description of the method and full benchmark results including a convergence analysis are provided in Appendix C and all data to reproduce the benchmarks are included in ASPECT.In Figure 4 we plot the heat flux postprocessing results for identical models computed with the CBF and traditional gradient based computations and compare them to benchmark results from the literature given above.It is very clear that while both methods converge to the reference results at high resolutions, the CBF method does so much faster and at much coarser resolutions (we show in the supplementary information that CBF's convergence order is 1-3 orders higher than the gradient based method).Additionally, the convergence of the CBF method is more consistent in the sense that it mostly approaches the reference value from one direction, while the gradient-based method consistently tends to underestimate heat flux at coarse resolutions and overestimate heat flux at intermediate resolutions.This behavior makes extrapolation of under-resolved model results-a common challenge in numerical geodynamics-much less accurate.Therefore, CBF heat flux represents a significant improvement in the accuracy of heat flux computations in geodynamic modeling studies.

RESULTS
Our models compute the heat flux distribution across the CMB throughout the last 1 billion years of Earth's history.We will first compare the evolution of average properties (heat flux, temperature) and the present-day state of our models to available observations to show that our models are Earth-like.
In a second step, we will analyze the computed heat flux patterns and constrain the possible spatial and temporal variations.

Thermal evolution and average CMB heat flux
In our reference (Thermochemical, Table 2) model setup, both the thermal evolution and the heat flux across the core-mantle boundary are consistent with observations.The average temperature (Figure 5a) conductivity reduces heat loss at the Earth's surface (due to the lower conductivity at low pressures, see Figure 2) and increases heat flow across the CMB (due to the higher conductivity at high pressures), making mantle cooling less efficient.The presence of a weak post-perovskite layer enhances convective heat transport away from the CMB, resulting in an increased heat flux.On the other hand, more rigid piles do not appear to substantially affect the amplitude of the CMB heat flow.
In our computations, all models featuring a p-T-dependent conductivity-which is likely the better approximation of Earth's mantle-display a thermal evolution inconsistent with observations (Figure 5a).This suggests that our models underestimate the heat flux at the Earth's surface.One possible reason is that the data the plate reconstruction is based on becomes more sparse going back in time, and therefore some-especially intraoceanic-subduction zones may be missing in our boundary conditions.In addition, our models do not account for mechanisms such as melt extraction that facilitate more efficient heat transport across the lithosphere.However, since our focus is primarily on the general patterns of heat transport throughout the supercontinent cycle and their changes over time, small variations in average or upper mantle temperature do not substantially impact our results.Hence, we consider our models to be a reasonable approximation of the lowermost mantle processes under investigation.
Note that the heat flux density (Figure 8) goes through a phase of unrealistically strong variations in the first ≈200 Myrs of model evolution before featuring smaller variations around a quasi-steadystate.This effect is caused by the initialization of the models, which initially do not feature subducted slabs or rising plumes.Therefore, the heat flux first strongly decreases to lower values than expected on Earth as the core-mantle boundary region heats up, growing a thick thermal boundary layer that insulates the core.When the first cold slabs reach the CMB, they trigger the ascent of the first plumes and abruptly increase the amplitude of the heat flux.After this "spin-up" phase of approximately 200 Myrs, the modeled heat flux stabilizes and its variations reflect changes in the lower mantle structure due to the subduction history.We therefore do not interpret these first 200 Myrs of model evolution.We here present only a qualitative comparison, since our objective is not to to achieve the best possible fit of lowermost mantle structure to present-day Earth but to show that our models are characteristic for Earth's changing CMB heat flux patterns in general.In the Thermal model, seismically slow regions are not matched well by regions of high temperature, with hot plumes being thin and roughly evenly spaced along the core-mantle boundary.However, all of our models that include an intrinsically dense layer show high-temperature anomalies at the base of the mantle in roughly the same regions as in the seismic tomography, i.e. below Africa and the Pacific (see Figure 6).Similarly, areas of low temperatures are located below the ring of subduction surrounding the Pacific, both in the geodynamic models and seismic tomography.The exact geometry of the hot and cold regions varies between the different models and also differs from the shape of fast and slow anomalies in the tomography.Specifically, the hot regions in the geodynamic models are larger than the seismically slow regions in the tomography and in some models are broken up by colder patches that do not have seismically fast counterparts.However, the shape of the LLSVPs is reproduced fairly well in the p-Tdependent model, with the African LLSVP being elongated in North-South direction and the Pacific LLSVP in East-West direction.
We emphasize that a more quantitative comparison would require a conversion from temperature Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 17 and composition to seismic velocities and applying a tomographic filter to make the amplitudes and gradients of the anomalies comparable.But even our more qualitative approach already shows that the match between geodynamic model and tomography is not as good as in some other recent models (i.e., Flament et al. 2017) whose plate motion history did not go as far back in time.This suggests that some of the older subduction zones (>250 Ma) that are less well-constrained in plate reconstructions deviate from their actual past locations in the Earth and that this old, cold material still affects the present-day lowermost mantle structure (i.e., visible as cold patches beneath Africa in the Thermochemical and Strong basalt models).
Seismic tomography models also suggest that the LLSVPs are not flat features, but extend several hundred kilometers above the CMB (e.g., Cottaar & Lekic (2016) estimate the vertical extents of different parts of LLSVPs between 300 km and 1200 km).This behavior is not reproduced in our Thermal model, but all models with an intrinsically dense layer feature high-temperature regions (> 750 K hotter than the mantle adiabat) at the base of the mantle that are at least 500 km high.The more effectively cold subducted slabs can displace this hot material, the more topography can be created, i.e., colder/denser/thicker slabs can generate taller piles of hot material.Therefore, the Thermochemical model, which has a higher conductivity at the surface leading to thicker subducted slabs, features the tallest piles (600-800 km, compared to 500-700 km in the p-T-dependent model).In addition, weak thermochemical structures can not maintain a high topography, leading to a relatively low height of piles above the CMB (500-600 km) in the model with weak post-perovskite.However, the model with the lowest height of piles (400-500 km) is the one with a higher viscosity of basaltic material.
Since convection within the strong piles is very sluggish, their tops cool down more efficiently than in the other models and hot regions do not extend as far upwards.Note that all the height values given represent the larger-scale structure of the piles, with individual plumes featuring high temperatures up to much shallower depth.Given the uncertainties of seismic tomography models, any of the models with an intrinsically dense layer could be representative of the Earth's lowermost mantle structure.q * = q max − q min q ave (7) where q max , q min and q ave are the maximum, minimum and spatially averaged heat flux across the CMB at a given point in time.We note that some studies follow a slightly different definition with an additional factor of 2 in the denominator (i.e., Olson & Christensen 2002), so care needs to be taken when comparing absolute values of q * .Figure 8 shows q max , q min and q ave , with the resulting q * for all our models displayed in Figure 9.The minimum heat flux is close to zero in all models, representing areas at the base of the mantle with similar temperatures as the core and therefore a vanishing CMB heat flux.Conversely, the maximum heat flux varies over time much more substantially and is different between the different models.Note that our models include radiogenic and frictional heating so that locally, a small negative heat flux (out of the mantle into the core) is possible.
In the Thermal model, the minimum and maximum heat flux differ from the average by almost the same amount, with only weak variations in the maximum heat flux.This model also features a relatively low area fraction of the CMB with the heat flux being close to zero (i.e., yellow areas in Figure 10).Conversely, in all thermochemical models, the average heat flux is much closer to the minimum than to the maximum heat flux, with the latter also undergoing more substantial changes over time.The underlying cause is the low heat flux at the base of the hot and dense thermochemical structures, which cover a large fraction of the CMB area.The presence of these hot regions reduces the average CMB heat flux compared to an isochemical model that otherwise uses the same parameters, while the minimum and maximum heat flux remain almost unaffected.The heat flux amplitude at the base of these piles is also different between the models.The more efficiently heat is both conducted and convected within the piles, the higher the CMB heat flux at their base.Consequently, the Strong basalt model-where piles are harder to deform and convect more sluggishly-exhibits the lowest CMB heat flux within the piles (Figure 10).A pressure-and temperature-dependent conductivitywhich increases the conductivity at the base of the mantle (see Figure 2)-amplifies the CMB heat sphere in these models leads to a thinner top thermal boundary layer, and therefore to thinner plates.
After these thinner plates are subducted and sink down to the lowermost mantle, they also heat up more rapidly because of the increase in conductivity with pressure.Both effects lead to a decreased volume and a shorter preservation of cold anomalies at the CMB.Therefore, the fraction of CMB area covered by high-heat flux patches is larger in the Thermochemical model-which features a constant conductivity-compared to the other models with chemical heterogeneities (Figure 10).
The Weak post-perovskite model exhibits the most extreme heat flux variations and the highest heat flux values out of all of the models (both locally and globally).This model features cold slabs that are substantially easier to deform compared to the other models, allowing cold material to spread out along the core-mantle boundary more quickly and facilitating the formation of localized cold "puddles" (see Figure B1).Consequently, a strong thermal gradient develops, leading to a large heat flux out of the core.However, this increased heat conduction also causes rapid heating of the slabs, reducing the thermal gradient and resulting in a faster drop in the maximum heat flux.In the Earth, convection would be faster within the thermochemical piles as well, but due to our limit on the minimum viscosity (which is due to our model resolution) this effect is not included in our models.
The different physical behavior of models with different material properties results in different characteristics of the amplitude of heat flux heterogeneities (as expressed by q * in Figure 9).The Thermal model features the lowest value, approximately 2, while the thermochemical models yield q * values ranging from 3 to 4. The only exception is the Weak post-perovskite model where q * values range between 5 and 9.Note that these values usually need to be adjusted in geodynamo simulations by subtracting the adiabatic core heat flux, as discussed in Section 4.
The frequency of temporal variations is primarily influenced by the viscosity of the lowermost mantle.The Weak post-perovskite model exhibits faster changes over time compared to the other models since material can deform more quickly.On the other hand, the heat flux changes more slowly in the Strong basalt model since it takes longer to deform the dense piles at the base of the mantle.2012) for suspected activity up to 600 Ma (lighter shades).Before 600 Ma, geomagnetic polarity is poorly known.Olson et al. 2010;Olson & Amit 2014).However, we see no clear relationship between q * and reversal frequency for any of our models.
Since the equatorial heat flux has been suggested to affect the magnetic field, and in particular its reversal frequency (for example Olson et al. 2010), we have computed the averaged heat flux within ±20 • latitude from the equator (shown in Figure 9, and as thick colored lines partly overlapping with the average heat flux in Figure 8).In all models, the equatorial heat flux and the average heat flux follow similar trends, with both curves being particularly close in the Thermal model (which features the lowest heat flux heterogeneity).When the equatorial heat flux deviates from the average, it generally exhibits lower values because the hot thermochemical structures are located close to the equator throughout most of the last billion years (see Figure 10).The different models do not follow the same evolution (Figure 9), highlighting the dependence of the timing of minima and maxima in both equatorial and average heat flux on lowermost mantle material properties.However, some trends are consistent in all models: There is a maximum in the equatorial heat flux around 600 to 550 Ma, which is related to cold remnants of slabs subducted during Rodinia break-up.This cold material is advected along the CMB into the equatorial region from further South at ∼600 Ma due to the large-scale downwelling in the Southern hemisphere during Gondwana assembly.Another, smaller maximum can be seen in the Thermochemical, p-T-dependent and Strong basalt models at ∼200 Ma.

Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 23
In addition, several models feature an equatorial heat flux minimum around 400 Ma, related to very little cold material being present near the equator (see Figure 10).Earlier modeling (Zhang et al. 2010) has suggested that minima in the equatorial heat flux at 270 and 100 Ma may have been responsible for the Kiaman and Cretaceous Superchrons (which are marked as green bars in Figure 9).These minima are not consistently present in all of our models, with only the Strong basalt model featuring low equatorial heat flux during both periods and none of the models showing a clear relation between equatorial heat flux and the observed reversal frequency.

Temporal changes in core-mantle boundary heat flux patterns
While the amplitude of heat flux heterogeneities is likely to affect the geodynamo, the distribution of these heterogeneities could have an important effect as well.Figure 10   when comparing the regions of high heat flux (blue colors in the first 5 columns in Figure 10) to the trench locations (black lines in the last column) 100-200 Myr earlier (one or two rows further up).For example, the girdle of high heat flux separating the two low heat flux piles between Africa and Pacific at the present-day (last row in Figure 10) corresponds to the ring of subduction zones surrounding Pangea as it was breaking apart (the black lines surrounding the continents in the last column and third-/second-to-last row in Figure 10).Conversely, during the time frame 600-500 Ma, the later stages of the Gondwana assembly, most of the subduction zones were located in one hemisphere (centered around where the Southern Atlantic Ocean is located today).This is reflected in the heat This connection between the supercontinent cycle and CMB heat flux patterns becomes even clearer in a more quantitative analysis of the temporal variations of their characteristic wavelengths using spherical harmonics.To evaluate which spatial pattern is most prevalent, we assess the power spectrum of the different spherical harmonics degrees over time (Figure 11).Degree 1 dominating the spectrum indicates a large-scale heat flux difference between the two hemispheres, one featuring high heat flux, the other featuring low heat flux.Conversely, a high power in degree two corresponds to two large regions of low heat flux separated by a band of low heat flux (or the other way around).If higher degrees are dominant, spatial variations on a smaller scale are more prevalent.
Figure 11 shows how the prevalence of the different degrees changes throughout the supercontinent cycle.At around 800 Ma, and then again around 200-0 Ma, there is a high power in degree 2 in all models, corresponding to two high-temperature (low heat flux) structures in the lowermost mantle that are split by a band or several patches of cold remnants of subducted slabs (high heat flux; see Figure 10).In both instances, this pattern succeeds the start of supercontinent (Rodinia or Pangea) break-up, where subduction zones surround the supercontinent.As this girdle of cold subducted slabs sinks down to the lowermost mantle, it displaces hot material, splitting it into two structures, one beneath the (disassembling) supercontinent, one beneath the superocean.This is also the stage that we observe for the present-day Earth.The later stages of the break-up and the transition to the assembly of the next supercontinent requires subduction zones to be more spread out.The Merdith et al. (2021) plate reconstruction also features a rapid reorganization of plate boundaries during this time frame (around 850-700 Ma).This is reflected in a low power of both degree 1 and 2 as this subducted material reaches the lowermost mantle (around 700-600 Ma), with higher spherical harmonics degrees Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 27 being more dominant and the corresponding smaller-scale alternation between hot and cold material along the CMB.The only time where degree 1 is dominant in some of our models (and the power of degree 2 is low across all models) is following the assembly of a supercontinent (i.e., around 450 to 300 Ma).Because the subduction zones are distributed predominantly in the Southern hemisphere in Merdith et al. (2021) at 600-450 Ma, where Gondwana is being assembled, cold material also reaches the lowermost mantle only in one hemisphere.This leads two one hemisphere with predominantly cold, subducted material (high heat flux) below the supercontinent and one hemisphere with predominantly hot material (low heat flux) below the superocean.Due to the delay of 150-200 Myrs between subduction zones in the plate reconstruction being able to affect CMB heat flux, the prevalence of degree 1 then occurs after the assembly and during the supercontinent stage.Note that the subducted material can easily be pushed along the CMB once it has reached the lowermost mantle because of the negligible friction at the CMB.Because cold material is pushed northwards in our simulations (which is mostly an effect of a net rotation of the whole mantle between 450 and 400 Ma), the location of the high heat flux hemisphere (centered around present-day Africa) does not exactly correspond to the hemisphere where subduction occurred (centered around the South Pole).
The alternation between a degree 1 and 2 lower mantle structure has been debated over the last decades.While some studies argue for the prevalence of a degree-1 structure before Pangea formation (Zhang et al. 2010;Zhong & Rudolph 2015) and the aggregation and dispersal of basal mantle structures over time (Flament et al. 2022), other studies have suggested that lower mantle structure is dominated by spherical harmonic degree 2 most of the time (Cao et al. 2021b) and that the two LLSVPs observed today have remained close to their present-day positions for at least the past 410 Myr (Bull et al. 2014) to be adjusted by subtracting the conductive heat flux along the core adiabat from the CMB heat flux provided by the mantle convection model.This adiabatic core heat flux depends on the thermal conductivity and the material properties that determine the temperature gradient along the core adiabat, so the value that needs to be subtracted will depend on the specific setup chosen for an individual geodynamo simulation.However, we here also want to demonstrate how this conversion would affect the amplitude of spatial heat flux variations as seen in the geodynamo model, i.e., the q * parameter.
We choose three different values of 5, 10, and 15 TW for the adiabatic core heat flux, within the range of recent estimates (2.3 to 16 TW, Pozzo et al. 2012;Nimmo 2015;Davies et al. 2015;Silber et al. 2019;Mound & Davies 2023).In addition, we assume that the q max , q min and q ave cannot change sign from subtracting the adiabatic heat flux, i.e., heat will not flow out of the mantle into the core because the CMB heat flux computed in the mantle model has a smaller amplitude than the conductive heat flux along the core adiabat.
q * = max(q max − q adi , 0) − max(q min − q adi , 0) max(q ave − q adi , 0) with q adi being the adiabatic core heat flux density.We ignore results where the average CMB heat flux is smaller than the adiabatic core heat flux since in these cases, heat loss to the overlying mantle would likely not be able to drive the geodynamo.
The resulting values for q * are shown in Figure 12.In all cases (including the case where q adi is not subtracted), q * ≥ 2, imposing a minimum on the expected amplitude of heat flux heterogeneity at the core-mantle boundary.Since estimates for the adiabatic core heat flux are similar to the estimates for the total heat flux out of the core, it is impossible to provide an upper limit to q * , and the case with the 15 TW adiabatic core heat flux illustrates that q * can easily reach values > 30.If we further assume that the adibatic core heat flux is at least 10 TW (based on the more recent higher estimates for the thermal conductivity of the core) and that there is a layer of dense material at the base of the mantle, then we would expect values of q * ≥ 6.

MODEL LIMITATIONS
We have already outlined some of the limitations of our models above, such as the uncertainty in the plate reconstructions for earlier times in Earth's history and the relatively large total CMB heat flux in some of our models.Another model simplification is that our CMB temperature remains constant over time rather than evolving based on the amount of heat extracted from the core.While these factors impact the evolution of the total CMB heat flux and the manifestation of specific CMB heat flux patterns, we do not expect them to significantly affect the amplitude of spatial and temporal CMB heat flux variations, which is the focus of our study.Below, we discuss some additional factors that are important for understanding how our results can be applied to gain insights about outer core convection.
While our models include a pressure-and temperature-dependence of the thermal conductivity, they do not explicitly incorporate how mineral phase changes affect thermal conduction.In particular, the thermal conductivity of post-perovskite has been estimated to be 20%-60% higher than that of bridgmanite (Ohta et al. 2012;Okuda et al. 2020;Wang et al. 2023).Post-perovskite is expected to predominantly be stable in the colder parts of the lowermost mantle due to the large and positive Clapeyron slope of the bridgmanite to post-perovskite transition.The higher conductivity therefore has the potential to not only enhance the total CMB heat flux, but also its spatial variations.
Another factor affecting the use of our results in geodynamo simulations is the balance of the thermal and compositional driving forces for outer core convection.The present-day geodynamo is in part driven by the release of light elements at the inner-core boundary (Nimmo 2015;Landeau et al. 2022).These light elements are not effectively accommodated by the mantle, making a vanishing flux across the CMB the most realistic boundary condition for chemical convection (Wicht & Sanchez 2019).In other words, the heterogeneous CMB heat flux would only affect thermal and not chemical buoyancy, reducing the impact on outer core convection (and therefore the effective value of q * ) the more of the driving force is contributed by chemical convection.
Finally, it is important to discuss the reference frame of our models, specifically when discussing the equatorial heat flux or the distribution of high and low heat flux patches.We here provide the CMB heat flux pattern in the reference frame of the plate reconstruction (a palaeomagnetic reference Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 31 frame derived from Tetley 2018).On geological timescales, the geodynamo coincides with the Earth's spin axis, (i.e.van Hinsbergen et al. 2015).Not taking into account small deviations on the order of several degrees between the magnetic north pole and the spin axis derived from non-dipole field components, our models using this reference frame therefore provide CMB heat flux patterns with respect to the Earth's spin axis.Note, however, that paleolongitude is generally not well-constrained in plate reconstructions.

CONCLUSIONS
In this study, we quantify the spatial variability and absolute heat flux changes at the CMB on long (billion-year) time scales using state-of-the-art numerical methods that improve the accuracy of heat flux computations in geodynamic models.We find that there are characteristic changes in the heat flux pattern throughout the supercontinent cycle, and that subduction history drives the changes in these patterns.As long as there is a layer of intrinsically dense material at the base of the mantle, cold material accumulates below the location of subduction zones at the Earth's surface-causing a large heat flux out of the core-and pushes hot material to the areas between these cold zones, which are characterised by a CMB heat flux close to zero.This behavior occurs for all combinations of material properties we have tested, including an increased viscosity of the intrinsically dense material.
The number and morphology of these hot and intrinsically dense thermochemical structures depends on subduction history and the stage of the supercontinent cycle.Stable subduction zones such as in Earth's immediate past lead to coherent and stable piles.Conversely, at times in Earth's history when subduction location changes frequently, thermochemical piles can fork or split up.If subduction zones are located predominantly in the supercontinent hemisphere during supercontinent assembly, a degree-1 pattern can develop, with high CMB heat flux (cold material) located beneath the supercontinent and low heat flux (hot thermochemical structures) located beneath the superocean.In return, a ring of subduction zones around the supercontinent during its break-up tends to lead to a degree-2 pattern, with two hot thermochemical structures separated by this girdle of cold material.Only models without a dense basal layer feature a significantly different heat flux pattern with areas of high and low heat flux alternating on much shorter spatial scales.This is because without a dense layer, cold remnants of subducted slabs can spread out along the core-mantle boundary more easily.
The amplitude of spatial heat flux variations at a given point in time is primarily affected by the material properties of the mantle.A thermal conductivity that increases with depth, a viscosity reduction in the lowermost mantle as expected for the post-perovskite phase, and the presence of an intrinsically dense layer all increase the amplitude of spatial heat flux variations.Together, these factors can increase this amplitude by a factor of 3.For a given set of material properties, the amplitude of heat flux variations (as characterized by q * ) only varies moderately (by 30-50%) over time.Our minimum estimate is q * ≥ 2 (and likely q * ≥ 6), but q * might be > 30 depending on the adiabatic heat flux out of the core.Lowermost mantle material properties, in particular the viscosity, strongly affect the timing of maxima and minima in the total heat flux, equatorial heat flux and the amplitude of spatial heat flux heterogeneity.Therefore, the current uncertainties both in plate motion history and lowermost mantle properties prohibit better constraints on the temporal evolution of the CMB heat flux and connecting them to paleomagnetic observations.However, we hope that our results will serve as a tool in future studies to better quantify the effect of CMB heat flux heterogeneity on the geodynamo and improve our understanding of the connection between past mantle flow and changes in the magnetic field behavior.
T (x, 0) = g(x) and boundary conditions: Here, ∂Ω 1 is the part of the boundary where Dirichlet boundary conditions are prescribed (prescribed temperature), ∂Ω 2 is the part of the boundary where Neumann boundary conditions are prescribed (prescribed heat flux), q 2 is any prescribed boundary heat flux on ∂Ω 2 , and n is the unit vector normal to the boundary.T is temperature, k thermal conductivity, ρ density, C p specific heat capacity, and g(x) is the initial temperature in dependence of the location x.Note that Gresho et al. (1987) and ASPECT define the normal of the boundary n in outwards direction, so that in our case heat flux into the mantle and out of the core results in negative values for the heat flux.However, we have flipped the sign of the output in all figures and in the main text for ease of understanding.
Reformulating equation ( 30) of Gresho et al. (1987) into the weak form of the consistent boundary flux method for our problem yields: In this equation we use the notation of (3), Γ i are those discretized temperature basis functions that are non-zero on the boundary, u velocity, S all heat sources (all right-hand side terms of (3)), and q 1 the heat flux through boundaries with Dirichlet conditions (which we want to solve for).Compared to Gresho et al. (1987) we additionally simplify the equation to assume we solve the Stokes equation without spurious numerical velocity divergence, and only considering Dirichlet or Neumann boundary conditions.In order to solve this equation we utilize the same type of finite elements that are used for the temperature equation (Q 2 ) on the boundary and expand Making use of a Gauss-Lobatto-Legendre quadrature for the quadrature points on the boundary (which are colocated with the finite element support points of our chosen element) ensures that the resulting linear system has only entries on the main diagonal and is therefore easily inverted.
We measure the relative error ϵ of our CBF implementation as the difference between the computed Nusselt number Nu and a reference value Nu ref computed using a Richardson extrapolation of Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 43 the model results of increasing resolution.The Richardson extrapolation value is assumed to be close to the exact value, which is supported by the observation that it is always very close (with a relative difference of 10 −3 − 10 −4 with only a few exceptions around 10 −2 ) to the values reported in Blankenbach et al. (1989) and King et al. (2010), see Table C1 and Table C2.Comparing the relative error of CBF (ϵ CBF ) to the relative error of a heat flux computed as the thermal conductivity times the temperature gradient at the model boundary (ϵ grad ) we notice that CBF results in a significantly smaller error and a better convergence rate with increasing resolution.In particular, at very coarse resolutions of 8 or 16 cells the CBF error is significantly smaller than the gradient error for well resolved simulations (Blankenbach case 1a), while for underresolved simulations both errors are of comparable magnitude.
For fine resolutions (64 cells) the CBF error is uniformly smaller by 1-4 orders of magnitude.We observe that the convergence rate of the CBF method (right panels of Figure C1) is better by at least one order (in low Ra number cases or when the flow is sufficiently resolved) and up to three orders (in high Ra number cases).In our opinion this variability in the convergence rate of the CBF method is caused by the larger number of input properties that are used to compute the heat flux.While the gradient based method only depends on the derivative of the temperature (and therefore converges with one order less than the temperature solution) the CBF heat flux depends on temperature, velocity, and energy source and sink terms and depending on which of these terms dominate, the remaining error will converge with the convergence rate of that property.C1 and C2.
Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 45 Table C1.Nusselt number and root mean square velocity of benchmark cases in Blankenbach et al. (1989)

Figure 1 .Figure 2 .
Figure 1.Density in dependence of temperature and pressure for harzburgite (left), MORB (center), and the density difference between the two compositions (right).
2014), 3430±130 K(Kim et al. 2020), or ≈3950 K with an uncertainty of 200-300 K(Pierru et al. 2022); and recent data on the melting curve of iron yield an estimate of 3760±290 K(Sinmyo et al. 2019).Our value of 3700 K leads to a jump of approximately 1200 K across the thermal boundary layer.Uncertainties in the exact value of the CMB temperature should predominanly affect the average CMB heat flux and should only have a minor effect on the pattern of heat flux variations.

Figure 3 .
Figure 3. Visual representation of the mantle convection model, with slabs (regions at least 200 K below the adiabatic temperature, blue-to-white colors indicating pressure, only shown below 670 km depth) and thermochemical piles (regions with a basalt fraction > 0.5, yellow-to-red colors indicating pressure) interacting.The African hemisphere is shown on the left and the Pacific hemisphere is shown on the right.The outer model boundary is indicated by the gray sphere, the inner boundary (where not covered by slabs or thermochemical piles) is represented by the dark yellow sphere, with the computational mesh shown in the bottom right quadrant.

Figure 4 .
Figure 4. Nusselt number (boundary heat flux) over the number of cells (resolution) for the benchmark models of Blankenbach et al. (1989) (left panel) and King et al. (2010) (right panel).We plot results of the CBF heat flux method (solid lines) and a gradient-based heat flux method (dashed lines) compared to the reference results

Figure 5 .
Figure 5. Change in average mantle temperature (a) and total heat flux out of the core (b) over time in all models.Estimates for the Earth's minimum, maximum and preferred cooling rate (from Jaupart et al. 2015) are marked as dashed gray lines for reference.Estimates for the Earth's CMB heat flux (based on Jaupart et al. 2015; Nimmo 2015) are marked with a gray background.

3. 2
Fit to present-day lowermost mantle structure Numerous studies have investigated the fit of lowermost mantle structure predicted from geodynamic models incorporating Earth's plate motion history to the LLSVPs observed in seismic tomography(McNamara & Zhong 2005;Bull et al. 2009;Zhang et al. 2010;Davies et al. 2012;Bower et al. 2013;Flament et al. 2017; Cao et al. 2021a,b;Flament et al. 2022;MacLeod et al. 2023), (see McNamara 2019, for a review).All of these studies show that cold downwellings can push aside hot material residing at the base of the mantle-both in the form of dense thermochemical piles and in form of plume clusters-shaping it into a geometry roughly matching that of the LLSVPs.The fit between regions of high temperature in such a geodynamic model to low velocities in seismic tomography models can therefore serve as model validation.

3. 3 Figure 6 .Figure 7 .
Figure 6.Temperature distribution (shown as the difference compared to the initial mantle adiabat) in a slice in 2600 km depth through all models at the end of the simulation (present-day).Bottom right panel shows the Swave velocity anomaly of the seismic tomography model SEMUCB (French & Romanowicz 2014) in the same depth slice for comparison.Yellow colors in the other panels indicate seismically slow regions (dv/v < 0.1%) in SEMUCB.

Figure 8 .
Figure 8. Changes in minimum, maximum, average (black line) and equatorial (darker colored thick line) heat flux at the core-mantle boundary over time for all models.The equatorial heat flux is defined as the average within ±20 • latitude of the equator (as in Zhang et al. 2010).All panels use the same scale, except for the bottom right panel, which is re-scaled to show the full heat flux variations in the Weak post-perovskite model.
Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 21 flux both within piles and in regions where slab remnants accumulate compared to the Thermochemical model (which has a constant conductivity).The p-T-dependent, Weak post-perovskite and Strong basalt models therefore exhibit both a larger maximum and average heat flux.At the same time, the minimum heat flux remains near zero and does not undergo significant changes because the thermal gradient in the hottest regions remains low.The p-T-dependent conductivity also impacts the shape and distribution of hot piles and cold subducted material.The lower conductivity within the litho-

Figure 9 Figure 9 .
Figure 9 also shows periods of particularly high or low reversal frequency, since geodynamo studies suggest a connection to the amplitude of the CMB heat flux hererogeneity (i.e.,Glatzmaier et al. 1999; shows the changes in these CMB heat flux patterns throughout the supercontinent cycle.All models show low CMB heat flux (yellowish colors) where hot material is located: Localized at the base of plumes in the Thermal model; more spread out within the thermochemical piles in all models with chemical heterogeneities.Since these piles convect internally (highlighted by the honeycomb-like pattern within these regions) but are too dense to advect this heat further upwards, they maintain relatively high temperatures and insulate the core.Conversely, regions where subducted slabs reach the CMB and increase the local thermal gradient feature high heat flux (bluish colors).Variations between areas of high and low heat flux occur on several different characteristic spatial scales: In the Thermal model, patches of high and low heat flux alternate on intermediate length scales (∼1000-2000 km), with cold subducted material and hot upwellings distributed relatively evenly along the whole CMB.This is in contrast to the models with dense thermochemical structures, which occupy large areas between the subducted slabs, preventing them from spreading out across the whole CMB and leading to larger-scale variations in heat flux, with large high-and low-heat-flux patches located in specific parts of the globe.In addition to these global-scale heterogeneities, there is an additional, superimposed small-scale pattern within the piles reflecting their internal convection.The characteristic length scale of these convection cells is controlled by the viscosity within the piles, with a low viscosity (as in the Weak post-perovskite model) leading to smaller-scale structures, and higher viscosity (such as in the Strong basalt model) leading to larger-scale structures.In all models, the heat flux patterns show characteristic changes throughout the supercontinent cycle, controlled by the subduction history in the plate reconstruction (Figure10, right column).Wherever subduction zones are located at the surface, cold material sinks downwards, reaching the CMB about 150-200 Myrs later and leading to an increased CMB heat flux.This becomes apparent

Figure 10 .
Figure 10.Evolution of the heat flux at the core-mantle boundary.Panels further up are further back in time.Large negative values (blue colors) indicate a large heat flux into the mantle and represent cold regions, values close to zero (yellow colors) indicate low heat flux into the mantle and represent hot regions, i.e. plumes or the hot thermochemical structures.The right column shows the location of plates (gray lines) and continents (green) in the plate reconstruction, with subduction zones marked in black.The complete evolution of all models is also displayed in Videos S1 to S5.
Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 25 flux patterns, especially in the Thermochemical, p-T-dependent and Strong basalt models at 400 Ma, which all feature one hemisphere with several high heat flux patches, while the other hemisphere consistently exhibits much lower heat flux.In particular the subduction between Australia-Antarctica and the Proto-Pacific Ocean at 600 Ma is clearly visible as a three-forked high heat flux patch between ∼60-120 • W (approximately across the area where North and South America are located today) at 400 Ma.Instead of a hemispherical distribution, at earlier times around 600 Ma, patches of high heat flux are spread out across the whole globe.This pattern matches the trench locations at 800-700 Ma, which are similarly distributed, with additionally several subduction zones forming, disappearing or changing location in the transition from Rodinia break-up to Gondwana assembly.Finally, the girdle of high heat flux surrounding the single thermochemical pile in the 800 Ma panel in the models with chemical heterogeneities corresponds to the ring of subduction zones surrounding Rodinia ∼150-50 Myr prior.How long exactly it takes for material subducted at the surface to affect the CMB heat flux, and for how long the resulting cold anomalies are preserved at the CMB depends on the material properties of the individual model.For example, in the Weak post-perovskite model the cold remnants of subducted slabs are deformed more easily, heat up faster, and impact the CMB heat flux for a shorter time, whereas convection at the CMB is more sluggish in the Strong basalt model so that high-CMB heat flux patches are preserved for a longer time.

Figure 11 .
Figure 11.Changing power of the first 4 spherical harmonics degrees of the core-mantle boundary heat flux over time.The scale at the top indicates the supercontinent cycle (times taken from Merdith et al. 2019).

Figure 12 .
Figure12.Evolution of q * for all models when subtracting the the heat flux conducted along the outer core adiabat, assumed to be 5 TW (left), 10 TW (middle), 15 TW (right).Since 15 TW is larger than the average heat flux in the Thermochemical model, this model would not have a net heat flux out of the core and is not shown in the right panel.

Figure C1 .
Figure C1.Nusselt number (boundary heat flux) over the number of cells (resolution) for the benchmark models ofBlankenbach et al. (1989) (top panels) andKing et al. (2010) (bottom panels).We plot results of the CBF heat flux method (solid lines) and a gradient based heat flux method (dashed lines) compared to the reference results reported in the publications (gray lines).Note that the bottom row only shows the case of Ra = 10 5 and the anelastic liquid approximation for different dissipation numbers to improve the visibility of the figure and that we used the case "RefVT" as reference value.The right column shows the respective convergence rate of the relative error ϵ N u over number of cells.The numeric values of the results are given in TablesC1 and C2.
(B89) using the consistent boundary flux method (Nu CBF ) and a temperature gradient based computation (Nu grad ).Note that both methods use the same solution to compute the postprocessed heat flux.The ASPECT reference result (AS) was computed using a Richardson extrapolation and the relative difference between each resolution and the reference value ϵ is computed as ϵ = |(Nu − Nu AS )/Nu AS |.The values are plotted in Figure C1 (top row).

Table 1 .
Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle 7 Model parameters.

Table 2 .
List of model variations.

Table C2 .
Nusselt number and other benchmarks properties of King et al. (2010) (models CU, UM, VT) using the consistent boundary flux method (Nu CBF ) and a temperature gradient based computation (Nu grad ).Note that both methods use the same solution to compute the postprocessed heat flux.The ASPECT reference result (AS) was computed using a Richardson extrapolation and the relative difference between each resolution and the reference value ϵ is computed as ϵ = |(Nu − Nu AS )/Nu AS |.The values are plotted in Figure C1 (bottom row).