{'fraction_non_alphanumeric': 0.034099739452392926, 'fraction_numerical': 0.020203867075010283, 'mean_word_length': 4.958961321953877, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Terrestrial Planet Evolution in the Stagnant-Lid Regime: Size Effects and the Formation of Self-Destabilizing Crust

The ongoing discovery of terrestrial exoplanets accentuates the importance of studying planetary evolution for a wide range of initial conditions. We perform thermal evolution simulations for generic terrestrial planets with masses ranging from that of Mars to 10 Earth-masses in the stagnant-lid regime, the most natural mode of convection with strongly temperature- dependent viscosity. Given considerable uncertainty surrounding the dependency of mantle rheology on pressure, we choose to focus on the end-member case of pressure-independent potential viscosity, where viscosity does not change with depth along an adiabatic temperature gradient. We employ principal component analysis and linear regression to capture the first-order systematics of possible evolutionary scenarios from a large number of simulation runs. With increased planetary mass, crustal thickness and the degree of mantle processing are both predicted to decrease, and such size effects can also be derived with simple scaling analyses. The likelihood of plate tectonics is quantified using a mantle rheology that takes into account both ductile and brittle deformation mechanisms. Confirming earlier scaling analyses, the effects of lithosphere hydration dominate the effects of planetary mass. The possibility of basalt-eclogite phase transition in the planetary crust is found to increase with planetary mass, and we suggest that massive terrestrial planets may escape the stagnant-lid regime through the formation of a self-destabilizing dense eclogite layer.


Introduction
Plate tectonics is only observed on Earth and is likely important to Earth's uniquely clement surface conditions (e.g., Kasting and Catling, 2003). Other terrestrial planets in the Solar System (i.e., Mercury, Mars, and Venus) are generally considered to feature a rigid spherical shell encompassing the entire planet, with hot mantle convecting beneath the shell (e.g., Schubert et al., 2001). This mode of mantle convection is known as stagnant-lid convection. In fact, stagnant-lid convection may be most natural for planetary mantles because the viscosity of constituent materials is strongly temperature-dependent (Solomatov, 1995). The discovery of many extrasolar terrestrial planets with mass 1 to 10M ⊕ (e.g., Rivera et al., 2005;Udry et al., 2007;Queloz et al., 2009;Mayor et al., 2009;Léger et al., 2009;Charbonneau et al., 2009;Borucki et al., 2011) makes understanding planetary evolution in the stagnant-lid regime especially critical.
Parametrized models of stagnant-lid convection have long been applied to planets in our Solar System in an effort to infer likely planetary evolution scenarios from limited observational constraints (e.g., Stevenson et al., 1983;Spohn, 1991;Hauck and Phillips, 2002;Fraeman and Korenaga, 2010). Previous studies of massive terrestrial planets are more theoretical in nature, focusing on two broad questions. First, the effects of planetary mass on the likelihood of plate tectonics have been studied through scaling analyses and simple parametrized convection models (Valencia et al., 2007;O'Neill and Lenardic, 2007;Korenaga, 2010a; van Heck and Tackley, 2011). Second, the evolution of planets in the stagnant-lid regime has been contrasted with evolution with plate tectonics in the hope of identifying atmospheric signatures that would indicate the regime of mantle convection for a distant planet (e.g., Kite et al., 2009). Mantle dynamics in the stagnant-lid regime, however, can be more complex than previously thought owing to the effects of mantle processing and crustal formation, and the scaling law of stagnant-lid convection that takes such complications into account has been developed only recently (Korenaga, 2009). It is thus warranted to take a fresh look at the fate of massive terrestrial planets in the stagnant-lid regime and to explore the general effects of initial conditions including planetary mass.
This study extends a parametrized model of stagnantlid convection recently applied to Mars (Fraeman and Korenaga, 2010) to terrestrial planets of various masses, including massive planets that evolve in the stagnantlid regime that are termed "super-Venus" planets. This model incorporates the effects of compositional buoyancy and dehydration stiffening on mantle dynamics (Korenaga, 2009), which are rarely accounted for ex-cept in simulations of plate tectonics. Unlike in previous studies, sensitivity analyses are extensively performed to quantify the relationship between initial conditions and modeling results. Principal component analysis is used to simplify the interpretation of a large number of simulation results. Simple scaling analyses are also conducted to derive a theoretical basis for major modeling results. Moreover, the likelihood of plate tectonics is quantified by tracking the viscosity contrast across the lithosphere during each simulation.
The purpose of this study is to investigate paths along which generic terrestrial planets may evolve and to estimate whether massive terrestrial planets are relatively more or less likely to escape the stagnant-lid regime. Throughout this paper, "Mars" and "Venus" should be considered shorthand for generalized 0.107M ⊕ and 0.815M ⊕ terrestrial planets, respectively. The evolution of a particular planet is likely to diverge from the predictions of these simple parametrized models. Few constraints are available beyond planetary mass and radius for extrasolar terrestrial planets. But for terrestrial planets in our Solar System, more data are available from decades of observations and spacecraft visits. Here, we explore hypothetical planetary evolution with the simplest assumptions on mantle dynamics, thereby serving as a reference model on which additional complications may be considered if necessary.

Theoretical Formulation
Parametrized convection models are used to simulate the evolution of Mars, Venus, and putative super-Venus planets for a wide range of initial conditions. Equations used to track the thermal and chemical evolution of terrestrial planets are taken from Fraeman and Korenaga (2010) with some modifications. Earth-like, peridotite mantle compositions are used to parametrize melting behavior. Although continuous evolution in the stagnant-lid regime is assumed, a simple model of lithospheric weakening is also considered to evaluate the likelihood of plate tectonics occurring at some point during planetary evolution.

Governing Equations
Mars, Venus, and super-Venus planets are assumed to begin as differentiated bodies with a mantle and core. Energy conservation yields two governing equations. First, the energy balance for the core is (Stevenson et al., 1983) [4πR 2 i ρ c (L c +E g ) where R c and R i are the radii of the core and inner core, respectively; ρ c is the density of the core; L c is the latent heat of solidification associated with the inner core; E g is the gravitational energy liberated per unit mass of the inner core; η c is the ratio of T cm , the temperature at the core side of the core/mantle boundary, to the average core temperature; C c is the specific heat of the core; and F c is the heat flux out of the core. The formulation of core cooling is identical to that of Stevenson et al. (1983).
Second, the energy balance for the mantle is (Hauck and Phillips, 2002) 4π where R m is the radius of the mantle; Q m is the volumetric heat production of the mantle; ρ m is the density of the mantle; C m is the specific heat of the mantle; η m is the ratio of the average temperature of the mantle to T u , the potential temperature of the mantle (a hypothetical temperature of the mantle adiabatically brought up to the surface without melting); f m is volumetric melt production with associated latent heat release, L m ; and F m is the heat flux across the mantle/crust boundary. Some of the above parameters are universal constants, but most are planet-specific. Many important parameters are also time-varying. In particular, mantle melt is extracted to form crust, causing R m to decrease with time. Likewise, Q m decreases with time because of radioactive decay with some approximated average decay constant, λ (Stevenson et al., 1983), and extraction through mantle processing. Figure 1 illustrates the assumed thermal and chemical structure in our model. Over time, melting processes an upper region of the original primitive mantle (PM) to form the crust and the depleted mantle lithosphere (DML). In parallel, the mantle lithosphere (ML), which is always thicker than the DML, develops as a conductive thermal boundary layer underlying the crust. As part of the DML can potentially delaminate and be mixed with the convecting mantle, the composition of the convecting mantle can be more depleted than that of the PM. The mantle below the DML is thus referred to as the source mantle (SM), the composition of which is initially identical to the composition of the PM but can deviate with time. The history of these layers strongly depends on convective vigor, effects of mantle melting, and initial conditions.

Stagnant-Lid Convection with Mantle Melting
Standard parameterizations are used for mantle rheology and the vigor of convection. Mantle viscosity is a function of mantle potential temperature and the degree of hydration as (Fraeman and Korenaga, 2010) where A is a constant factor calculated using a reference viscosity η 0 at a reference temperature T * u = 1573 K; E is the activation energy; R is the universal gas constant; and ∆η w is the viscosity contrast between wet  Figure 1: Cartoons showing the assumed thermal and chemical structure of terrestrial planets taken from Fraeman and Korenaga (2010). In general, terrestrial planets are divided into a crust, mantle, and core, as shown in the left panel. Mantle that has been processed by melting and stays in the thermal boundary layer is depleted mantle lithosphere (DML). The thickness of the DML must always be equal to or less than the thickness of the mantle lithosphere (ML). The section of the mantle below the thermal boundary layer is the sublithospheric mantle. The right panel shows the horizontally-averaged temperature distribution. Key model parameters are also indicated.
and dry mantle. We use an activation energy of E = 300 kJ mol −1 , which is appropriate for diffusion creep and dislocation creep within a Newtonian approximation (Christensen, 1984;Karato and Wu, 1993;Korenaga, 2006). The normalized water concentration in the source mantle, C W SM , has an initial value of one and decreases toward zero as mantle melting causes dehydration. To write Eq. 3, we make the major assumption that mantle viscosity is not strongly pressuredependent, which is consistent with early studies of the evolution of large rocky planets (e.g., Valencia et al., 2006) and some theoretical predictions (Karato, 2011) but in contrast to recent work (e.g., Papuc and Davies, 2008;Stamenkovic et al., 2011Stamenkovic et al., , 2012. Our choice is thus further explained in the discussion section. Two non-dimensional parameters characterize thermal convection with the above viscosity formulation (Solomatov, 1995). First, the internal Rayleigh number serves to quantify potential convective vigor (Fraeman and Korenaga, 2010) where α is the coefficient of thermal expansion; κ is the thermal diffusivity; T c and T u are, respectively, the temperature at the bottom of the crust (called "Moho tem-perature") and the mantle potential temperature defined at the top of the mantle; and h m is the thickness of the mantle. Second, the Frank-Kamenetskii parameter is defined as (Solomatov, 1995;Fraeman and Korenaga, 2010) With these two parameters, the average convective velocity beneath the stagnant-lid may be calculated as (Solomatov and Moresi, 2000) To include the effects of compositional buoyancy and dehydration stiffening, the Nusselt number, which is a non-dimensional measure of convective heat flux, must be calculated with a local stability analysis at each time step (Korenaga, 2009). The symbolic functionality may be expressed as where ∆ρ and ∆η m are the density and viscosity contrasts between the source mantle and depleted mantle, respectively, and h l is the thickness of the depleted man-tle lithosphere. The thickness of a thermal boundary layer in the mantle is then easily calculated using The chemical evolution of the mantle strongly affects terrestrial planet evolution. To first order, partial melting of the mantle can be considered to begin at a depth where the temperature exceeds the solidus of dry peridotite, as long as the mantle is not significantly wet (Hirth and Kohlstedt, 1996). The initial pressure of melting is (Korenaga et al., 2002) where (dT /dP ) S is the adiabatic mantle gradient, which is roughly constant for the pressure range relevant to mantle melting. Therefore, P i should be approximately constant for any terrestrial planet with Earth-like mantle composition. Melting stops when the convective upwelling reaches the base of the mantle lithosphere. That is, the final pressure of melting is given by where h c is the thicknesses of the crust; g is gravitational acceleration; and ρ L is the density of the lithosphere. For convenience, we use the Martian ρ m as ρ L , noting that ρ L should remain roughly constant whereas ρ m , an averaged mantle parameter, increases with planetary mass because of pressure effects. If P f < P i , then melting occurs in the melting zone between P i and P f , with thickness d m and average melt fraction equal to where (dφ/dP ) S is the melt productivity by adiabatic decompression. Volumetric melt productivity is finally parametrized as where χ ∼ 1 if the upwelling mantle is cylindrical and all downwelling occurs at the cylinder's edge (Solomatov and Moresi, 2000;Fraeman and Korenaga, 2010). The crustal temperature profile is calculated as in Fraeman and , with the modification that crustal material with a temperature above T crit = 1273 K is considered to be buoyant melt that migrates within one time step immediately below the planet's surface, producing a relatively cooler crust and a larger mantle heat flux. This modification is not important for Martian cases, because the Moho temperature does not reach the threshold except for some extreme cases, but becomes essential to achieve a realistic crustal thermal profile for larger planets.

Likelihood of Plate Tectonics
Thermal evolution models featuring stagnant-lid convection are not applicable to planets on which plate tectonics occurs. If a suitable weakening mechanism exists, the lithosphere may be broken into plates and recycled into the mantle. Many aspects of plate tectonics on Earth, however, are not captured in current mathematical models (Bercovici, 2003). Quantifying the conditions under which plate tectonics is favored over stagnant-lid convection is likewise difficult, and the effect of planetary mass on the likelihood of plate tectonics has been controversial (Valencia et al., 2007;Korenaga, 2010a; van Heck and Tackley, 2011). Recent studies suggest, however, that the effects of planetary mass on yield and convective stresses may be dominated by uncertainties in other important planetary parameters, such as internal heating and lithosphere hydration (Korenaga, 2010a;van Heck and Tackley, 2011).
This study uses a simple scaling that is consistent with current understanding of rock mechanics (Korenaga, 2010a), though the possibility of different lithosphere weakening mechanisms (e.g., Landuyt et al., 2008) cannot be excluded. We assume that plate tectonics can occur if convective stress exceeds the brittle strength of lithosphere given by where c 0 is the cohesive strength, µ is the effective friction coefficient, and z is depth (Moresi and Solomatov, 1998). Experimental data indicate that the cohesive strength is negligible under lithospheric conditions, i.e., c 0 /(µρz) 1 (Byerlee, 1978). We use another nondimensional parameter (Korenaga, 2010a): where the relevant temperature difference is the difference between the mantle potential and surface temperatures. In the parameterized convection model formulated in the previous section, we separately consider the crust and the mantle, but when discussing the likelihood of plate tectonics using the scaling of Korenaga (2010b), it is more convenient to treat the crust and mantle together, assuming that crustal rheology is similar to mantle rheology. Detailed scaling analyses (Korenaga, 2010a,b) show that the effective viscosity contrast across the lithosphere can be parameterized as where θ tot is the Frank-Kamenetskii parameter defined using the total temperature difference explained above, i.e., A transition from plate-tectonic to stagnant-lid convection can take place if the above viscosity contrast ex-ceeds a critical value ∆η L,crit = 0.25Ra where Ra i,tot is defined to incorporate surface temperature as For each simulation, if ∆η L /∆η L,crit ≤ 1 at any time, then plate tectonics may have been favored at some point during the evolution of a given planet. The satisfaction of this criterion may strongly depend on the value of µ, so a wide range of values should be tested. For silicate rocks, plausible values of µ range from 0.6 to 0.7, according to both laboratory studies (Byerlee, 1978) and measurements of crustal strength on Earth (e.g., Brudy et al., 1997). Surface water, however, may lower these values substantially via thermal cracking and mantle hydration (Korenaga, 2007).

Numerical Models
The parametrized model described above was used to calculate thermal histories of Mars-and Venus-like planets and 1 to 10M ⊕ super-Venus planets, where the ⊕ subscript denotes parameters for Earth, for a duration of 4.5 Gyr using numerical integration with a time step of 1 Myr. A wide parameter space was explored by varying the initial mantle potential temperature, T u (0); the initial core/mantle boundary temperature, T cm (0); the initial volumetric heat production, Q 0 ; the reference mantle viscosity, η 0 ; and the viscosity contrast between dry and wet mantle, ∆η w . Previous work for Mars demonstrated that simulation results were not very sensitive to the degree of compositional buoyancy and other parameters (Fraeman and Korenaga, 2010). Table 1 lists model constants common to all simulations, and Table 2 lists planet-specific ones.

Application to Mars-and Venus-like Planets
For Venus, the following sets of initial conditions were used: Initial mantle potential temperature, T u (0) = 1400, 1550, 1700, 1850, and 2000 K; initial core/mantle boundary temperature, T cm (0) = 3500, 4000, and 4500 K; reference viscosity, η 0 = 10 18 , 10 19 , and 10 20 Pa s; and dehydration stiffening, ∆η w = 1, 10, and 100. For Mars, initial core/mantle boundary temperatures were 2250, 2500, and 3000 K. In addition, five different values were tested for the amount of internal heating Q 0 . Compositional buoyancy was set as (dρ/dφ) = 120 kg m −3 for all simulations. The fraction of light elements in the core was fixed at 0.2 for all simulation runs to avoid inner core solidification (Schubert et al., 1992;Fraeman and Korenaga, 2010). Simulations were performed for all permutations of the above initial conditions, although unrealistic simulation results were discarded in the following way: For both Venus and Mars, inner core growth was disallowed and total surface heat flux at the present was required to be positive. Furthermore, the condition h c (t p ) < 500 km was imposed to disregard results with unrealistic crustal growth. None of the 675 simulations for Venus failed these criteria, but 30 of the 675 simulations for Mars were discarded. The appropriate magnitude of radiogenic heating is poorly constrained in general, especially for terrestrial exoplanets. Even for Earth, the abundance of radiogenic heating is controversial. Geochemical constraints support a low Urey ratio, the ratio of internal heat production to surface heat flux, but this is known to conflict with the cooling history of Earth unless a nonclassical heat-flow scaling for mantle convection is assumed (Christensen, 1985;Korenaga, 2008). A Urey ratio close to one has thus long been preferred from a geophysical perspective (Davies, 1980;Schubert et al., 1980Schubert et al., , 2001 and can be used to provide an upper bound for the initial concentration of radioactive elements in Earth's chemically undifferentiated mantle. Assuming a present-day surface heat flux of 46 TW (Jaupart et al., 2007), an extreme upper bound for Earth is Q 0 ≈ 3.5×10 −7 W m −3 .
A recent petrological estimate on the thermal history of Earth is actually shown to favor a low Urey ratio (∼0.3) with a non-classical heat-flow scaling (Herzberg et al., 2010), indicating that geochemical constraints on the heat budget may be robust. In the thermal evolution models of Kite et al. (2009), for example, concentrations of 40 K, 232 Th, 235 U, and 238 U taken from Ringwood (1991) and Turcotte and Schubert (2002) were considered, corresponding to values for Q 0 between 1.2 × 10 −7 and 8.7 × 10 −8 W m −3 for Venus. We thus chose to use the following values for initial volumetric radiogenic heating: Q 0 = 0.5, 0.75, 1.0, 1.25, and 1.75   (2010) and references therein. Venusian values can be found in Spohn (1991) and Noack et al. (2011). Super-Venus values were calculated in this study from simple interior models following Seager et al. (2007).
ets, Q 0 was multiplied by ρ m /ρ m,♀ , where ♀ indicates the Venusian value, to maintain constant element abundances in more or less compressed mantles.

Application to Super-Venus Planets
One-dimensional profiles of massive terrestrial exoplanets were generated to calculate planet-specific constants used in the above stagnant-lid convection model. Many interior structure models exist for massive solid exoplanets, ranging from simple to very complex (Valencia et al., 2006;Seager et al., 2007;Sotin et al., 2007;Wagner et al., 2011). To study the first-order effects of planetary mass on stagnant-lid convection, a relatively simple structure with an Fe( ) core and a MgSiO 3 mantle is assumed, as in Seager et al. (2007) and Kite et al. (2009). The resulting interior density and pressure distributions neglect several obvious factors such as temperature effects, but yield results remarkably similar to those from more complex models.
Three equations are solved to calculate m(r), the mass contained within radius r; P (r), the pressure distribution; and ρ(r), the density distribution. A selfconsistent internal structure must satisfy the material specific equations of state the equation of hydrostatic equilibrium and the conservation of mass equation for a spherical mass distribution where G is the gravitational constant; T (r) is the radial temperature profile; and f EOS represents a materialspecific equation of state (Seager et al., 2007). Equations of state are numerically calculated using constants from Table 3 to sufficient resolution so that ρ(r) can be determined to within ±1 kg m −3 . For a desired M P , equations (20) and (21)  The material equations of state and the equations of conservation of mass and hydrostatic equilibrium were numerically integrated to build simple planets, and the interior boundary condition was adjusted until the resulting planet had the desired mass. These simple models are used to calculate averaged values for mantle density, ρm; core density, ρc; surface gravity, g; central pressure, Pc; and core/mantle boundary pressure, Pcm. For comparison, the density distribution for Earth from the Preliminary Reference Earth Model (PREM) of Dziewonski and Anderson (1981) is also plotted.
boundary conditions M (0) = 0 and P (0) = P c , where P c is a guessed central pressure. The outer boundary condition is simply P (R P ) = 0. Errors associated with ignoring temperature effects are limited to a few percent (Seager et al., 2007). With this method, the choice of P c determines the R P at which the outer boundary condition is satisfied. These calculations are iterated with the bisection method until P c is found such that m(R P ) = M P to within 0.1%. The equation of state for Fe( ) is used until m(r) = 0.325M P , mandating a 32.5% core mass fraction. The MgSiO 3 perovskite to enstatite phase transition is assumed to occur at 23 GPa, although the transition pressure increases to ∼25 GPa at ∼800 K (e.g., Ita and Stixrude, 1992). Neglecting this phase transition would produce unrealistically high near-surface densities.
Pressure, mass, and density distributions were cal- 3.97 -0.016 4100 BME4 MgSiO 3 (en) 125 5 N/A 3220 BME3 culated for planets with M P = 1, 2, 4, 5, 6, 8, and 10M ⊕ . From these models, averaged densities for the core and mantle were calculated, along with P c and P cm . Surface gravitational accelerations were calculated using g = GM P /R 2 P for each planet. These constants are reported in Table 2. The density distributions for these planets are shown in Fig. 2, along with the density profile for Earth from the Preliminary Reference Earth Model (PREM) of Dziewonski and Anderson (1981). Compared to PREM, the scheme for calculating internal structures used in this study overestimates the density of the core and underestimates the radius of the core of an Earth-mass planet, but returns R P (M ⊕ ) ≈ R ⊕ despite ignoring details of mineral composition, phase transitions, and temperature effects. In more massive planets, the enstatite to perovskite phase transition occurs at much shallower depths because a higher surface gravity causes a greater increase in pressure with depth. Furthermore, P c increases much more rapidly than P cm with increasing planetary mass.
Most of the Martian and Venusian initial conditions can be used for super-Venus thermal evolution models, but some must be modified appropriately. For example, the core/mantle temperature for super-Venus planets should increase along the mantle adiabatic temperature gradient: for 5 and 10M ⊕ super-Venus planets, initial core/mantle boundary temperatures are increased by roughly 350 and 900 K, respectively, from the initial conditions for Venus. These temperatures still correspond to a so-called "hot start," which is likely for terrestrial planets because of the large magnitude of gravitational potential energy released during accretion (Stevenson et al., 1983). Only three simulations for the 5M ⊕ planet, and no simulations for the 10M ⊕ planet, failed the requirements on crustal thickness, surface heat flux, and inner core growth.

Results
Thermal evolution simulations were performed for Mars-and Venus-like planets and super-Venus planets. The following sections summarize the results, beginning with a few representative examples for Mars and Venus. Principal component analysis, as described in the appendix, was applied using all simulation results to identify major model behaviors. We also tried to quantify relations between input and output parameters and, despite the complexity of our model formulation, a linear function of initial conditions is found to reasonably approximate many output parameters of interest.

Sample Thermal Histories for Mars-and Venuslike Planets
Sample thermal histories for Mars and Venus are shown in Fig. 3. These models span the entire range of initial radiogenic heating values with all other initial conditions set to intermediate values. In particular, T u (0) = 1700 K, µ 0 = 10 19 Pa s, and ∆µ w = 100. For Venus and Mars, respectively, T cm (0) = 4000 and 2500 K and T s = 220 and 730 K. Initially very hot cores are assumed here because core segregation is expected to release a large amount of gravitational potential energy. This excess heat is released into the mantle during the first hundred million years of planetary evolution. Thereafter, mantle dynamics controls core cooling. Whereas internal heating has a great effect on surface heat flux, mantle temperatures only differ to within ±200 K for the sampled range of Q 0 . Mars evolves with a consistently lower potential temperature than Venus. Because Mars also has a relatively shallow mantle, the Martian core is cooled down more efficiently.
Crustal thickness is an important, potentially observable constraint for planetary evolution models. Mars and Venus, with different magnitudes of radiogenic heating, have very different crustal formation histories. Both planets start with no initial crust, but quickly produce some through mantle melting. For Venus, Moho temperatures quickly reach the melting point of basalt for all initial internal heating choices. Crustal production occurs for the first ∼1 Gyr of evolution, with thicker crust for higher internal heating. For Mars, crustal production is gradual and crustal temperatures are much lower, with increased internal heating causing an longer period of crustal formation and increased total crustal production. Both Mars and Venus undergo substantial mantle processing, indicating that the deep interior serves as a significant source of endogenous water, especially during the first ∼1.5 Gyr of their evolution.

Sensitivity Analyses for the Evolution of Mars-and Venus-like Planets
Figures 4 and 5 summarize the results of 1320 simulations for Venus and Mars, respectively. Present-day values for selected output parameters are plotted against crustal thickness for both planets. Several correlations Principal component analysis facilitates the interpretation of the correlations between output parameters. For Venus, two principal components account for most (>65%) of the variance of the planetary parameters after 4.5 Gyr of thermal evolution. Calculations of the principle components returns coefficients with values between -1 and 1 that are associated with each model output parameter. Appendix A contains a table of principal component basis vectors for Venus. Comparing numerical values of select coefficients may reveal correlations with physical explanations. Arrows representing the eigenvectors associated with these principal components are plotted in Figs. 4 and 5. These arrows indicate the axes along which the majority of the variance in the model output primarily lies. No preferred polarity exists for the principal component eigenvectors; plotting these arrows with a 180 • rotation would be equally valid.
The first principal component represents the most dominant correlations among present-day planetary parameters, which are characterized mainly by the thicknesses of the crust and mantle lithosphere layers, as they are associated with large coefficients: h c (0.34), h l (-0.36), and h M L (-0.39). Because the sign of the coefficient for h c is opposite to the sign of the other two coefficients, the thicknesses of the crust and mantle lithosphere are anti-correlated. In other words, thick crust is associated with thin depleted mantle lithosphere and a thin thermal boundary layer and vice versa, since principal components have no preferred polarity. An initially hotter mantle produces thicker crust and thicker depleted lithosphere. Because mantle viscosity is lower for hotter mantle, however, the depleted lithosphere is more likely to be destabilized, resulting in a thinner lithosphere (and thus thermal boundary layer) for thicker crust. Other coefficients in the first principal component indicate the effects of crustal thickness on other model parameters, including the first-order correlations observed during inspection of Fig. 4. For instance, thick crust is associated with high Moho temperature and high surface and mantle heat fluxes. Thick crust also indicates a high degree of mantle processing and a corresponding low present-day mantle water content. Finally, the large negative coefficients for both t c,10% and log 10 (∆t c,tot ) indicate that thick crust tends to form early and quickly.
The second principal component elucidates the effects of planet temperatures on other model parameters, since large coefficients are associated with T u (0.45) and T cm (0.44). Unsurprisingly, high mantle potential and core/mantle boundary temperatures are associated with high Moho temperature, since T c has a coefficient of 0.25. Moreover, high interior temperatures correspond to thick crust and a high degree of mantle processing, which would cause the present-day mantle water con- centration to be very low. Note, however, that high (present-day) interior temperatures do not correspond to thick crust in case of the first principal component (Fig. 4b). In this space of crustal thickness and upper mantle temperature, the first and second principal components are nearly orthogonal, thus explaining the overall spread of simulation results. With principal component analysis, we can visualize how the most dominant trend (represented by the first principal component) is affected by secondary factors and how these secondary factors manifest in different parameter spaces. An important point is that the overall variability of planetary evolution can be compactly represented by a small number of principal components; that is, the effective dimension of the model space is actually small.
The principal components for Mars are very similar to those for Venus, with some notable exceptions. The first principal component again represents the effects of strongly correlated Moho temperature and crustal thickness, and thus explains the largest portion of the variance in the model output. As for Venus, a thin, cold crust is associated with thick depleted mantle lithosphere, a thick thermal boundary layer, a low surface heat flux, and a low degree of mantle processing. Unlike Venus, however, the surface and mantle heat fluxes in the first principal component are anti-correlated (see also Fig. 5).
Despite the complexity of our thermal evolution model, some present-day parameters are found to be predicted with reasonable accuracy for Venus and Mars using a linear function. A general formula for this function is B i = A i,0 + A i,1 T u,n (0) + A i,2 T cm,n (0) + A i,3 (log 10 (η 0 )) n + A i,4 (log 10 (∆η w )) n + A i,5 Q 0,n , (22) where B i is the value of the desired output parameter after 4.5 Gyr, constants A i,0 through A i,5 are estimated using the least-squares method for each B i , and the subscript n indicates that the input parameters are normalized and mean subtracted. Table 4 lists constants for Venusian B i that have relatively high correlation coefficients between predicted and actual simulation results. Figure 6 shows contour plots with predicted values of mantle potential temperature, crustal thickness, and duration of crustal formation for given initial internal heating and mantle potential temperature. While present-day mantle potential temperature and crustal thickness depend strongly on both initial mantle potential temperature and the magnitude of internal heating, the total duration of crustal formation is primarily a function of initial mantle potential temperature (see Table 4 for more complete information on parameter sensitivity). Figure 6 also demonstrates a reasonable correspondence between the predicted and actual values of these model outputs for all of the simulations. This way of summarizing simulation results allows us not only to see the sensitivity of model outputs to initial parameters but also to quickly reproduce major modeling results without redoing simulation.

Evolution of Super-Venus Planets
We investigate the evolution of super-Venus planets to explore the effects of planetary mass on stagnant-lid convection. For simplicity, surface temperatures for all super-Venus planets are assumed to be 300 K, though in reality this temperature may vary with time and is highly dependent on atmospheric composition and on the distance to and the luminosity of the central star.

Sample Thermal Histories
Super-Venus planets with M P = 1, 5, and 10M ⊕ were evolved to study the effects of increasing planetary mass on a variety of parameters, particularly crustal production. For all three planets, Q 0 was scaled to the Venusian value of 1.0×10 −7 W m −3 , T u (0) = 1700 K, and T cm (0) = 4000, 4350, and 4900 K, respectively. Dehydration stiffening and compositional buoyancy were both incorporated as usual. Figure 7 shows the results of these simulations. As with Venus and Mars, the transient "hot start" in the core is lost in the first ∼100 Myr. After this initial cooling, mantle dynamics controls core cooling. Because the mantle heats up for the first ∼1 Gyr and then cools only very slowly, core cooling is precluded for the first ∼2 Gyr. As suggested by simple scaling laws (Stevenson, 2003), mantle cooling paths for massive super-Venus planets are roughly parallel. Figure 7 also shows how the thicknesses of the crust, mantle lithosphere, and depleted mantle lithosphere vary with time. With increasing planetary mass, crustal thickness decreases. The simple scaling analyses below indicate that more massive planets have greater melt production. The observed increase in mantle potential temperature with planetary mass only accentuates this effect. The increased melt volume, however, is not sufficient to create a thicker crust on a larger planet. The 1M ⊕ planet in the stagnant-lid regime ceases crustal production soon after 1 Gyr as mantle potential temperature drops below a critical value. The increased interior temperatures for the more massive planets allow longer durations of crustal production. For the first ∼2 Gyr of thermal evolution, the thickness of the depleted mantle lithosphere is close to that of the mantle thermal lithosphere, reflecting the continuous delamination of excess depleted mantle lithosphere. Decreased crustal production with increasing planetary mass corresponds to a smaller degree of mantle processing and a higher content of residual mantle water.

Sensitivity Analyses
The output of 1347 simulations for 5 and 10M ⊕ super-Venus planets are shown in Figs. 8 and 9. Three simulations for the 5M ⊕ super-Venus planet were excluded because they did not meet the requirements that h c < 500 km and that inner core growth did not occur. As for Mars and Venus, present-day parameters of  Figure 6: Top panels show predicted values for Venus after 4.5 Gyr given initial volumetric radiogenic heating and initial mantle potential temperature, the initial conditions to which the output is most sensitive. Bottom panels show the correspondence between predicted and actual simulation results. The dashed line represents perfect predictive power. Default initial conditions are Tcm(0) = 4000 K, η 0 = 10 19 Pa · s, ∆ηw = 100, and (dρ/dφ) = 120 kg m −3 . Panels show predicted values of (a) mantle potential temperature, (b) crustal thickness, and (c) total time for crust to grow from 10% to 95% of its present thickness.
interest are plotted against present-day crustal thickness. The principal component eigenvectors, explained below, are projected onto each plot, emanating from the average simulation output. The table in the appendix contains the principal component basis vectors for the 10M ⊕ planet.
These scatter plots reveal similarities between the evolution of both massive planets. For instance, Moho temperature increases with crustal thickness in a linear fashion before reaching the critical value for basalt melting. With increasing planetary mass, the critical crustal thickness at which this transition occurs decreases. For relatively thick crust, Moho temperatures remain near the critical value for basalt melting. For both super-Venus planets, an increase in crustal thickness is associated with an increase in present-day mantle potential temperature, mantle heat flux, and degree of mantle processing. The total duration of crustal formation decreases with increasing present-day crustal thickness. Again, correlations between model parameters may be studied in more detail with principal component analysis.
For both planets, as for Venus and Mars, the first principal component is characterized by a strong correlation between Moho temperature and crustal thickness, explaining the general trends observed in Figs. 8 and 9. A decrease in both quantities is associated with an increase in the thicknesses of the depleted mantle lithosphere and the thermal boundary layer, a decrease in surface and mantle heat fluxes, an increase in the duration of crustal formation, and a decrease in the degree of mantle processing. The second principal component illuminates the effect of correlated interior temperatures. As expected, increasing mantle potential and core/mantle boundary temperatures causes an increases in crustal thickness, the total duration of crustal formation, and the degree of mantle processing.
Many present-day model parameters of interest can be represented as a linear function of initial conditions (Table 5). Compared to the case of Venus, a greater number of parameters are found to be approximated reasonably well by this approach. The effects of melting at the base of the crust undoubtedly remain a large source of nonlinearity in the model output for all terrestrial planets more massive than Mars. A more elaborate numerical implementation to deal with exceedingly high crustal temperatures may reduce such nonlinearity, though we did not explore this possibility.

Scaling of Crustal Thickness and Mantle Processing
We conduct simple scaling analyses to better understand the cause of decreasing crustal thickness and a decreasing degree of mantle processing with increasing planetary mass.

Crustal Thickness
A number of parameters govern the scaling of crustal thickness with planetary mass. Increased melt production, for instance, is the first requirement for thicker  Table 4: Coefficients for the best-fit linear function (Eq. 22) relating parameter values after 4.5 Gyr for parameters with correlation coefficients > 0.70 to a given set of initial conditions for Venus. Correlation coefficients quantifying the correspondence between the actual and predicted output parameters were calculated using normalized and mean subtracted input and output parameters. The average values of the input parameters are Tu(0) = 1700 K, Tcm(0) = 4000 K, log 10 (η 0 ) = 19, log 10 (∆ηw) = 1, and Q 0 = 1.05 × 10 −7 W m −3 . For the best-fit function, the input parameters are mean subtracted and normalized by 212 K, 408 K, 0.82, 0.82, and 4.30 × 10 −8 W m −3 , respectively.  Table 5: Coefficients for the best-fit linear function (Eq. 22) relating parameter values after 4.5 Gyr for parameters with correlation coefficients > 0.70 to a given set of initial conditions for a 5M ⊕ super-Venus. Correlation coefficients quantifying the correspondence between the actual and predicted output parameters were calculated using normalized and mean subtracted input and output parameters.
crust. From Eq. 12, volumetric melt production for a planet may scale as where the subscript ⊕ denotes values for an Earth-mass planet and A m stands for the mantle surface area.
We can approximate δ using the representative interior models of Valencia et al. (2006), for which R ∝ M 0.262 , ρ m ∝ M 0.196 , and g ∝ M 0.503 . First, consider the thickness of a melting region, d m = z i -z f . Since z f = P f /(ρ L g) is approximately constant for any planet, where a roughly constant mantle to core thickness ratio is assumed, although planetary mantles grow slightly more than cores with increasing planetary mass. Next, The Rayleigh number for a massive planet scales as (26) Assuming that the first and second terms on the right hand side are roughly equal to unity, we have Because h c is usually much smaller than R P , (29) Finally, the rest of the scaling relations may simply be assumed as and Hence, δ ≈ 0.240 and (f m /f m, , where ∆t is the duration of crust growth, an increase in melt productivity with mass does not guarantee an increase in crustal thickness with mass. As planetary mass, and thus radius, increases, a larger volumetric melt production is required to produce a certain crustal thickness. Specifically, crustal thickness would only increase with mass for δ > 0.524, assuming that ∆t is roughly con-stant. Therefore, although melt productivity increases with planetary mass, this simple scaling analysis indicates that crustal thickness should decrease with scaling Panel (a) of Fig. 10 is a plot of model output presentday crustal thickness as a function of planetary mass for simulations of Mars, Venus, and seven super-Venus planets. While initial conditions strongly affect simulation results, the model outputs generally follow this simple scaling. Smaller planets can have thicker crust though they tend to be characterized by lower mantle temperatures.

Mantle Processing
The scaling of mantle processing with planetary mass follows easily from the above analysis. A simplified equation for the volume of processed mantle is where ∆t is a duration for crustal growth. Thus, the amount of processed mantle scales with planetary mass as so ξ ≈ δ ≈ 0.240. The volume of a super-Venus planet scales as so ζ = 0.786. Therefore, (V proc /V ) ∝ (M/M ⊕ ) −0.546 . Although the amount of processed mantle material increases with planetary mass, the fraction of processed mantle decreases with increasing planetary mass because the mantle volume increases more rapidly than the amount of processed material. Panel (b) in Fig. 10 confirms that the fraction of processed mantle does indeed decrease with increasing planetary mass according to this scaling law, although initial conditions strongly affect the simulation results.

Viscosity Contrasts During Stagnant-Lid Convection
The viscosity contrast across the lithosphere is tracked during each thermal evolution simulation, along with the critical viscosity contrast above which a planet is locked in the stagnant-lid regime. Figure 11 shows the output of 595 simulations for Mars, Venus, and two super-Venus planets for which Q 0 , T u (0), and µ were varied over a wide range. In particular, all permutations of Q 0 = 0.5, 1.0, and 1.75 ×10 −7 W m −3 (scaled as usual with ρ m ) and T u (0) = 1400, 1700, and 2000 K were considered for a range of µ between 0.0 and 0.9.
From this plot, several conclusions may be drawn. First, for values of the frictional coefficient associated with dry silicate rocks, µ ∼ 0.7 to 0.8, plate tectonics is never favored. Second, increasing planetary mass does not substantially affect the likelihood of plate tectonics. Third, the effects of choosing different initial conditions are amplified for greater planetary mass. Finally, although choosing extreme initial conditions can change the viscosity contrast by orders of magnitude, the effect of the friction coefficient is far more important.

Formation of an Eclogite Layer
At depth, crustal rock may undergo a phase transition to eclogite. To extend the simple analysis from earlier, we write the thickness of the crust in the eclogite stability field as h e = h c -d e , where d e is the depth of the phase boundary. Likewise, we consider δ e and δ a , the fractions of the crust in and above, respectively, the eclogite stability field. Because δ e +δ a = δ e,⊕ +δ a,⊕ = 1, we may write δ e = δ e,⊕ + δ a,⊕ 1 − δ a δ a,⊕ . (f) Figure 9: Summary of parameter values at the present for 675 simulations of the thermal evolution of a 10M ⊕ super-Venus. Arrows are projections of the principal component basis vectors that emanate from a point representing the averaged simulation results, indicating axes that account for the vast majority of the data set's variance. The red arrow represents a larger percentage of cumulative variance (42%) than the green arrow (30%). Panels show (a) Moho temperature, (b) mantle potential temperature, (c) surface heat flux, (d) mantle heat flux, (e) fraction of mantle processed by melting, and (f) total time for crust to grow from 10% to 95% of its present thickness as functions of crustal thickness. Because crustal melting causes highly discontinuous surface and mantle heat fluxes, the model output is the averaged values for the final 100 Myr of planetary evolution.
The fraction of the crust above the eclogite stability field may scale as (36) If we assume that pressure increases hydrostatically with depth and that the critical pressure below which the phase transition occurs is a constant, then d e ∼ 1/g. So, = 0.284 -0.503 = -0.219. Therefore, the fraction of crust in the eclogite stability field should increase with planetary mass as In thermal evolution models, the heat conduction equation is numerically solved to calculate crustal temperatures. An approximate temperature profile can also be calculated using a steady-state approximation as (Turcotte and Schubert, 2002) T where Q c , the volumetric crustal heat production, is calculated as where t p is 4.5 Gyr and V c is the volume of the crust. The boundary condition T (h c ) = T c is used to calculate surface heat flux for specified Moho and surface temperatures and magnitude of internal heat production. Finally, Eq. 38 is used to calculate the temperature profile throughout the entire thickness of the crust. Representative temperature profiles for a planet can be used to approximate the fraction of crust that lies within the eclogite stability field. Figure 12 shows representative temperature profiles for Venus and 5 and 10M ⊕ super-Venus planets, calculated using representative crustal thicknesses, degrees of mantle processing, and Moho temperatures from the previous sensitivity analyses. A range of internal radiogenic heating was also considered. The stability field of eclogite is taken from Philpotts and Ague (2009) and is drawn assuming a hydrostatic pressure increase with depth. Panel (d) in Fig. 12 summarizes the effects of initial conditions on the fraction of crust in the eclogite stability region after 4.5 Gyr and shows the scaling from Eq. 37. For Mars, Venus, and seven super-Venus planets, 54 thermal evolution simulations were run to study all permutations of the initial conditions T u (0) = 1700 and 2000 K and Venus-equivalent Q 0 = 0.5, 1.0, and 1.75 ×10 −7 W m −3 . As predicted, the fraction of eclogite crust increases with planetary mass.  (c) 10M ⊕ super-Venus planets, calculated assuming representative crustal thicknesses, degrees of mantle processing, and Moho temperatures. The green shaded area is the approximate stability field for eclogite, drawn using the phase diagram from Philpotts and Ague (2009) and assuming a hydrostatic pressure increase with depth. Black solid and red dashed lines represent Venus-equivalent Q 0 = 1.0 × 10 −7 and 1.75 × 10 −7 W m −3 , respectively. Panel (d) shows the fraction of crust in the eclogite phase for Mars, Venus, and seven super-Venus planets. Circles and triangles represent Tu(0) = 1700 and 2000 K, respectively. Blue, black, and red symbols represent Venus-equivalent Q 0 = 5.0 × 10 −8 , 1.0 × 10 −7 , and 1.75 × 10 −7 W m −3 , respectively.

Pressure Effects on Mantle Rheology
The rheological behavior of the mantles of large rocky planets is difficult to predict.
While the core/mantle boundary pressure for Earth is ∼135 GPa (e.g., Dziewonski and Anderson, 1981), pressures within the silicate mantles of large rocky planets likely exceed 1 TPa (e.g., Valencia et al., 2006). Above the transition to post-perovskite at ∼120 GPa, Earth's mantle is primarily made of MgSiO 3 -perovskite and MgO (e.g., Murakami et al., 2004). In contrast, much of the mantles of super-Venus planets will be dominated by post-perovskite and perhaps, above 1 TPa, a mixture of MgO and SiO 2 (Umemoto et al., 2006). Unfortunately, we lack experimental measurements of the properties of planetary materials under these extreme conditions. Until such data are available, conjectures about the rheology of the silicate mantles of large rocky planets will remain controversial. Thermal evolution simulations are very sensitive to assumed rheological behaviors, so investigating the implications of various possible assumptions is essential.
The viscosities of most planetary materials increase with pressure when examined at relatively low pressures (e.g., Karato, 2008). Simple extrapolation of this trend predicts extreme increases in viscosity within massive terrestrial planets. Extensions of known perovskite rheology, for instance, imply an increase of >15 orders of magnitude as pressure increases to 1 TPa in an adiabatic mantle (Stamenkovic et al., 2011). Specifically, a viscosity profile may be calculated as (Stamenkovic et al., 2012) where η 0 is a reference viscosity at the reference temperature T * = 1600 K and V * is an activation volume. Figure 13 shows calculated viscosity profiles for η 0 = 10 21 Pa s and V * = 2.5, 1.7, and 0.0 cm 3 mol −1 for the convecting, adiabatic mantle within a 10M ⊕ super-Venus planet, following Stamenkovic et al. (2012). Our parameterized formulation for stagnant-lid convection is based on numerical modeling with the incompressible fluid approximation using temperature-dependent but pressure-independent viscosity (e.g., Solomatov and Moresi, 2000;Korenaga, 2009). Therefore, we are assuming pressure-independent constant potential viscosity, where the effect of temperature increase along an adiabatic gradient exactly balances the effect of pressure on viscosity, which requires a non-zero, positive activation volume. With the linear temperature gradient in Fig. 13, the assumption of constant potential viscosity corresponds to V * = 0.22 cm 3 mol −1 .
Strongly pressure-dependent viscosity can cause dramatically different behavior to emerge from parameterized convection models, including sluggish lower mantle convection and even the formation of a conductive lid above the core/mantle boundary (Stamenkovic et al., 2012). If convection were effectively suppressed in the The black curve shows internal temperature as a function of depth in the convecting mantle. Green, red, and blue dashed lines represent viscosity profiles calculated using Eq. 40 for V * = 0, 0.22, 1.7, and 2.5 cm 3 mol −1 , respectively.
lower mantle of large rocky planets, melt production would be significantly decreased and the likelihood of plate tectonics might decrease along with convective vigor as planetary mass increased. The thermal conductivity and expansivity of the mantle are also predicted to increase and decrease, respectively, with depth because of increasing pressure, but the effects of these changes on mantle dynamics are dwarfed by the putative increase in mantle viscosity (Stamenkovic et al., 2011(Stamenkovic et al., , 2012. Although an increase in viscosity under greater pressures seems intuitive, the straightforward application of limited, low-pressure experimental data may not accurately describe the rheology of massive terrestrial planets. In fact, four mechanisms, including a transition from vacancy to interstitial diffusion mechanisms, may cause a viscosity decrease with depth above a pressure of ∼0.1 TPa (Karato, 2011), as post-perovskite and additional high-pressure mineral phases dominate mantle rheology. According to this study, viscosities in the deep interiors of super-Venus planets may be less than the viscosity of Earth's lower mantle, potentially by as much as 2-3 orders of magnitude. Moreover, the depth-dependence of viscosity within Earth's perovskite-dominated mantle is still debated because Earth's viscosity profile has not been well-constrained.
Despite basic consensus that Earth's lower mantle is probably more viscous than the upper mantle, the magnitude of the viscosity contrast remains controversial. Early studies of Earth's topography and geoid suggested a ∼300-fold increase in viscosity between Earth's upper and lower mantle (Hager and Richards, 1989). On the other hand, analyses of post-glacial rebound predict an order of magnitude less of viscosity increase (e.g., Kaufmann and Lambeck, 2002) and gravity data are consistent, albeit loosely, with uniform or only slightly depth-dependent mantle viscosity (Soldati et al., 2009). A joint inversion of these data sets predicts that Earth's internal viscosity increases by ∼2 orders of magnitude throughout the mantle (Mitrovica and Forte, 2004), but it is noted that such inversions are known to suffer from severe nonuniqueness (e.g., King, 1995;Kido and Cadek, 1997). In any case, considering the tremendous uncertainty surrounding viscosity profiles within Earth and putative super-Venus planets, it remains legitimate to investigate how large rocky planets might evolve in the limiting case of pressure-independent potential viscosity.

Escaping the Stagnant-Lid Regime
Terrestrial planet evolution strongly depends on the regime of mantle convection. Assuming that brittle failure limits the strength of the lithosphere, our simulations indicate that the effects of lithosphere hydration dominate the effects of planetary mass on yield and convective stresses. That is, the increase in convective vigor with planetary mass only makes plate tectonics marginally more likely. Modeling results for super-Venus planets, however, suggest two additional mechanisms for escaping the stagnant-lid regime. First, massive terrestrial planets in the stagnant-lid regime feature crustal temperature profiles that enter the stability field of eclogite after crust grows beyond a critical thickness. If a sufficiently large fraction of the total crustal thickness is composed of eclogite, the entire crust could be gravitationally unstable and susceptible to foundering because eclogite is intrinsically denser than mantle peridotite.
On Earth, the phase transition from (metamorphosed) basalt to eclogite primarily occurs in subduction zones. Hydration may thus be important to allowing this phase transition to occur relatively rapidly (e.g., Ahrens and Schubert, 1975), although this type of kinetic calculation strongly depends on diffusion data that are not well-constrained (e.g., Namiki and Solomon, 1993). Eclogite is also formed during continent-continent collisions such as the Eurasian and Indian plate collisions (Bucher and Frey, 2002). Furthermore, the high density of eclogite is theorized to have caused delamitation, foundering, and recycling of relatively thick oceanic lithosphere on Earth during the Archaean (Vlaar et al., 1994). Finally, eclogite may be produced in large mountain ranges through magmatic differentiation (Ducea, 2002) and pressureinduced phase transition (Sobolev and Babeyko, 2005) in thick continental crust. Evidence for the strong influence of recent eclogite production and foundering on the topography of the central Andes Mountains has been gathered through geodynamics, petrology, and seismology (e.g., Kay and Abbruzzi, 1996;Beck and Zandt, 2002;Sobolev and Babeyko, 2005;Schurr et al., 2006), as synthesized in a numerical study (Pelletier et al., 2010). Because massive terrestrial planets have relatively high surface gravity, the phase transition to eclogite will occur at a comparatively shallow depth, making eclogite the stable mineral phase for a large fraction of the crust. The formation of a thick eclogite layer then could cause lithosphere foundering or intermittent plate tectonics, as has been proposed in episodic subduction mechanisms for Venus (Turcotte, 1993;Fowler and O'Brien, 1996).
Representative temperature profiles for massive terrestrial planets pass through the eclogite stability field for plausible initial conditions. In fact, radiogenic heating and thus crustal temperatures should be greater than calculated with Eq. 38, which would increase the speed of the phase transition to eclogite, because crustal construction mostly occurs early in planetary history. For super-Venus planets with masses greater than ∼4M ⊕ , eclogite may be the stable phase for the majority of the crust unless the initial mantle potential temperature or magnitude of internal heating is very low. So, crust material may undergo a phase transition to eclogite at relatively shallow depths as the crust grows during thermal evolution in the stagnant-lid regime, forming a thick eclogite layer that could subsequently founder. The buoyant stress from the presence of eclogite scales as ∆ρgh e , where ∆ρ ∼ 100 kg m −3 is the difference between the densities of eclogite and the mantle. On the other hand, the lithospheric strength scales as µρ L gd e , which can only be overcome when the depth scale of the eclogitic layer becomes large enough (at least locally, for example, by foundering). As long as crustal production continues on large rocky planets, eclogite formation and foundering could occur periodically, possibly yielding a regime of mantle convection resembling intermittent plate tectonics. Although we suggest that this process is plausible, pursuing its dynamics in detail is left for future studies.
High surface and crustal temperatures may also cause periodic transitions from the stagnant-lid regime to a form of mobile-lid convection. In this work, massive terrestrial planets in the stagnant-lid regime with surface temperatures held constant at 300 K tend to have very hot crusts. If high surface temperatures exist alongside high crustal temperatures, a transition from stagnantlid convection to a mobile-lid regime can occur (Reese et al., 1999). Feedback between a changing mantle convection regime and a periodic atmospheric greenhouse effect driven by varying amounts of volcanism, for instance, may be very important to the evolution of Venus (Noack et al., 2011). As surface temperature depends on atmospheric mass and the composition and luminosity of the central star, however, this possibility of escaping the stagnant-lid regime may not be as robust as the first mechanism based on the formation of selfdestabilizing crust

Limitations of Parameterized Models
Any parameterized model suffers shortcomings. Steady-state evolution is assumed, for instance, which poorly captures transient events that occur early in planetary evolution such as large impacts (Agnor et al., 1999) and the crystallization of a magma ocean (e.g., Solomatov and Stevenson, 1993). Fundamental assumptions such an adiabatic temperature gradient in the mantle and pressure-independent potential viscosity are controversial, and different approaches such as mixing length theory (Wagner et al., 2011) may be necessary to calculate the thermal structure of planetary interiors if they are not valid. However, our simplified simulations only aim to illuminate the first-order, relative effects of planetary mass on terrestrial planet evolution. Recreating the thermal history of a particular planet would require the introduction of many additional complications. One-dimensional models only return globally averaged values for important quantities, for instance, but mantle plumes, which may upwell from the core/mantle boundary when the core heat flux is positive, are likely important to local magmatism and surface features on terrestrial planets like Mars (e.g., Weizman et al., 2001) and Venus (e.g., Smrekar and Sotin, 2012). Furthermore, applying a parameterized approach to Venus, where the precise quantity of magmatism is an key output, requires more computationally intensive simulations to benchmark the relevant scaling laws.

Conclusions
Terrestrial planet evolution is complicated. Although plate tectonics is observed on Earth, the stagnant-lid regime of mantle convection may be most natural for terrestrial planets; at least, it is most common in our Solar System. Thermal evolution models in this study yield first-order, relative conclusions about the evolution of generic terrestrial planets in the stagnant-lid regime. Principal component analysis of simulation results conducted with a wide range of initial conditions captures the relationships between the large number of parameters that describe the interior of a planet. Depending on initial conditions, these planets may have evolved along a variety of paths, featuring different crustal thicknesses and temperatures, interior temperatures, and degrees of mantle processing. To produce specific histories consistent with spacecraft data obtained from Mars and Venus, complications must be added to these simple models.
Properties of massive terrestrial exoplanets are poorly constrained, so questions about the effects of planetary mass on the likelihood of plate tectonics and other important planetary parameters await definitive answers. In this study, we explored what might happen if internal viscosity is not strongly-pressure dependent, the alternative to which has been explored previously using parameterized models. Although convective vigor increases with planetary mass, the likelihood of plate tectonics is only marginally improved. Simple scaling analyses indicate that mantle melt productivity should increase with planetary mass. Because the increase in mantle processing is slow, however, crustal thickness and the relative fraction of processed mantle actually decrease with increasing planetary mass, as thermal evolution simulations confirm. Surface gravity increases with planetary mass, so pressure in the crust of massive terrestrial planets increases relatively rapidly with depth. Plausible temperature profiles favor a phase transition to gravitationally unstable eclogite during normal crustal formation, whereas the basalt to eclogite transformation rarely occurs aside from subduction on Earth. Therefore, thick eclogite layers, along with mobile, hot crustal material, may be important to the evolution of massive terrestrial planets.