Spatiotemporal Variations in Surface Heat Loss Imply a Heterogeneous Mantle Cooling History

Earth's heat budget is strongly influenced by spatial and temporal variations in surface heat flow caused by plate tectonic cycles. Here, we use a novel set of paleo‐seafloor age grids extending back to the mid‐Paleozoic to infer spatiotemporal variations in surface heat loss. The time‐averaged oceanic heat flow is 36.6 TW, or ∼25% greater than at present‐day. Our thermal budget for the mantle indicates that 149 K/Gyr of cooling occurred over this period, consistent with geochemical estimates of mantle cooling for the past 1 Gyr. Our analysis also suggests sustained rapid cooling of the Pacific mantle hemisphere, which may have cooled ∼50 K more than its African counterpart since 400 Ma. The extra heat released from the Pacific mantle may have been trapped there by the earlier long‐lived supercontinent Rodinia (∼1.1–0.7 Ga), and the Pacific mantle may still be hotter than the African mantle today.

and the construction and analysis of paleo-seafloor age grids (Karlsen et al., 2020;Müller et al., 2016). Thus, in order to derive a representative value of the long-term oceanic heat flow, it may be necessary to compute a time-averaged value across a wider swath of time, approximating a full supercontinent cycle. Using the first set of oceanic lithosphere age grids that extend backward to the mid-Paleozoic (400 Ma, Karlsen et al., 2020), we here reevaluate the time-dependence of oceanic heat flow and consider its implications for mantle cooling history.
We also reconstruct spatial variations in mantle heat loss across the interval spanning Pangea's assembly, duration and breakup. Such variations should cause some parts of the mantle to cool faster than others, producing spatial variability in mantle temperature. There is some geochemical and geophysical evidence for such variability. Major element geochemistry analyses of lavas from the Pacific and Atlantic ocean basins suggest that the central Atlantic upper mantle was hotter than the Pacific mantle during the 170-100 Ma KARLSEN ET AL.  Meridians (60°W and 120°E) drawn in magenta indicate the separation of the Pacific and African hemispheres. (c) Seafloor age-area distributions based on the paleo-seafloor age girds generated after Karlsen et al. (2020). interval and may have cooled by up to 150 K until the present (Brandl et al., 2013). This cooling released heat that was likely trapped by continental insulation of the central Atlantic by Pangea (∼320-180 Ma). Increases in mantle potential temperature (T p ) of similar magnitude, caused by supercontinent insulation, have also been reported in numerical and laboratory experiments (Lenardic et al., 2011). For the present-day, Brandl et al. (2013) present thermometric evidence (temperature of average zero-age mid-ocean ridge basalts (MORBs)) that the Pacific mantle domain is 25-30 K hotter than the African domain (the Atlantic and Indian ridges). This finding agrees with Dalton et al.'s. (2014) combined analysis of geophysical and geochemical data, which indicates a higher mantle potential temperature in the Pacific. Such hemispheric variability in mantle temperature should be linked to spatial variations in the mantle's cooling rate. Here, we use models for Earth's surface heat flow during the past 400 Myr to estimate individual cooling histories for the Pacific and African mantle domains, and briefly entertain some simple scenarios that could explain their potentially distinct thermal evolution.

Estimating Mantle Heat Loss Variations
For the present-day, Earth's total surface heat loss is estimated to be 46 TW. Of this, 29 TW is attributed to the cooling of oceanic crust (hereafter "oceanic heat loss"), 3 TW is delivered to the surface via plumes and 14 TW represents heat loss from the continents (Jaupart et al., 2015). The continental component of 14 TW includes 8 TW of radiogenic heat produced within the continental crust and lithosphere. Here, we estimate changes to the heat flow conducted from the mantle through the oceanic and continental lithosphere (currently 35 TW; the 3 TW plume heat flux is treated separately according to the description in the extended methods section, supporting information Text S1). We focus particularly on the oceanic heat loss, both because it is the largest component, it exhibits large variations in space and time, and it can be estimated directly from tectonic reconstructions of past seafloor.
We compute oceanic heat loss (Q ocean ) from the mantle by applying the age-heat flow relation (see Text S1) of Hasterok (2013) to paleo-seafloor age grids that extend back to 400 Ma, generated using the method of Karlsen et al. (2020). The plate reconstruction that was used to generate these age grids was the updated version of Matthews et al. (2016), which includes corrections for the Pacific after Torsvik et al. (2019). Mantle heat loss through the continental lithosphere (Q cont ) is taken to be 6 TW, evenly spread over the continental area. Note that this excludes the 8 TW of radiogenic heat production in the crust and lithospheric mantle (Jaupart et al., 2015). Although the mantle heat loss through the continents (∼6 TW) might vary with time, the amplitudes of its variations are likely small compared to the more than four times larger oceanic heat flow (∼29 TW at present-day). Snapshots of the resulting global heat loss grids covering the past 400 Myr are shown in Figure 1a. Note the strong time-dependence, and that present-day heat loss from Earth's mantle represents an absolute minima for the 400-0 Ma time period, having been up to ∼15 TW higher in the past ( Figure 1b). These variations are driven by changes to the age-area distribution of the seafloor (Figure 1c), which are caused by fluctuations in the rates of seafloor spreading and subduction during the supercontinent cycle (e.g., East et al., 2020;Karlsen et al., 2019a). This emphasizes the need for a long time-series (ideally spanning one complete supercontinent cycle) of heat flow to robustly infer long-term average values. Based on the 400 Myr reconstruction of seafloor ages, we estimate the time-averaged value for oceanic heat loss to be 36.6 TW (Figure 1b), which is ∼25% higher than the present-day value of ∼29 TW (Figure 1b). The present-day low in oceanic heat flow has also been noted before, although based on significantly shorter time-histories (Crameri et al., 2019;Loyd et al., 2007).
To investigate spatial variations in time-integrated surface heat loss, we computed the total heat loss accumulated over the past 400 Myr (Figures 2a and 2b). Our results show that even over the time scale of several hundred million years, heat loss from the mantle is far from uniform. For instance, the region of the Pacific overlying the large low-shear-velocity province (LLSVP) "Jason" has lost 2-3 times as much heat as the region above "Tuzo", the African LLSVP ( Figure 2b). This is partly due to the time-dependent distribution of continental masses and the assembly of Pangea over Tuzo, which insulated the mantle beneath Pangea but left the Pacific "exposed" to more rapid cooling via cycling of oceanic lithosphere, and partly due to the persistence of fast-spreading ridge systems (now the East Pacific Rise) positioned over Jason. Note, however, that the position of ridges in the Pacific-Panthalassic Ocean is entirely inferred before ∼150 Ma (Torsvik et al., 2019); we will return to this limitation later.

Implications for Mantle Cooling History
To estimate the impact of spatiotemporal variations in surface heat loss on the mantle's thermal evolution, we formulated a heat budget for the mantle, which we further break down into individual heat budgets for the Pacific and African mantle domains (here approximated as two hemispheres, Figures 1a and 2, delineated as described in supporting information S1). These two mantle domains are considered separate and persistent degree-two convection cells , each bounded, at least since Pangea formation, by the long-lived circum-polar belt of subduction and persistent downwelling (e.g., Torsvik et al., 2016). The fact that slabs tend to continually sink between the LLSVPs ( Figure S1) indicates little mixing, and thus limited heat transport, between the Pacific and African mantle domains. Assuming negligible heat exchange between these domains (thermal conduction across such large spatial scales is small, Figure S2), we can treat the thermal history of each domain separately. Indeed, limited thermal mixing by convection, combined with an uneven distribution of continental landmasses, has been shown to cause large lateral temperature variations in numerical and laboratory experiments (Lenardic et al., 2011).
We construct a mantle heat budget following Jaupart et al. (2015), taking mantle heat sources to be 11 TW (Q radio ) from radiogenic heat production (which excludes the 8 TW from the continental crust and lithosphere), plus 11 TW (Q core ) from the core. For heat sinks we include our estimates of oceanic KARLSEN ET AL.  (Q ocean ) + continental (Q cont ) heat loss (Figure 1b, black line) and an additional 3 TW (Q hotspots ) heat output from hotspots (Arnould et al., 2020;Jaupart et al., 2015). This value might be relatively conservative, considering a recent estimate based on excess topography (including a magmatic swells) suggesting that the present-day plume heat flux may be as high as ∼6.3 TW (Hoggard et al., 2020). The heat sources and plume heat flux are assumed to be constant through time and evenly distributed between the Pacific and African mantle domains. This assumption allows us to isolate the effect of variations in oceanic heat loss on mantle cooling. A relatively stable plume heat flux with time, around 3 TW, is supported by state of the art mantle convection models with Earth-like tectonics on the surface (Arnould et al., 2020), while time-variations in the core heat flux were probably less than ∼3 TW during the Phanerozoic (T. Nakagawa & Tackley, 2011;Olson, 2016;N. Zhang & Zhong, 2011). By comparison, the continuous surface heat loss estimated herein has varied by ∼15 TW since 400 Ma (Figure 1b, black line).
Based on the mantle heat budget adopted, we can calculate rates of mantle cooling for the entire mantle and for each of the two domains separately using: where, c = 1250 J/(K kg) and m = 4.01•10 24 kg (or m/2 for each hemisphere; see supporting information Text S1 for more details). The resulting estimates indicate that since 400 Ma the whole mantle has cooled ∼60 K (Figure 3c), or at a time-averaged rate of 149 K/Gyr, with a maximum of ∼200 K/Gyr at 350 Ma, and a minimum of ∼100 K/Gyr at the present-day (Figure 3b). However, that cooling was unevenly distributed with significantly more heat having left through the Pacific hemisphere. Our calculations suggest that the Pacific cooled at a time-averaged rate of 208 K/Gyr since 400 Ma, more than twice that of the 90 K/Gyr time-averaged rate for the African hemisphere (Figure 3b). In the context of cumulative heat loss, the Pacific and African domains cooled by 83 K and 36 K, respectively since 400 Ma, implying a differential in mantle cooling of almost 50 K since at least the mid-Paleozoic (Figure 3c).

Implications for the Whole-Mantle: Heat Budget and Cooling History
The mantle heat budget and its related cooling history are often inferred based on present-day values of surface heat loss (e.g., Jaupart et al., 2015). Here we have shown, using a novel set of seafloor age grids that extend back to the mid-Paleozoic (400 Ma), that the present-day likely represents a time of anomalously low surface heat flow caused by an "all time old" configuration of ocean basins (see e.g., Müller et al., 2016 for a review). Considering that the time-averaged value of mantle heat loss from cooling oceanic lithosphere might be closer to 36.6 TW, rather than ∼29 TW (Figure 1b), the long-term rate of planetary cooling may be significantly higher than previous estimates based on the present state of the Earth.
Using source temperatures from non-arc basaltic rocks covering the past 4 Gyr (Herzberg et al., 2010), and assuming a present-day mantle potential temperature of a little more than 1,400 °C, Jaupart et al. (2015) estimated the mantle's long-term cooling rate to be 50-100 K/Gyr. However, if we instead consider the data from Herzberg et al. (2010) only for the past 1 Gyr and use a present-day potential temperature of 1,350 °C, the least squares solution suggests a more recent mantle cooling rate of 155 K/Gyr ( Figure S3). This rate is comparable to cooling rate estimates for the past 1 Gyr from Korenaga (2008aKorenaga ( , 2008b, and is close to our time-averaged whole-mantle cooling rate of 149 K/Gyr for the past 400 Myr (Figure 3b). An alternative explanation, which would reconcile our high time-averaged oceanic heat flow (Figure 1b) with the slower cooling rates of Jaupart et al. (2015) (50-100 K/Gyr), is that the poorly constrained internal heat sources are larger than the typical "preferred values" that we used in our mantle heat budget (Section 3). Indeed, we assumed a core heat flux of 11 TW and a radiogenic mantle heat production of 11 TW, but Jaupart et al. (2015) suggest that these could be as high as 17 TW and 17 TW, respectively. We therefore explore the effect of higher (and lower) values for these heat sources on time-averaged mantle cooling rates and total mantle cooling over the past 400 Myr in Figure S4. Note that the nearly 50 K of accumulated temperature difference between the Pacific and African mantle domains since 400 Ma is independent of the magnitude of these heat sources, as long as they are evenly distributed between the domains (Figures S4e and S4f).

Is the Pacific Mantle Hotter than the African Mantle?
Although our calculations suggest that the Pacific hemisphere has been cooling at a higher rate, this does not necessarily mean that it is currently colder than the African hemisphere. The currently higher heat flow in the Pacific may be the manifestation of a hotter mantle beneath the Pacific (Brandl et al., 2013;Dalton et al., 2014). Given our rates of cooling (Figure 3b), this would imply that the Pacific has been hotter for more than 400 Myrs, and remains so, despite its faster cooling. The initial heating may have been caused by supercontinental thermal insulation (e.g., Lenardic et al., 2011) by the Precambrian supercontinent Rodinia, which may have covered the paleo-Pacific (Panthalassic) mantle domain (e.g., Tegner et al., 2019). In the case of an "ideal" supercontinent configuration (i.e., at ∼250 Ma, Figure 1a), one hemisphere is stripped of continents and cools about 150 K/Gyr faster than the one occupied by the supercontinent (Figure 3b). This suggests that if Rodinia lasted for about 400 Myr (∼1.1-0.7 Ga, e.g., Torsvik & Cocks, 2017), the Pacific mantle domain may have accumulated heat corresponding to an extra 60 K compared to the African. Given that the rate of mantle heat loss in the Pacific hemisphere may not have increased above its African counterpart until the Pacific lost most of its continents, the heat accumulation within the Pacific domain may have persisted up until about ∼370-350 Ma (Figures 1a and 3a). In this scenario, much of that "residual" (Rodinian) heat has been lost in the past 400 Myr, but some still remains (as the Pacific mantle domain has cooled ∼50 K relative to the African over this period). Of course, this assumes that heat transfer between these hemispheres has been limited. A "long-lived hotter Pacific"-scenario could explain the consistently KARLSEN ET AL.
10.1029/2020GL092119 6 of 10 higher plate velocities of the Pacific hemisphere during the past 400 Myr (Figure S5), because a hotter mantle implies lower mantle viscosity and smaller resistance to plate motions.

Constraints on Mantle Temperature Variations
The notion of a persistently hotter Pacific mantle domain may seem partly contradictory to geochemical and crustal thickness constraints. The compositions of MORBs from the Atlantic and Pacific basins suggest that the mantle potential temperature (T p ) beneath the early central Atlantic (∼170 Ma) was up to 150 K hotter than at present (Brandl et al., 2013). Finding no evidence of comparably elevated mantle temperatures for older Pacific MORB within the 170-5 Ma age range, Brandl et al. (2013) attributed the central Atlantic thermal anomaly to continental insulation. By contrast, the average T p recorded by zero-age Pacific and central Atlantic basalts are about 1,413 °C and 1,385 °C, respectively (Brandl et al., 2013), suggesting a currently hotter Pacific domain. A slightly hotter present-day Pacific mantle is in agreement with our Rodinia scenario (Section 4.2). A regression of oceanic crustal thickness for the Pacific and African hemispheres over the past 170 Myr presented by Van Avendonk et al. (2017), also supports a currently hotter Pacific but a faster-cooling African domain. Their results indicate that the present average crustal thickness, broadly corresponding to T p , of the East Pacific Rise exceeds that of the Atlantic and Indian Ocean ridges, but that the crustal thickness in the African hemisphere has thinned considerably more than in the Pacific hemisphere during the past 170 Myr. The strong decrease in oceanic crustal thickness in the African domain corresponds to an estimated cooling of 390 K/Gyr. Such rapid cooling is also reflected by the T p decrease of 880 K/Gyr estimated for the central Atlantic by Brandl et al. (2013). Both estimates are for the last 170 Myr. The major caveat to these geothermometric and crustal thickness constraints, however, is that they may reflect more fertile source compositions combined with a hot asthenosphere shortly after continental rifting (supporting information Text S2).
The apparent faster cooling of the African hemisphere suggested by van Avendonk et al. (2017) and Brandl et al. (2013) may be limited to just the upper mantle. In this case, the geochemically constrained rapid cooling of the African domain could represent the release of heat trapped shallowly beneath Pangea. Slow thermal mixing within the African domain, as indicated by the slow plate speeds there ( Figure S5), would allow rapid cooling of the uppermost mantle without rapid cooling of the entire African domain (upper and lower mantles). In comparison, the Pacific mantle domain, with its fast spreading ridges and plates ( Figure S5), may be more vigorously convecting and well-mixed (King & Adam, 2014). Therefore, it is not straightforward to compare our hemispheric cooling rates (Figure 3b), which represent averages of the full mantle depth, with cooling rates inferred from analyses of oceanic crust. The thick oceanic crust and elevated T p estimates of the early Atlantic and Indian ocean basin might also be related to more fertile upper mantle sources resulting from the Pangea assembly and breakup processes, as well as from plume head materials supplied by a number of large igneous provinces (e.g., Doucet et al., 2020;Torsvik et al., 2020;supporting information Text S2). Alternatively, additional heat sources in the Pacific mantle, such as elevated core heat flow or extra heat-producing elements, could explain the more rapid long-term heat flow there, but there is currently no evidence of such a source (Olson, 2016).

Limitations
The quality of our surface heat loss calculations depends on paleo-seafloor age grids generated based on the plate tectonic reconstructions of Matthews et al. (2016), using the algorithm of Karlsen et al. (2020). This algorithm has been shown to introduce only negligible errors beyond those already present in the underlying plate tectonic model. Thus, the quality of the age-grids depends on the restoration of mid-ocean ridges and subduction zones, and the reconstruction of global plate kinematics, for which formal errors are difficult to impossible to estimate. Despite being unquantified, it is reasonable to assume that these spatial uncertainties are very significant for time-periods lacking preserved oceanic lithosphere (i.e., older than ∼200 Ma). Nevertheless, our main conclusions, namely (1) that the present-day represents a time of anomalously low surface heat flow and (2) that the Pacific mantle domain has been cooling more rapidly than its African counterpart, persist even if we consider only the last 200 Myr.
Notwithstanding that the spatial uncertainties may be very large for times lacking preserved oceanic crust, estimates of global surface heat flow for these times may still be reasonable. This is because global sea level variations on time-scales of ∼10 8 yr are dominated by ocean basin volume changes driven by the time-dependent seafloor age-area distribution (e.g., Conrad, 2013;Karlsen et al., 2019a;Müller et al., 2008). Thus, the sea level record serves as a first-order proxy for oceanic surface heat flow (Harrison, 1980; Figure S6). Karlsen et al. (2020) demonstrated that their seafloor age grids show good agreement with the history of sea level extending back to 400 Ma. With this in mind, we suggest that the temporal variations in global heat flow presented here are likely reasonable, but the spatial variations in oceanic heat flow for times older than ∼200 Ma should be cautiously interpreted.

Conclusions
Earth's thermal evolution is the result of an imbalance between surface heat loss and internal heat production. Here, we have shown that the former is strongly time-dependent, and far from spatially uniform. The assembly and breakup of supercontinents alters the age-area distribution of the ocean basins, resulting in uneven thermal insulation and cooling of the mantle. These variations create large lateral temperature variations between mantle domains. We summarize our main findings as follows: -Surface heat flow, just like the sea level (Conrad, 2013; Figure S6), shows peaks during the assembly and dispersal phases of the supercontinent cycle, and lows during the lifetime of the supercontinent (Pangea) and periods of maximally dispersed continents (such as present-day). Over the past 400 Myr, variations in oceanic heat flow ranged from a maximum of ∼44 TW during the assembly phase of Pangea (350 Ma), to an all-time low of ∼29 TW at present-day, with a time-averaged value over the supercontinent cycle of 36.6 TW (Figure 1b). This loss of heat exceeds that gained from mantle heat sources, and we estimate overall cooling of 149 K/Gyr since 400 Ma, consistent with observations of mantle potential temperature during the past 1 Gyr (Herzberg et al., 2010; Figure S3). -The proportion of fast-spreading ridges and area covered by cooling oceanic crust determines the relative heat loss between the Pacific and African hemispheres. We find that high-spreading activity in the Pacific hemisphere ( Figure 1a) led to greater heat flow there and drove ∼50 K more cooling in the Pacific mantle domain compared to the African domain during the past 400 Myr. This excess heat may have accumulated during an earlier period of slow cooling in the Pacific domain, caused by the presence of the longlived supercontinent Rodinia (∼1.1-0.7 Ga). The two hemispheres may thus have experienced alternating periods of slower and faster cooling, according to the formation and dispersal of supercontinents.
Looking ahead, we consider improved plate tectonic reconstructions, combined with better geophysical and geochemical constraints on large-scale temperature variations within the mantle, both past and present, to be crucial for a greater understanding of hemispheric differences in mantle temperature. We contend that such temperature differences result from tectonic cycles at the surface, and we expect them to have significantly impacted Earth's thermal evolution.

Data Availability Statement
All the data and software necessary to reproduce our results, obtained as described in Section 2, are publicly available. The plate model used can be downloaded from https://www.earthbyte.org/global-plate-boundary-evolution-and-kinematics-since-the-late-paleozoic/and the code used to calculate the seafloor age grids (Karlsen et al., 2019b)