Abstract The notion of self-regulating mantle convection, in which heat loss from the surface is constantly adjusted to follow internal radiogenic heat production, has been popular for the past six decades since Urey first advocated the idea. Thanks to its intuitive appeal, this notion has pervaded the solid earth sciences in various forms, but approach to a self-regulating state critically depends on the relation between the thermal adjustment rate and mantle temperature. I show that, if the effect of mantle melting on viscosity is taken into account, the adjustment rate cannot be sufficiently high to achieve self-regulation, regardless of the style of mantle convection. The evolution of terrestrial planets is thus likely to be far from thermal equilibrium and be sensitive to the peculiarities of their formation histories. Chance factors in planetary formation are suggested to become more important for the evolution of planets that are more massive than Earth.

Keywords

Heat flow

thermal budget

planetary evolution

INTRODUCTION The possibility of self-regulating mantle convection arises from the following negative feedback, which was originally suggested to H. Urey by R. Revelle (1). If the surface heat flux of a planet is higher than its internal heat production, the planet cools down, lowering the vigor of convection and thus heat flux. Conversely, if heat flux is lower than heat production, the planet heats up to enhance the vigor of convection until the balance is achieved between surface heat loss and internal heat production. Because radiogenic heat production decreases steadily with time, the convecting mantle must be able to adjust its internal temperature sufficiently quickly. The thermal adjustment time scale becomes shorter as the sensitivity of heat flux to temperature change increases, and the sensitivity can be quite high if one considers the temperature dependency of mantle viscosity. Tozer (2) was the first to quantify the thermal adjustment time scale using temperature-dependent viscosity, and his estimate of ~200 million years to achieve self-regulation was deemed sufficiently short. The validity of “exact” thermal equilibrium was later questioned by Schubert et al. (3), but their theoretical modeling of Earth’s thermal evolution, along with a similar study by Davies (4), still indicated that the surface heat flux should follow closely with the temporal variation of internal heat generation. An approximate version of self-regulating mantle convection has since been regarded as the standard theory for the thermal evolution of Earth and other terrestrial planets (5–7). One important corollary of self-regulation is that the thermal evolution of a planet becomes virtually insensitive to its initial condition. Whether a planet had a hot or cold start, its subsequent evolution would quickly adjust its internal temperature so that the vigor of mantle convection is adequate for the amount of heat provided by radioactive decay. However, despite its prevailing popularity, the actual efficacy of self-regulated mantle convection is debatable, at least for the case of Earth, which is undoubtedly the best understood terrestrial planet. First, the amount of internal heat production in the mantle is determined by the concentration of radioactive isotopes, such as 238U and 232Th, and geochemical estimates on the composition of the primitive mantle (8, 9), after correcting for continental crust, suggest that internal heat production in the convecting mantle accounts for only ~10 to 30% of convective heat loss (10, 11). Estimates on the composition of the convecting mantle, based on the petrology of mid-ocean ridge basalts, point to an even lower radiogenic contribution (12, 13). Many geophysicists seemed to have speculated that the present-day surface heat flux might be unusually higher than a long-term average, resulting in the apparent discrepancy between heat flux and heat production (14, 15), but such possibility has been shown to be inconsistent with the history of global sea-level change (16). At the same time, geochemical estimates are always associated with important assumptions and nontrivial uncertainty (17, 18), and some geochemists have used the notion of self-regulation to deduce the presence of a hidden geochemical reservoir enriched in heat-producing isotopes (19, 20). Though not widely appreciated, most of the geochemical estimates on the primitive mantle composition implicitly assume that the whole mantle is well mixed (21); violating this assumption would lead us to an inconvenient situation in which radiogenic heat production is geochemically unconstrained. Unfortunately, a potentially more direct approach with geoneutrinos continues to suffer from substantial uncertainty (22), and a notion akin to self-regulating mantle convection is still used to derive a geophysical estimate on internal heat production (23). The present-day heat production of what is referred to as the “reference undepleted mantle” in the widely used textbook by Turcotte and Schubert (23, 24) is based on the assumption that the radiogenic heat production in the convecting mantle is responsible for 80% of the convective heat flux, as opposed to only ~10 to 30% as mentioned above. Although the degree of self-regulation is still unclear for the present-day Earth’s mantle, we are now in a better position to evaluate under what condition mantle convection can be self-regulated, owing to recent progress on the heat flow scaling of mantle convection (25–27). Here, I conduct a critical appraisal of this long-standing issue by clarifying the underlying physics and by providing numerical examples that demonstrate the generality of my analysis. In particular, I show that a system’s propensity to self-regularization can be succinctly evaluated by calculating one nondimensional number, and I discuss its application to the thermal evolution of terrestrial planets at large.

BASIC THEORY Previous theoretical considerations on self-regularization have been biased to the condition in which self-regulation is possible (4, 11, 28). Therefore, to treat a variety of heat flow scaling laws proposed in recent years, I start with a more general analysis of the relation between heat flow scaling and the degree of self-regulation. Global energy balance for a convecting mantle may be expressed as (1)where c p and ρ are the specific heat of the mantle and its density, respectively; D is the depth of the mantle; T a is the average mantle temperature; h 0 is the initial radiogenic heat production per unit mass; λ e is the effective decay constant that best approximates the contribution of all of major radioactive isotopes; and q(T a ) is the surface heat flux per unit area, expressed as a function of T a . Although basal (core) heat flux is an important part of the global thermal budget (29), it is assumed to be zero for the following reasons. First, it does not play a role in the classical negative feedback mechanism for self-regulation. Second, the thermal evolution of the core, which affects the temporal evolution of core heat flux, is a complex problem itself (30, 31), especially given the recent development in the energetics of geodynamo (32). Third, using a few reasonable assumptions (10), the energy balance equations for the mantle and the core can be combined to a single equation similar to Eq. 1, with the mantle heat capacity term on the left-hand side (c p ρD) replaced with the whole-Earth heat capacity. Considering nonzero core heat flux increases the thermal inertia of the convecting system, helping it to deviate from the self-regulated state. The assumption of zero basal heat flux is therefore equivalent to the ideal condition to achieve self-regulation. This ideal setup will strengthen the main conclusion of this paper. The above equation can be nondimensionalized as (2)where temperature, time, heat production, and heat flux are normalized by ΔT, D2/κ, kΔT/(ρD2), and kΔT/D, respectively; ΔT is the temperature scale, κ is the thermal diffusivity [= k/(ρc p )], and k is the thermal conductivity. By first-order Taylor expansion, surface heat flux may be approximated as (3)where is the heat flux at an arbitrary reference temperature, . The (nondimensional) thermal adjustment rate, as introduced by Tozer (2), is given by the temperature derivative of heat flux, dq*/dT a *, which is also denoted here by the symbol γ* for short, and the ratio of the thermal adjustment rate over the decay constant, γ*/λ e *, will be called the Tozer number. With this linear approximation, Eq. 2 can be solved analytically, and in the case of initial surface heat flux being equal to initial heat production, that is, q*(T a *(0)) = h 0 *, the subsequent evolution may be expressed as (4)or in terms of heat flux as (5)where t′ is the reduced time defined as , is the internal temperature in equilibrium with internal heat production [= ], and Tz is the Tozer number. When Tz = 1, the right-hand side of Eq. 4 reduces to and that of Eq. 5 reduces to just t′. One reduced time corresponds to the e-folding time scale for internal heat production, and the Tozer number measures the rate of thermal adjustment with respect to that of radioactive decay. Starting with a different initial surface heat flux would merely introduce an additional term proportional to exp[−(Tz − 1)t′)]; thus, the above solutions capture the essence of temporal evolution. When Tz >> 1, Eq. 4 reduces to , that is, the perfect thermal equilibrium. The degree of self-regulation may be quantified by the Urey ratio, denoted by Ur (10, 33), which is defined as the contribution of internal heating to surface heat flux, that is, Ur(t) = h*(t)/q*(t), and Eq. 5 indicates that the Urey ratio converges to 1 − Tz−1 with an e-folding time scale of (γ* − λ e *)−1 when Tz > 1 and to zero with a time scale of (λ e * − γ*)− 1 when Tz < 1. The e-folding time scale may seem to diverge as Tz → 1, but when Tz = 1, the Urey ratio converges to zero with an e-folding time scale of (e − 1)/λ e *, which can serve as an upper bound on the time scale. Convection may be said to be self-regulated when Ur is ~1, but approach to such a condition is bounded by Tz−1 when Tz > 1, and self-regulation of any degree would become impossible when Tz ≤ 1. The debates over the present-day thermal budget of Earth mentioned in the Introduction may be summarized as follows: Geochemical models for the mantle composition indicate that Ur is ~0.1 to 0.3, but the standard geophysical models prefer that Ur be ~0.8, which is referred to here as an approximately self-regulated situation.

NUMERICAL EXAMPLES The generality of the above analysis may be understood by the following illustrative examples with fully dynamical systems. All cases considered here were run for a period of one reduced time, during which initial heat production decreases by a factor of e. The first example is with stagnant lid convection (Fig. 1). Stagnant lid convection is one type of thermal convection, and it occurs when viscosity is strongly temperature-dependent (34). The viscosity of the top thermal boundary layer is too high for the layer to appreciably deform, and convection takes place only beneath the boundary layer. Because the viscosity of silicate rocks, which constitute the mantle of terrestrial planets, is known to depend strongly on temperature (35), stagnant lid convection may be regarded as the most natural form of mantle convection, and Venus and Mars are in fact considered to be in this mode of convection (5). Stagnant lid convection, in its simplest form, requires only two nondimensional parameters, the internal Rayleigh number and the Frank-Kamenetskii parameter; the former indicates the vigor of convection, whereas the latter measures the temperature dependency of viscosity. For the example shown in Fig. 1, the internal Rayleigh number is set to 3 × 108 and the Frank-Kamenetskii parameter is set to 15. Both values are high enough to approach situations expected in actual planetary mantles (see Materials and Methods for further details). For purely internal heating (that is, no basal heat flux), time-averaged surface heat flux should be equal to internal heat generation at a statistically steady state or a quasi–steady state, and theoretical scaling exists for the relation between the internal temperature and the surface heat flux (25). For the cases shown in Fig. 1A, I use a quasi–steady-state solution corresponding to a q* of 7 (with T a * ~ 0.85) as the initial condition, and for those in Fig. 1B, I use a solution corresponding to a q* of 5 (with T a * ~ 0.76). The thermal adjustment rate γ* is ~17 at T a * ~ 0.8 (Fig. 1C). For both cases, internal heating starts with an h 0 * of 6, and is varied to generate cases with a Tz of 1 to 30. These two initial conditions are used to consider the approach to thermal equilibrium from both above and below. Fig. 1 Numerical simulation results for stagnant lid convection with decaying internal heat production. (A) The evolution of average temperature as a function of the reduced time t′, when starting with the initial temperature higher than the equilibrium temperature. The equilibrium temperature (black dashed line) is the temperature corresponding to the internal heat production according to the theoretical heat flow scaling law. Four cases with different (initial) Tz values are shown: 1 (blue), 3 (green), 10 (red), and 30 (magenta). Solutions from parameterized convection models, that is, solutions of Eq. 2 with the theoretical scaling are shown as colored dashed lines. (B) Same as (A), but starting with the initial temperature lower than the equilibrium temperature. (C) Covariation of average temperature and surface heat flux from the runs shown in (A) and (B). All cases are plotted right on top of each other and also on the theoretical scaling (black dashed line). Results from steady-state solutions are shown as solid squares (see Materials and Methods). Their temporal fluctuations are smaller than the size of the symbol. Note that in (A) and (B), the temperature deviation for the case of Tz = 30 starts to increase slightly after a t′ of ~0.6, because dq*/dT a * decreases at lower temperatures (C), making the actual Tz value lower than 30. The slight nonlinearity of the actual heat flow scaling results in small temporal variations in Tz. In Fig. 1 (A and B), the hypothetical internal temperature in perfect equilibrium with decaying internal heating is shown as a reference. Note that although internal heating follows exponential decay, this equilibrium temperature decreases almost linearly with time because of slight nonlinearity in the theoretical heat flow scaling (Fig. 1C). When Tz is as high as 30, thermal equilibrium is quickly achieved within a t′ of 0.1 in both cases. As Tz decreases, the degree of self-regulation deteriorates as expected. When starting from a low initial temperature (Fig. 1B), the system always passes through a thermal equilibrium at one point, but with low Tz, the system moves away from it because the Urey ratio should converge to 1 − Tz−1. Also shown in Fig. 1 (A and B) is the temperature evolution obtained by solving Eq. 2 with the theoretical heat flow scaling, which markedly agrees well with numerical simulation results, regardless of the value of Tz. This is because heat flux from numerical simulation closely follows the theoretical scaling of (quasi–)steady-state heat flux (Fig. 1C). With stagnant lid convection, only positive Tz can be tested. A more general dynamical system may be constructed with viscosity depending solely on average temperature, and by assuming that viscosity increases with higher temperature, steady-state heat flow scaling can exhibit even negative γ* (see Materials and Methods). The second example is made with such a system with a γ* of around −39 (Fig. 2C). Again, is varied to generate different Tz values. When Tz ≤ 1, the Urey ratio converges to zero, meaning that a system moves away from thermal equilibrium. More negative Tz leads to more rapid departure (Fig. 2, A and B), as predicted by the theoretical analysis. Surface heat flux from numerical simulation exhibits significant temporal fluctuations (Fig. 2C) because the top thermal boundary layer is not stabilized by higher viscosity, unlike in the previous case of stagnant lid convection. Note that these fluctuations are not due to the negative γ*; even with positive values of γ*, similar fluctuations take place with this dynamical system with average-temperature–dependent viscosity. Despite such fluctuations in surface heat flux, the temporal variation of internal temperature is smooth because of the thermal inertia of the convecting system. In any event, the relation between transient surface heat flux and average temperature still follows the theoretical steady-state scaling (Fig. 2C), and as a result, the solution of Eq. 2 provides a good approximation to the results of fully dynamic calculations (Fig. 2, A and B). Fig. 2 a *. Similar to Fig. 1 , but with a dynamical system characterized by a negative value of dq*/dT*. (A and B) The equilibrium temperature increases with time, whereas internal heating decreases with time, because of the negative dependency of surface heat flux on internal temperature. Three cases with different Tz values are shown: −1 (blue), −3 (green), and −10 (red). (C) Only the cases of Tz = −1 and −10 are shown for clarity. Solving Eq. 2 with steady-state heat flow scaling is also known as the parameterized convection approach (3, 4, 36), and the parameterized approach is valid as long as the assumed steady-state scaling can predict instantaneous heat flow reasonably accurately. As seen in the above examples, as well as in previous studies (37–39), this condition is usually satisfied for vigorously convecting systems, in which dynamical processes associated with the top thermal boundary layer take place more rapidly than the cooling of the entire system (36). Recently, Moore and Lenardic (40) suggested that parameterized convection models would be valid only when adjustment to thermal equilibrium could take place rapidly, that is, Tz >> 1. Their reasoning appears to stem from the misconception that the use of steady-state heat flow scaling requires a convecting system to reach thermal equilibrium sufficiently quickly. As seen in Figs. 1 and 2, the parameterized convection approach is valid for a wide range of Tz, including negative values. However, its applicability to more complex convection systems is still open to question, and future studies are warranted for this important issue.

DISCUSSION Figure 3A compares a few different heat flow scaling laws for Earth’s mantle, and corresponding Tozer numbers [with an effective decay constant of 2.5 billion years (Gy)] are shown in Fig. 3B (see Materials and Methods). The “classical” scaling, which predicts higher heat flux for a hotter mantle, is what has long been used in parameterized convection models (5). The “plate tectonics” scaling incorporates the effect of mantle melting (27); a hotter mantle starts to melt deeper beneath mid-ocean ridges, forming a thicker, dehydrated lithosphere, which is more difficult to deform at subduction zones and thus slows down plate motion (41). Also shown in Fig. 3A are scaling laws for (hypothetical) stagnant lid convection on Earth, using the same mantle properties assumed for the above scaling laws, with and without the effect of mantle melting (25, 26). Heat transport is much less efficient in stagnant lid convection than in plate tectonics because of the lack of subduction in the former. The Tozer number reflects the absolute value of heat flux as well as its gradient with respect to temperature, and whereas the classical scaling predicts Tz > 10 at T p > 1450°C, other scaling laws predict |Tz| < ~4 at T p < 1600°C. According to the plate tectonics scaling, surface heat flux is relatively insensitive to the mantle temperature, which also means that past plate tectonics was unlikely to be faster than the present one. There are a variety of observations in favor of this scaling, including the history of passive margins (42), the cooling history of the upper mantle (43), Earth’s degassing history inferred from the atmospheric composition (44), and the paleomagnetic reconstruction of continental motion (45, 46). The apparent lack of self-regulation in the present-day Earth’s mantle, as indicated by the geochemical thermal budget, can be understood as the consequence of low Tz values associated with the plate tectonics scaling. As mentioned in Basic Theory, the degree of self-regulation can be lowered further by considering basal heat flux (30). Fig. 3 Heat flux scaling and the Tozer number. (A) Some representative heat flow scaling laws for convection in Earth’s mantle are shown as a function of mantle potential temperature. Solid curves denote calculations based on the scaling of plate tectonics with pseudoplastic rheology (27), with an activation energy of 300 kJ mol−1, a reference viscosity of 1019 Pa s at 1350°C, and a viscosity contrast due to dehydration stiffening of 1 (that is, no melting effect; red) and 100 (blue). The effective friction coefficient is set to 0.025 (red) and 0.02 (blue) to reproduce the present-day average convective heat flux [~75 mW m−2 (10)] at the present-day mantle potential temperature of 1350°C (51). The black dashed line represents the scaling for stagnant lid convection (25), using the same viscosity parameters as stated above, and the pink dashed line shows the effect of mantle melting (26). (B) The Tozer number variations corresponding to the scaling laws shown in (A). A half-life of 2.5 Gy is used to calculate the effective decay constant. The potential temperature T p is equated with the average temperature T a , neglecting the effect of a thin, cold thermal boundary layer. The e-folding time scale for approach to or departure from thermal equilibrium is also denoted on the right. (C) The Tozer number variations as a function of planetary mass. Earth mass is used here as a reference. Different curves represent different heat flow scaling laws [classical, plate tectonics (PT), or stagnant lid] or different potential temperatures, as indicated by the labels. Mantle depth, mantle density, and surface gravity, which are used to calculate heat flux, are varied with planetary mass (50). The chemical composition of the mantle is assumed to be the same as that of Earth’s mantle; thus, the effective decay constant is the same as in (B). Similarly, the viscosity parameters are the same as those used in (A) and are assumed to be insensitive to planetary mass. However, even with the same chemical composition, more massive planets may give rise to higher average viscosity because of a pressure effect (52), which would lower surface heat flux; the decreasing trend of Tz with increasing planetary mass would thus be enhanced further. Using dimensional quantities, the Tozer number may be expressed as (6)Because surface heat flux is primarily controlled by the instability of the top thermal boundary layer, it is relatively insensitive to the total depth of the mantle D (47), and so is dq/dT a . Therefore, with other parameters being equal, the Tozer number decreases in general as the size of a planet increases (Fig. 3C). This trend is most clearly observed in the classical scaling. The variation of Tz with planetary mass is more complex for the plate tectonics scaling because the temperature at which the effect of mantle melting becomes significant also varies with planetary mass. However, the overall effect of increasing planetary mass remains to reduce the magnitude of the Tozer number, and if plate tectonics on Earth cannot achieve thermal equilibrium, then it would be more unlikely for super-Earths. Deviation from thermal equilibrium would be even more pronounced for the case of stagnant lid convection. A lower Tozer number also means a longer e-folding time scale (Fig. 3C), indicating that how a planet forms in the first few tens of million years could have a profound impact on its subsequent evolution over a few billion years. Parameterized convection models with the effect of mantle melting suggest that the influence of initial conditions on present-day observables is significant even for planets smaller than Earth (48). On the basis of self-regulated mantle convection, terrestrial planets have long been believed to have only short-term “thermal” memory. This has been the case even after the notion of the exact thermal equilibrium was called into question in the early 1980s (3). However, in the absence of self-regulation, chance factors in planetary formation, for example, the details of late-stage giant impacts, could lead to diverse planetary evolution.

MATERIALS AND METHODS Numerical simulations The nondimensionalized governing equations for thermal convection of an incompressible fluid with internal heat production was solved using the finite element method (49). To avoid wall effects, the aspect ratio of the convection models was set to 4, and the model domain was discretized with 256 × 64 uniform, two-dimensional quadrilateral elements. The surface temperature was fixed to zero, and the bottom boundary was insulated. The top and bottom boundaries were free-slip, and a reflecting boundary condition was applied to the side boundaries. Stagnant lid convection runs Before runs with decaying internal heat production, five runs were conducted with time-independent heat production to obtain steady-state solutions (or, more precisely, “quasi–steady-state” solutions because they are time-dependent) as well as to examine the magnitude of their temporal fluctuations. The Frank-Kamenetskii parameter was set to 15, equivalent to the activation energy of ~250 kJ mol−1 when scaled to Earth’s mantle condition, and the internal Rayleigh number was set to 3 × 108. Both parameter values are a priori nominal values, corresponding to the nondimensional temperature of 1.0. For convection with internal heating, the steady-state internal temperature is found out only after a convection model is run, and when the internal temperature is lower than 1.0, the internal Rayleigh number and the Frank-Kamenetskii parameter become lower accordingly (26). The five runs were started with a uniform temperature of 1.0, with internal heat production ranging from 3 to 7, and were run until t* = 1. Results shown in Fig. 1C were taken from the period of t* = 0.5 to 1. The theoretical heat flow scaling is based on the thick lid formula of Solomatov and Moresi (25). It can be shown that the internal temperature used in the formula and the average temperature are related as (7)where δ is the nondimensional thickness of the stagnant lid, which is given by , and Nu is the Nusselt number. The runs with decaying heat source shown in Fig. 1 (A and B) started with steady-state solutions for the internal heat production of 7 and 5, respectively, and the initial value of the decaying heat source was set to 6. The effective decay constant was calculated for a range of Tz, with γ* = 16.94. A dynamical system with average-temperature–dependent viscosity A convecting system with an arbitrary value of dq*/dT a * may be built on the scaling of isoviscous convection. Steady-state solutions were first obtained with internal Rayleigh numbers Ra 0 (defined with T* = 1) of 106, 3 × 106, and 107, and with time-independent internal heat production ranging from 10 to 20. From these solutions, the following heat flow scaling was obtained (8)where a = 0.264 and β = 0.336. The Nusselt number Nu here is defined as q*/ΔT, where the temperature scale, ΔT, is determined from the maximum temperature observed in the model domain, and Ra i is defined as Ra 0 ΔT. The temperature scale and the average temperature were found to be related as . If we assume that viscosity depends on the average temperature as , where θ is a measure of temperature dependency, the dynamics of such a system at any average temperature is identical to that of an isoviscous system with the same value of viscosity; thus, surface heat flux may be expressed as a function of average temperature as (9)For the example shown in Fig. 2, θ = −10 and Ra 0 = 105 were used. Initial conditions for the runs with decaying heat source were taken from quasi–steady-state solutions for isoviscous models: Ra 0 = 106.09 and h* = 25.5 for Fig. 2A and Ra 0 = 106.52 and h* = 29.4 for Fig. 2B. The effective decay constant is calculated for a range of negative Tz values, with γ* = −39.26. Heat flow scaling laws for mantle convection Among the four kinds of heat flow scaling laws shown in Fig. 3, the one for stagnant lid convection without the effect of mantle melting has been established by Solomatov and Moresi (25), and it is the simplest to reproduce. At the moment, the effect of mantle melting can be introduced only through local stability analysis (26), and the closed-form expression is not yet available [see figure 13 of Korenaga’s work (26) for the cases of Mars, Venus, and super-Earth]. The remaining two scaling laws (classical and plate tectonics) were built from a scaling law for mantle convection with pseudoplastic rheology (27), which is summarized in the following. The input parameters include the thermal conductivity (k), the mantle depth (D), the mantle density (ρ), the thermal expansivity (α), the thermal diffusivity (κ), the gravitational acceleration (g), the potential temperature contrast across the top thermal boundary layer (ΔT), the effective friction coefficient (μ), the viscosity at a reference temperature (η 0 = η(T 0 )), the activation energy for temperature-dependent viscosity (E), the viscosity increase owing to melting-induced dehydration (Δη), and the upper mantle density (ρ UM ). The parameters other than ρ UM represent mantle-wide averages. Surface heat flow is then given as (10)where Nu is the Nusselt number calculated as (11)Here, Ra i is the internal Rayleigh number, Ra c is the critical Rayleigh number, and Δη L,e is the effective viscosity contrast for lithosphere, which includes the effects of pseudoplastic rheology and mantle melting. The critical Rayleigh number is set to 103, and the internal Rayleigh number is defined as (12)and the internal viscosity, η i , is calculated as (13)where R is the universal gas constant and T s is the surface temperature (set to 273 K in all examples here). The effective viscosity contrast is calculated as (14)where Δη L,r is the reference viscosity contrast (without the effect of mantle melting), is the (nondimensional) thickness of the dehydrated lithosphere, and is the (nondimensional) thickness of reference thermal boundary layer. The reference viscosity contrast is calculated as (15)where ξ = μ/(αΔT), and θ is the Frank-Kamenetskii parameter given by (16)The thickness of the dehydrated lithosphere is determined as (17)where z 0 is the initial depth of mantle melting calculated as (18)and p 0 is the initial pressure of mantle melting (= (T s + ΔT − 1423) × 107 [Pa]). The thickness of the reference thermal boundary layer is given simply by 1/Nu r , where (19) For Fig. 3 (A and B), the following values were used: k = 4 W m−1 K−1, D = 2900 km, ρ = 4000 kg m−3, α = 2 × 10− 5 K−1, κ = 10−6 m2 s−1, g = 9.8 m s−2, η 0 = 1019 Pa s at a T 0 of 1350°C, E = 300 kJ mol−1, and ρ UM = 3300 kg m−3. For the classical scaling, Δη was set to unity, that is, no mantle melting effect, and μ was set to ~0.02 to reproduce a q of 74.5 mW m−2 at a T 0 of 1350°C. For the plate tectonics scaling, Δη was set to 100, and μ was set to ~0.025. The specific heat c p can be calculated as k/(ρκ), which is 103 J kg−1 K−1. For Fig. 3C, the mantle depth D, the gravitational acceleration g, and the mantle density ρ were assumed to be scaled as M0.3, M0.46, and M0.19, respectively, where M is the planetary mass, based on the work of Valencia et al. (50), and all other parameters were assumed to be the same as in the above.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

Acknowledgments: I would like to thank two anonymous reviewers for their constructive comments and suggestions that helped improve the clarity of the article. Funding: This article is based on work supported by the NASA through the NASA Astrobiology Institute under Cooperative Agreement No. NNA15BB03A issued through the Science Mission Directorate. Competing interest: The author declares that he has no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper. Additional data related to this paper may be requested from the author.