On the Relation Between Basal Erosion of the Lithosphere and Surface Heat Flux for Continental Plume Tracks

While hotspot tracks beneath thin oceanic lithosphere are visible as volcanic island chains, the plume‐lithosphere interaction for thick continental or cratonic lithosphere often remains hidden due to the lack of volcanism. To identify plume tracks with missing volcanism, we characterize the amplitude and timing of surface heat flux anomalies following a plume‐lithosphere interaction using mantle convection models. Our numerical results confirm an analytical relationship in which surface heat flux increases with the extent of lithosphere thinning, which is primarily controlled by the viscosity structure of the lower lithosphere and the asthenosphere. We find that lithosphere thinning is greatest when the plate is above the plume conduit, while the maximum heat flux anomaly occurs about 40–140 Myr later. Therefore, younger continental and cratonic plume tracks can be identified by observed lithosphere thinning, and older tracks by an increased surface heat flux, even if they lack extrusive magmatism.

Plain Language Summary Extra heat is transmitted through Earth's tectonic plates above hot upwellings called plumes, and contributes to the heat budget of glaciated regions such as Greenland. A thin oceanic plate moving over such an upwelling typically yields a chain of age-progressing volcanic islands such as Hawaii. However, such volcanism is often missing in continental regions such as Africa or Greenland, which have thicker plates. Nonetheless, the passage of the plate over the upwelling leaves a trace in the form of a reduced plate thickness and an increased amount of emitted heat even millions of years after the plume passage. In this study, we use numerical models of mantle convection to show that these two observations are directly linked to each other, with a larger heat flux anomaly observed for areas that have been thinned more extensively. Furthermore, we demonstrate that there is a significant time delay between the thinning, which happens during the passage of the upwelling, and the heat flux that is observed many millions of years later. As a consequence, past interactions of plates with plumes can influence today's heat output in continental regions.
• We used numerical and analytical approaches to characterize heat flux anomalies following plume impingement on the lithosphere • Surface heat flux anomalies follow an analytical relationship that predicts increasing heat flux with increasing basal lithospheric erosion • Lithosphere thinning is mostly controlled by viscosity structure, with a maximum surface heat flux anomaly following 40-140 Myr later

Supporting Information:
Supporting Information may be found in the online version of this article.
10.1029/2022GL098003 2 of 10 track in the eastern United States (Chu et al., 2013;Yang & Leng, 2014). All of these studies agree that plumes erode the lithosphere, but it is still unclear whether the cratonic root can heal (as suggested for the Slave Craton [Liu et al., 2021] or the Gondwana craton in Africa and South America [Hu et al., 2018]) or if the destruction is irreversible (as seen in Africa, Celli et al., 2020). Furthermore, the tectonic setting and age differ significantly between the examples given above, making it impossible to identify systematic trends in the amount of lithospheric thinning associated with plume-lithosphere interaction.
Finally, there are two examples of plume-craton interaction for which extensive focus has been placed on heat flux: Marie Byrd Land in Antarctica (Lösing et al., 2020;Maule et al., 2005;Seroussi et al., 2017;Shen et al., 2020), and the Iceland plume in Greenland (e.g., Colgan et al., 2021;Martos et al., 2018). In both locations, the heat flux associated with the plume can have significant effects on the melting rates of ice (Rogozhina et al., 2016;Rysgaard et al., 2018;Smith-Johnsen et al., 2020). However, these two scenarios are quite different: while Antarctica has been rather stationary over the plume (Seroussi et al., 2017), Greenland passed over the plume within about 30-90 Myr (for an overview of potential hotspot tracks, see Martos et al., 2018). For Greenland, seismic studies also suggest a thinned area above the plume track (e.g., Celli et al., 2021;Mordret, 2018), although the reported anomalies for heat flux, lithospheric thinning and reconstructed plume tracks vary significantly between different studies and often do not coincide. Therefore, a discrimination between the different scenarios requires a better and more systematic understanding of the interaction between plumes and cratonic lithosphere, and the effect of this interaction on geophysical observables. Therefore, in contrast to previous studies (e.g., Yang & Leng, 2014), we develop both numerical and analytical approaches that enable us to characterize the relationship between surface heat flux and lithospheric thinning.

Model Setup
We use 2-D and 3-D numerical models of mantle convection in Cartesian geometry to systematically investigate the connection between surface heat flux and lithosphere thinning. Modeling is done with the open source finite element code ASPECT (Bangerth et al., 2020;Heister et al., 2017;Kronbichler et al., 2012), and we use boxes of 1,200 km by 600 km or 5,500 km by 800 km (x by z dimensions, 5,500 km by 2,000 km by 800 km for 3-D) for stationary and moving plate models, respectively. We define the initial temperature field as a linear gradient from the isothermal surface (set to 273 K) to the lithosphere-asthenosphere boundary represented by the 1500 K (or in some models the 1623 K) isotherm at a chosen depth, underlain by an upper mantle with a temperature gradient of about 1 K/km and a temperature perturbation added at the base to generate a plume. For cases with stationary plates, all domain sides are free-slip, except for the surface which is set to no-slip, and plume sources are removed after 50, 100, or 200 Myr. We obtain cases with moving plates by imposing a velocity component parallel to the x-axis on the top and uppermost (i.e., lithospheric) side walls, and a velocity that decreases linearly through the asthenosphere to a free-slip (no influx) condition below the base of the asthenosphere. All the remaining sides are set to free-slip. We neglect compressibilty and depth dependence of parameters (except for viscosity), as we expect their impact to be negligible for the purpose of this study (Albers & Christensen, 1996).
Viscosity is implemented via a modified diffusion-dislocation creep law with Parameters in the equation above are the layer viscosity scaling η j (implementing a viscosity jump beneath the asthenosphere), the reference viscosity η ref , the prefactor A, the grain size d with the grain size exponent m, the strain rate with the stress exponent n, the activation energy E and volume V, the gas constant R, pressure P and temperature T. The subscript i refers to either diffusion or dislocation creep. Depending on the choice of parameters, the resulting viscosity can be simplified to a Newtonian (n = 1) and grain size independent (m = 0) viscosity, as we choose in most of our simulations. In these models, we further simplify the equation by setting V = 0, making the viscosity only temperature-dependent. A change of viscosity with depth is then added as a step-wise viscosity increase below the asthenosphere via η j . This means that we can control the temperature-dependence of viscosity via the activation energy, resulting in a viscosity increase with decreasing temperature. In some models, we employ a higher initial temperature for the lithosphere-asthenosphere boundary to decrease the lower lithosphere viscosity. All input parameters for the viscosity law are defined in Table S1 in Supporting Information S1.
Our setup gives us direct control over several parameters that may affect the heat flux and lithosphere thinning: lithosphere thickness and viscosity, asthenosphere thickness and viscosity, plume excess temperature, and the interaction time between the plume and the lithosphere (via plume lifetime or plate velocity). Plume excess temperature relates to the plume flux via the surface integral over a plume cross-section: B = ∫v r ραT P dS, with density ρ, thermal expansivity α, plume excess temperature T P , the vertical velocity v r , and surface increment dS (Farnetani et al., 2018;Heyn et al., 2020). Estimated buoyancy fluxes are about 1 Mg/s for our reference cases and 1.5 Mg/s for models with T P = 400 K, which falls within the range of observed fluxes on Earth (e.g., Hoggard et al., 2020, and references therein). However, as discussed in Text S2 in Supporting Information S1, plume excess temperature is a better parameter for estimating the effect of plume parameters on expected lithosphere thinning and heat flux anomalies. In order to quantify their impact on the predicted heat flux and lithosphere thinning anomalies, we test a wide range of parameter combinations (see Tables S2-S4 in Supporting Information S1), with the ranges given in Table S1 in Supporting Information S1. Due to computational costs, only parameters that are shown to be important in 2-D have been tested in 3-D geometry.
To quantify the effect of the plume on both lithospheric thickness and surface heat flux, we measure these plume-generated anomalies for every model relative to their average values over the first 300 km of x-direction (upstream of the plume) for each time step, which ideally represents a plume-unaffected lithosphere. This choice of reference value makes our results more easily comparable to heat flux and lithospheric structure observations since we only require information from one snapshot in time, and eliminate the effects of conductive lithosphere growth. However, the reference can be affected significantly by perturbations within the reference area, such as dripping instabilities that may temporarily appear to thicken the reference lithosphere, and thus cause an apparent increase in lithosphere thinning. We therefore use the 1400 K isotherm instead of 1500 or 1623 K as a reference for our calculation of the lithosphere thinning, as this isotherm is less affected by drip formation than the actual LAB, while still exhibiting significant thinning. Finally, since the reference values for both heat flux and lithosphere thickness vary within and between models, we express these anomalies as a percent of the reference value. For stationary plate cases, we only track the evolution of the central maximum above the plume, while for moving plate cases, we track the overall maximum in time and space as well as the anomaly values at positions fixed in space (i.e., fixed x-positions) and fixed to the plate (i.e., moving with the plate over time).

Relationship Between Heat Flux and Lithospheric Thinning
As soon the rising plume reaches the lithosphere-asthenosphere boundary, it starts to spread laterally beneath the lithosphere, either symmetrically or asymmetrically depending on whether the plate is stationary (Figure 1a) or moving (Figures 1b and 1c). The temperature and viscosity anomalies related to the emplacement of hotter than ambient material beneath the lithosphere erode the base of the lithosphere both mechanically via the formation of drips (visible as cold extensions below the lithosphere, expressed as negative values of lithospheric thinning, Figure 1), and by locally uplifting isotherms via steepening the lithospheric temperature gradient (green lines in Figures 1a and 1b and isosurface in Figure 1c, lower panels). Drips form predominantly next to the plume track, where the plume pushes the eroded lithosphere, but some also form in the plume-affected area downstream of the plume. The respective lithosphere thinning and surface heat flux anomalies are shown in Figure 1 above the snapshots of the temperature field. As can be seen, the patterns of heat flux and lithosphere thinning are directly comparable to each other for stationary plate cases (Figure 1a), while the anomalies differ significantly at any given time for moving plate models (Figures 1b and 1c). The short wavelengths of lithosphere thinning are not reflected in the heat flux in either case.
In order to better understand the timing and relation between the thinning and heat flux anomalies, we track their evolution for two positions that move along with the plate. Figure 2 shows the relative lithosphere thinning (blue lines) and relative surface heat flux (red lines) for two 3-D cases with different plume excess temperatures. In both cases, lithosphere thinning starts significantly earlier than the onset of the surface heat flux, and may even reach a maximum value before a detectable heat flux anomaly has evolved. For most models, the maximum heat flux anomaly is delayed by about 40-140 Myr compared to the maximum lithosphere thinning, with (initially) thinner or more thinned lithosphere resulting in smaller delays. Apart from the timing, Figure 2 also shows that more extensive thinning, for example, due to a hotter plume (right panel), results in a larger heat flux anomaly. The same trend can be observed by comparing different positions on the same plate, where stronger thinning is followed by a larger heat flux anomaly (dashed and solid lines in Figure 2).

Parameter Dependence of Anomalies
The amplitudes of plume-induced anomalies vary significantly among different model runs (Figures 1 and 2). Although all reference cases have exactly the same parameters (except whether the plate is moving or not), predicted anomalies are largest for stationary plates due to the prolonged time of plume-lithosphere interaction. As a consequence, a shorter plume lifetime or a faster plate velocity reduces heat flux and thinning anomalies, while a longer plume lifetime or a lower plate velocity have the opposite effect. Apart from plate velocity, several other parameters affect how much the lithosphere is thinned by the plume, and therefore directly affect the predicted surface heat flux anomaly. Figure 3 shows the temporal evolution of relative surface heat flux anomalies and relative lithosphere thinning for 6 different cases of stationary plates. As can be seen, a lower asthenosphere viscosity facilitates lithosphere erosion ( Figure 3b) and results in larger surface heat flux anomalies, as can also be seen in Figure 2. Lower lithosphere viscosity produces an even stronger trend, with a weaker lithosphere being eroded more easily (compare Figure 4 and Tables S2-S4 in Supporting Information S1). A stress-dependent (non-Newtonian) viscosity can also significantly reduce the viscosity locally above the plume due to high strain-rate (Equation 2), resulting in stronger lithosphere thinning and elevated heat flux ( Figure 4). In contrast, the plume excess temperature has a relatively minor impact on the anomalies (Figures 3c and 3d), especially for moving plate cases. The same holds true for the initial lithosphere or asthenosphere thickness, although the absolute thinning is smaller/bigger for thinner/thicker plates (Tables S2-S4 in Supporting Information S1). The maximum heat flux anomaly occurs many 10s of Myr after the maximum lithospheric thinning for all parameter combinations (Figure 3), consistent with our previous findings (Figure 2).
Even though the maximum anomalies of lithosphere thinning and heat flux are not observed at the same time, they are directly related. Figure 4 shows the predicted relative maximum surface heat flux anomaly as a function of maximum relative lithosphere thinning for all tested parameters and geometries. As can be seen, all models fall along the same trend, independent of the geometry (2-D or 3-D) or whether the plate is stationary or moving, with most of the tested stationary plate cases exhibiting significantly larger values for both anomalies compared to the corresponding moving plate cases (see also Figure 1). Models with a stationary plate can reach heat flux anomalies of up to about 52% (about 16.3 mW/m 2 ) for an asthenosphere with η Asth = 1e18 Pa⋅s (with about 39% or 49.6 km of local thinning), or even up to 80% (31.3 mW/m 2 , with 50% or 80 km thinning) for a model with diffusion-dislocation creep. The highest absolute heat flux of 43.7 mW/m 2 (60% relative anomaly) is observed for a weak lower lithosphere with activation energy E = 100 kJ/mol, but this model features ample small-scale convection with unstable lithosphere and thus a rather unreliable lithosphere thinning (130.4 km or 58%). In contrast, 2-D and 3-D moving plate cases mostly cluster below heat flux anomalies of 30% (about 10 mW/m 2 ), with only very few models reaching higher values. Heat flux anomalies of 20% or more (or at least about 6 mW/m 2 ) can only be reached for asthenosphere viscosities of η Asth ≤ 1e18 Pa⋅s and/or a weak lower lithosphere, which might be convectively unstable anyway.

Analytical Solution
We can compare our models to an analytical solution for the time-dependent thermal structure of a stationary system. If we instantaneously offset the lithosphere-asthenosphere boundary by Δh (inlay in Figure 4) and let the system equilibrate thermally, the surface heat flux will develop a time-dependent anomaly Δq(t). Starting from Carslaw and Jaeger (1959) Chapters 3.3 and 3.4, we can express the time-dependent temperature profile in the thinned lithosphere as with the surface and LAB temperatures T 0 and T m , the initial and reduced lithosphere thicknesses L and l = L − Δh, the thermal diffusivity κ and the differential temperature gradient = ( − 0) ( 1 − 1 ) . Using the (1−Δℎ ) . A detailed derivation is given in the Supporting Information S1.
In the time limits of zero (t → 0, no time to equilibrate) or infinity (t → ∞, fully equilibrated), the heat flux anomaly approaches Δq rel = 0 and Δ = Δℎ 1−Δℎ , respectively. The red line in Figure 4 represents the maximum heat flux of a fully-equilibrated system (t → ∞), while the black and gray lines in Figure 4 are obtained for different values of equilibration time t assuming our reference L = 137.8 km. As can be seen, the gray lines approach the maximum heat flux for larger values of t. However, none of our dynamic models reach that maximum, because the plume dynamic changes and the lithosphere starts to regrow following maximal thinning (see Figure 2). Most of our models result in anomalies similar to a stationary model equilibrating for about 40-100 Myr, with an optimal value of t = 60 Myr.

Discussion
Our results show a clear causal relation between a plume eroding the lithosphere at the base, and a surface heat flux anomaly associated with the passage of a plume. In order to produce a significant heat flux anomaly (>30% or >10 mW/m 2 ), the lithosphere must be thinned by more than 30% (about 40 km). The trend we identify is independent of the chosen geometry or whether the plate is moving or not, and corresponds to the analytical solution for instantaneous thinning of a thermally conducting layer (Equation 3). However, due to the limited interaction time between plume and lithosphere, the system in our models (and in Earth) usually cannot fully adjust thermally to the emplaced temperature anomaly before the lithosphere starts to regrow (thicken) again. As a consequence, the analytical solution can be used to estimate expected heat flux anomalies, and to place an upper bound on the possible surface heat flux for any value of lithosphere thinning. Based on the identified trends, our results indicate that the lithosphere in dynamic models may stay close enough to maximum thinning for about 40-100 Myr (using Equation 4 and L = 137.8 km) to generate a corresponding maximum heat flux anomaly before it starts to regrow again. This time frame is also approximately reflected in the delay times of heat flux anomalies (40-140 Myr), although these delay times are more variable because they depend on the effective lithosphere thickness and the time-integrated dynamics of the system.
Of course our analysis is simplified in several ways. One important simplification is the absence of melt, which may alter the local heat flux and its delay time significantly. If melt intrudes into the lithosphere, hot material would infiltrate up to much shallower depths, rapidly increasing local heat flux without affecting lithosphere thinning (Von Herzen et al., 1989). However, for continental and especially cratonic lithosphere, intrusive volcanism appears to be limited (e.g., Chu et al., 2013;Yang & Leng, 2014), and extrusive volcanism is only observed in localized areas (e.g., Knott et al., 2020). As the absence of continuous volcanism indicates, the ability of melt to intrude into the lower cratonic lithosphere might be limited (Aulbach et al., 2017). As a consequence, we may expect a primarily conductive heat flux anomaly along most of a continental plume track (Martos et al., 2018). However, future work should include melt and melt dynamics in order to understand how melt can locally and regionally affect the surface heat flux.
We also did not include a variable radiogenic heat or any compositional differences within the lithosphere. A higher radiogenic heat production in the upper crust (Martos et al., 2018) would increase overall heat fluxes, and thus reduce relative heat flux anomalies, while lateral variations in radiogenic heat production can cause apparent anomalies. Otherwise, we do not expect the shallow lithospheric structure, for example, the Moho discontinuity between the crust and the lithospheric mantle (Mordret, 2018), to have any significant effect on the predicted anomalies, except for melt intrusions. Thus, only lateral heterogeneities in the lithospheric structure, which may include local weak zones or a pre-existing variations in lithosphere thickness, are expected to have an impact on lithospheric thinning and associated surface heat flux. In addition, anisotropic lithosphere viscosity may affect the position and dynamics of dripping instabilities (Király et al., 2020;Lev & Hager, 2008).
Finally, we can apply the identified trends to examples on Earth. Based on visual estimates of lithosphere thinning for central Greenland of about 27% (or 40 km, Celli et al. (2021)) and for southern South Africa of about 47% (about 70 km, Celli et al. (2020)) for 150 km thick cratons, we would expect maximum relative heat flux anomalies of about 20% and 60%, respectively ( Figure 4). In fact, Martos et al. (2018) infer a relative steady state heat flux anomaly of about 25% (assuming a maximum value of 70 mW/m 2 and a reference value of 56 mW/m 2 ), which represents the expected maximum heat flux anomaly after full equilibration. Note, however, that due to the age of the plume-lithosphere interaction (about 50-100 Ma for Greenland [Martos et al., 2018], and 100-130 Ma for Africa [Celli et al., 2020]), the heat flux anomaly in Greenland is likely still increasing, while the heat flux anomaly in Africa should be approximately at its maximum value or already decreasing. Furthermore, inferred heat fluxes along the hotspot track in Greenland (Martos et al., 2018) require more lithosphere thinning in the southeastern part, consistent with seismic tomography (Celli et al., 2021).

Conclusions
The interaction between a plume and continental or cratonic lithosphere is usually less apparent than it is for oceanic lithosphere, since extrusive volcanism is less common. However, the plume still leaves a trace in form of a thinned lithosphere and a (potentially) increased surface heat flux. The emplacement of hot plume material rapidly thins the lithosphere locally, and increases the lithospheric temperature gradient over time. As a consequence, positions and amplitudes of the two anomalies are directly linked and follow the same path along the plate. However, while lithosphere thinning occurs throughout the plume-lithosphere interaction, and in most cases maximum thinning is observed within a few Myr after the plume has passed, the maximum heat flux anomaly occurs approximately 40-140 Myr later due to the slow process of heat conduction. This has to be taken into account when interpreting geophysical data, which provide only a snapshot in time. The extent of lithosphere thinning, and thus also the amplitude of surface heat flux anomalies, is most sensitive to the viscosity of the lower lithosphere and asthenosphere, with plume excess temperature and plate velocity (for moving plate cases) or plume life time (stationary plate cases) having secondary influence. The thicknesses of the lithosphere and asthenosphere only play a minor role.
The relation between relative surface heat flux and relative lithosphere thinning is independent of the chosen geometry, and can be approximated by an analytical expression for the time-dependent thermal structure of a stationary system. However, the lithosphere in dynamic models (and probably Earth) typically does not have time to fully equilibrate thermally, and thus does not achieve the maximum possible heat flux. In fact, most models exhibit anomalies analogous to about 40-100 Myr of stationary thermal evolution, allowing us to predict expected and potential heat flux anomalies if the lithospheric thinning is known. This spatio-temporal relation of plume-induced lithospheric thinning and associated surface heat flux has important implications for understanding the potential of geothermal energy sources and for estimating glacial dynamics in polar regions.

Data Availability Statement
All data used in this numerical modeling study can be reproduced using the information given in the text, the Text S2 in Supporting Information S1 and the parameters in Table S1 in Supporting Information S1. An overview of the reproducible results from the numerical simulations with given parameters as shown in Figure 4 is provided in Tables S2-S4 in Supporting Information S1. The generated data sets are reproducible using the given parameters and the software ASPECT v.2.2.0, which is openly available from Bangerth et al. (2020).