Blue Marble, Stagnant Lid: Could Dynamic Topography Avert a Waterworld?

Limited by computing power, early constructions of a scaling function for dynamic topography have used a constant viscosity for the convecting region (Parsons & Daly 1983; Kiefer & Hager 1992). Under this isoviscous paradigm, a single dimensionless parameter, the Rayleigh number, describes the convective vigor of the system:

$$\begin{eqnarray}&&\mathrm{Ra}=\displaystyle \frac{\alpha \rho g{\rm{\Delta }}{{Td}}^{3}}{\kappa \eta },\end{eqnarray} \tag{ 1 }$$

where α is the thermal expansivity of the material in K⁻¹, ρ is its density in kg m⁻³, g is the surface gravity in m s⁻², ΔT is the temperature contrast across the layer in K, d is the depth of the convecting region in m, κ is the thermal diffusivity in m² s⁻¹, and η is the dynamic viscosity in Pa s—in isoviscous convection these parameters are all constant. The Ra number can act as a useful independent scaling variable for many convection phenomena because the vast majority of temperature variations in a convecting cell occur in its boundary layers (McKenzie et al. 1974). Boundary layer theory justifies a power-law relationship between Ra and the thickness of the upper thermal boundary layer. Hence, these previous works on isoviscous dynamic topography supposed scaling relationships of the form h/(αΔTd) ∼ Raⁿ (the αΔTd term ensures that both sides of the proportionality are dimensionless and n is uniquely defined).

In rocky planets, however, η changes with temperature (Karato & Wu 1993); steep viscosity gradients across the mantle are a defining trait of natural stagnant lid convection in that the cold surface is too viscous to flow (Davaille & Jaupart 1993; Solomatov 1995). Scalings based on Equation (1) defined using constant viscosity will not necessarily provide an optimal fit to the topography of stagnant lid bodies (Sembroni et al. 2017; Bodur & Rey 2019). In identifying a convecting system whose viscosity decreases quickly with temperature, we need a second dimensionless parameter in addition to a Rayleigh number: the viscosity contrast across the layer, Δη = η₀/η₁, where η₀ is the viscosity at the top and η₁ the viscosity at the bottom. A nonuniform viscosity profile implies many possible thermal Rayleigh numbers. Here Ra₁ denotes Equation (1) evaluated at η = η₁. In simple numerical models, viscosity is often assumed to follow an exponential law, $\eta (T)={\eta }_{0}\,\exp (-{bT})$, where the temperature prefactor $b=\mathrm{ln}\left({\rm{\Delta }}\eta \right)$ is a constant (Solomatov 1995).

Further, any Ra scaling function will only apply to its intended convection regime. Canonical studies of temperature-dependent viscosity convection distinguish between at least two series of regimes. These regimes have their own transitions in Ra₁-Δη space, which would manifest as discontinuities in the scaling function. A first series concerns the mobility of the surface: as Δη increases, a convecting system will move from a small viscosity contrast regime (similar to the isoviscous case) to a stagnant lid regime, via an intermediate regime of a sluggish lid (Moresi & Solomatov 1995; Solomatov 1995; Kameyama & Ogawa 2000). In a second series of transitions, the so-called stationarity of convection changes. As Ra₁ increases, the system will move roughly from a steady-state regime to a chaotic time-dependent regime, again through a transitional regime (Dumoulin et al. 1999; Solomatov & Moresi 2000). For either series, the regime boundaries are not sharp in Ra₁-Δη space but depend on this parameter space in a complex way via the aspect ratio of convection and the initial conditions. While ascribing Ra₁ presumes a bottom-heated convection cell, different modes of heating may also affect dynamic topography scaling relationships in ways we do not yet understand.

A waypoint objective of this work is therefore to develop a preliminary dynamic topography scaling relationship for the stagnant lid regime. While the topography of stagnant lid bodies has indeed been modeled numerically before (Moresi & Parsons 1995; Solomatov & Moresi 1996; Vezolainen et al. 2003, 2004; Orth & Solomatov 2011; Golle et al. 2012; Huang et al. 2013), the majority of this literature is directed at producing geoid-to-topography ratios to invert for interior properties of Venus or Mars, as opposed to fully exploring parameter space with forward models. As such, we are aware of no published scalings as explicit functions of the relevant convective parameters. Given the scope of our work here, we do not attempt to characterize the scaling behavior near the regime discontinuities (which would require a much finer grid of models in Ra₁-Δη space). Instead, we restrict ourselves to the chaotic time-dependent stagnant lid regime located in Moresi & Solomatov (1995) and Orth & Solomatov (2011) and simulated previously with Venus- and Mars-like parameters (Solomatov & Moresi 2000; Hauck 2002; Reese et al. 2005; Orth & Solomatov 2011). As such, we are assuming that chaotic stagnant lid convection will apply to most geodynamically active rocky planets—an assumption that may be tested in the future when detailed characterization of rocky exoplanets becomes possible.

In the second part of this study, we convert scaled dynamic topographies into the corresponding volumes of the largest possible ocean basin. The key product here is a spherical harmonic expansion of this scaled topography onto a Cartesian grid, as a synthetic elevation map. Yet all that our stagnant lid convection scaling law provides is a scalar height value. With some convenient assumptions about dynamic topography's spectral properties, it is in fact straightforward to find a power spectrum that is consistent with both the scalar height we have and some set of spherical harmonic coefficients we need.

Initial observations of Venus's, Earth's, and Mars's (total) topographies suggested a remarkably loglinear relationship between the 1D power spectral density, ${\phi }_{h}^{\mathrm{PSD}}$ in m³, and the wavenumber, k in m⁻¹. From spherical harmonic degree l = 5 to at least l = 100, the available spectra appeared consistent with a slope dlog${\phi }_{h}^{\mathrm{PSD}}$/dlog k ∼ −2 (Turcotte 1987; Rapp 1989; Balmino 1993). This precise slope value was predicted earlier still by Vening Meinesz (1951) and appears to be physically motivated (Sayles & Thomas 1978; Lovejoy et al. 1995)—perhaps emerging from sediment transport laws (Pelletier 1997, 1999; Roberts et al. 2019), although we will not be considering topography's modification by erosion explicitly. Statistically, a slope of −2 corresponds to red noise, the noise associated with a random walk process.

The convenient consequence of a loglinear spectral model—with a predetermined slope—is that it would let us approximate the shape of any planetary surface given just one free parameter, i.e., the y-intercept of ${\phi }_{h}^{\mathrm{PSD}}(k)$. As for dynamic topography in particular, models and Earth observations have indicated a shallower spectral slope roughly consistent with pink noise ∝k⁻¹, up to l ∼ 30 (Hoggard et al. 2016, 2017; Davies et al. 2019). However, there is no evidence that this spectral structure should characterize dynamic topography under all tectonic regimes. Hence, we extract the spectral structure of our own numerically modeled stagnant lid topography profiles. We will see that our rudimentary analysis again produces constant dlog${\phi }_{h}^{\mathrm{PSD}}$/dlogk values, albeit ones more strongly negative than −2. Observations of real stagnant lid bodies in the solar system could then suggest an empirical modification of this purely dynamic spectral model.

Our methods are described in Section 2. The approach we take is outlined as follows: We begin by extracting scaling relationships for the rms amplitude of dynamic topography from 2D numerical mantle convection simulations with temperature-dependent viscosity (Section 3.1.1). Second, we embed these scaling relationships in a suite of 1D parameterized convection models, allowing us to explore the sense of change of rms dynamic topography across a wide model parameter space and planet age distribution (Section 3.2.2). For this parameter study we focus on the planet mass, age, radiogenic element abundance, and core mass fraction, all relevant to the cooling history and Rayleigh number of a planet. We focus on these four parameters because they may be amenable to being observationally constrained for exoplanets, at least in principle. Third, we synthesize 2D maps from the projected rms amplitudes to see how the maximum capacity of ocean basins, and hence the minimum elevation gain needed for dry land, might trade off with planet size (Section 3.3). We end with a discussion of the study's limitations and applicability (Section 4).

Numerical computations were performed using the ASPECT code version 2.2.0 (Kronbichler et al. 2012; Heister et al. 2017; Bangerth et al. 2020). For each case we systematically varied two key input parameters: Ra₁ and Δη. Although we originally explored Ra₁ varying from 1 × 10⁷ to 3 × 10⁸, we found that simulations below Ra₁ = 1 × 10⁸ were not in the chaotic time-dependent regime and showed characteristically different topography scaling behavior. Because the present study was not designed to precisely locate these transitions, we focused only on the chaotic time-dependent regime. Simulations above Ra₁ = 3 × 10⁸ were found to be computationally impractical.

Our ASPECT implementation results in dimensionless temperature and velocity fields, denoted by the prime symbol. These and their derivative quantities can be dimensionalized as, e.g.,

$$\begin{eqnarray}\begin{array}{rcl}T & = & {\rm{\Delta }}T\,{T}^{{\prime} }+{T}_{0},\\ u & = & \displaystyle \frac{\kappa }{d}\,{u}^{{\prime} },\\ (x,y) & = & d\,({x}^{{\prime} },{y}^{{\prime} }),\\ {\delta }_{\mathrm{rh}} & = & d\,{\delta }_{\mathrm{rh}}^{{\prime} },\\ {\delta }_{\mathrm{lid}} & = & d\,{\delta }_{\mathrm{lid}}^{{\prime} },\\ h & = & \alpha {\rm{\Delta }}{Td}\,{h}^{{\prime} },\end{array}\end{eqnarray} \tag{ 2 }$$

where ${T}^{{\prime} }$ is the dimensionless temperature, ${u}^{{\prime} }$ is the horizontal component of the dimensionless velocity, ${\delta }_{\mathrm{rh}}^{{\prime} }$ is the dimensionless thickness of the upper thermal boundary layer, ${\delta }_{\mathrm{lid}}^{{\prime} }$ is the dimensionless thickness of the stagnant lid, ${h}^{{\prime} }$ is the dimensionless height of topography, T₀ is the dimensional surface temperature, ΔT is the dimensional temperature difference from bottom to surface, and the other (dimensional) parameters are defined under Equation (1) above.

All simulations involve a 2D rectangular box with fixed top and bottom temperatures, ${T}_{0}^{{\prime} }=0$ and ${T}_{1}^{{\prime} }=1$, respectively, and no internal heating. Free-slip boundary conditions are ascribed to the top and bottom surfaces, while reflecting boundary conditions are ascribed to the sides. We use a wide box with a nondimensional depth ${Y}^{{\prime} }$ of unity and a nondimensional width of ${X}^{{\prime} }=8{Y}^{{\prime} }$ to minimize the influence of the side walls. We assume an incompressible, infinite-Prandtl-number fluid and use the Boussinesq approximation. Viscosity is Newtonian and varies with temperature according to an exponential rheology law, ${\eta }^{{\prime} }=\exp (-b\,{T}^{{\prime} })$, where $b=\mathrm{ln}({\rm{\Delta }}\eta )$. We use the coarsest mesh size still able to resolve the lower thermal boundary layer; this varies for different Ra₁. Table 1 lists the relevant details of the model setup.

Each experiment is deemed to have reached quasi-steady-state when both its rms velocity stabilizes to within 0.002% and its top and bottom heat fluxes converge. All time steps prior to this point are discarded, and the models are then allowed to run for long enough such that the distribution of rms dynamic topography is well characterized. All cases are confirmed to be in the stagnant lid mode of convection based on the surface mobility criterion, $S={({\delta }_{0}^{{\prime} })}^{2}{u}_{0}^{{\prime} }\ll 1$, where ${\delta }_{0}^{{\prime} }={\delta }_{\mathrm{lid}}^{{\prime} }+{\delta }_{\mathrm{rh}}^{{\prime} }$ is the dimensionless thickness of the lithosphere and ${u}_{0}^{{\prime} }$ is the dimensionless surface velocity (Solomatov & Moresi 1997).

2.1.1. Extraction of Parameters From the Temperature and Velocity Profiles

2.1.2. Fitting a Topography Scaling Relationship

The ASPECT code calculates the horizontal profile of the surface dynamic topography via a stress balance at the center of each cell on the top boundary,

$$\begin{eqnarray}&&{\sigma }_{{yy}}=-\rho {gh},\end{eqnarray} \tag{ 3 }$$

where σ_yy is the vertical component of the stress imparted by convection, g is the gravity, and ρ is the mantle density. Equation (3) assumes mechanical equilibrium between the surface topography and the interior density structure, a safe assumption for the long timescales of convection (e.g., Ricard 2015). At each time step, we first normalize the dimensionless topography profile to ensure that its mean is zero, and then we find its rms value, ${h}_{\mathrm{rms}}^{{\prime} }$.

We choose the rms amplitude of topography as the representative scalar quantity to fit, rather than the peak amplitude. This choice is based on the reasoning that the rms value may be less sensitive to the model geometry—crucial for inferring 3D behavior from 2D models, as we will be doing. As such, we ran preliminary isoviscous convection simulations to confirm that neither Cartesian nor cylindrical 2D geometries show the same peak topographies as the equivalent 3D spherical experiments from Lees et al. (2020), whereas, for all three setups, the rms topographies align well. That the same result holds for non-isoviscous simulations is an outstanding caveat of this study.

Earlier in Section 1.1, we motivated the need for two parameters, a Rayleigh number and viscosity contrast, to fully describe stagnant lid convection. These will serve as the independent variables in the scaling function. We define an interior Rayleigh number,

$$\begin{eqnarray}&&{\mathrm{Ra}}_{i}=\frac{\alpha \rho g{\rm{\Delta }}{{Td}}^{3}}{\kappa \eta ({T}_{i})}={\mathrm{Ra}}_{1}\frac{\eta ({T}_{1})}{\eta ({T}_{i})},\end{eqnarray} \tag{ 4 }$$

that is, evaluating Equation (1) using the "interior" viscosity at ${T}_{i}^{{\prime} }$ (Solomatov & Moresi 2000). This formulation of the Rayleigh number is easily transferable to 1D convection models that predict a single mantle temperature and sidesteps any problems with predicting lower mantle viscosities (where pressure effects are important). Also with an eye toward 1D model integration, we use the exponential temperature prefactor $b=\mathrm{ln}({\rm{\Delta }}\eta )$ as the second variable. We anticipate a power-law relationship and thus fit a linear model to b, $\mathrm{log}({\mathrm{Ra}}_{i}$), and $\mathrm{log}({h}_{\mathrm{rms}}^{{\prime} }$), with an interaction term between b and $\mathrm{log}({\mathrm{Ra}}_{i}$):

$$\begin{eqnarray}&&\mathrm{log}{h}_{\mathrm{rms}}^{{\prime} }=A+{Bb}+C{\mathrm{logRa}}_{i}+D\left(b{\mathrm{logRa}}_{i}\right),\end{eqnarray} \tag{ 5 }$$

where ${T}_{i}^{{\prime} }$ in Equation (4) is determined from the horizontally and time-averaged temperature profile as per Section 2.1.1, ${h}_{\mathrm{rms}}^{{\prime} }$ is taken as the mean of the rms value over all time steps, and the log notation refers to the base-10 logarithm here and throughout. Thus, each experiment provides one (b, Ra_i , ${h}_{\mathrm{rms}}^{{\prime} }$) coordinate. While these data have some distribution due to the chaotic time dependence of convection, we found that including the standard error of the mean of $\mathrm{log}{h}_{\mathrm{rms}}^{{\prime} }$ has negligible effect on the regression results (for simplicity we do not consider the uncertainty on Ra_i ).

Coefficients A, B, C, D, and their covariance matrix are estimated using orthogonal distance regression. The interaction term, $D\left(b\mathrm{log}\,{\mathrm{Ra}}_{i}\right)$, accounts for cross-effects between b and Ra_i . Although including the interaction term adds an extra parameter, we will see that we need this term to properly capture the observed effect of Ra_i on ${h}_{\mathrm{rms}}^{{\prime} }$, which has magnitude and direction depending strongly on b as the data will show; the presence of the fourth term decreases the residual variance of the fit threefold compared to its absence.

In a fraction of the CPU time of a full dynamical convection simulation, parameterized convection models can result in similar temperatures to numerical models by tracking heat fluxes across the two thermal boundary layers (Thiriet et al. 2019). Parameterized convection can also produce a thermal history of the planet, from which we can extract a self-consistent evolution of dynamic topography. Further, such low-cost models invite parameter studies, which naturally we conduct in this segment. Important caveats are discussed in Section 4.

We will be exploring how topography changes with planet age, τ, mass, M_p , core mass fraction, CMF, and radiogenic heating expressed as an abundance of U and Th relative to the Sun, χ_rad. As such, these four parameters are independently and systematically varied between experiments. Meanwhile, we anticipate that some of the biggest uncertainties lie in the unknown mantle rheology. To see how these uncertainties would propagate, rather than testing their effect on h_rms explicitly, we will treat the parameters in the viscosity law as uniform random variables. In addition to the viscosity parameters, we also account for model uncertainty by drawing the topography scaling coefficients in Equation (5) from a multivariate normal distribution whose mean and covariance are given by the results of the regression from Section 2.1.2. Table 2 lists all dimensional input parameters used in the 1D model, which the remainder of this section describes.

Table 2. Dimensional Parameters Used in the 1D Thermal History Model

Symbol	Description	Value	Units	Reference
Constant Bulk Properties for All Planets
ρm	Mantle density	3500	kg m−3	Thiriet et al. (2019)
cm	Mantle specific heat	1142	J kg−1 K−1	Thiriet et al. (2019)
cc	Core specific heat	840	J kg−1 K−1	Thiriet et al. (2019)
km	Mantle thermal conductivity	4	W m−1 K−1	Thiriet et al. (2019)
αm	Mantle thermal expansivity	2.5 × 10−5	K−1	Thiriet et al. (2019)
κm	Mantle thermal diffusivity	1 × 10−6	m2 s−1	Thiriet et al. (2019)
Ra${}_{\mathrm{crit}}^{u}$	Critical Rayleigh number	450	⋯	Thiriet et al. (2019)
arh	Viscosity temperature scale coefficient	2.44	⋯	Thiriet et al. (2019)
β	Heat flow scaling exponent	1/3	⋯	Solomatov (1995)
Ts	Surface temperature	273	K
Variables Tested in the Parameter Study
τ	Planet age	2–4.5, baseline: 4.5	Gyr
Mp	Planet mass	0.1–5.0, baseline: 1.0	M⊕	Rogers (2015); Zeng et al. (2016)
CMF	Core mass fraction	0–0.4, baseline: 0.3	⋯	Zeng et al. (2016)
χrad	U and Th budget relative to solar	0.3–3.0, baseline: 1.0	⋯	Nimmo et al. (2020)
Unknown Random Variables
Ea	Viscosity activation energy	${ \mathcal U }(200,300)$	kJ mol−1	Karato & Wu (1993); Zhang et al. (2017)
η0	Viscosity prefactor	${ \mathcal U }(2.6\times {10}^{10},5.3\times {10}^{13})$	Pa s	See Section 2.2.3 in text
A, B, C, D	Topography scaling coefficients	${ \mathcal N }({\boldsymbol{\mu }},{\boldsymbol{\Sigma }})$ a	⋯	This work

Notes. The top section lists parameters that are constant in all runs. The middle section lists those parameters that are systematically varied in certain sections of the study and held constant at the baseline value where noted. The bottom section lists the unknowns, treated here as random variables distributed as given, such that a distribution of output parameters is obtained.

^a With mean μ and covariance Σ given by the results of the linear regression (see Section 2.1.2 and Table 3).

Download table as:  ASCIITypeset image

2.2.1. Governing Energy Balances

The approach outlined here closely follows that of Thiriet et al. (2019). The mantle and core temperatures are governed by the 1D energy balances,

$$\begin{eqnarray}\begin{array}{rcl}{M}_{m}\,{c}_{m}\displaystyle \frac{{\rm{d}}{T}_{m}}{{\rm{d}}t} & = & -{q}_{{\rm{u}}}\,{A}_{{\rm{u}}}+{q}_{\mathrm{rad}}\,{M}_{m}+{q}_{c}\,{A}_{c},\\ {M}_{c}\,{c}_{c}\displaystyle \frac{{\rm{d}}{T}_{c}}{{\rm{d}}t} & = & -{q}_{c}\,{A}_{c},\end{array}\end{eqnarray} \tag{ 6 }$$

where t is time in s, M_m is the mass of the convecting part of the mantle in kg, c_m is the mantle specific heat capacity in J kg⁻¹ K⁻¹, q_rad is the radiogenic heat flux in W kg⁻¹, q_u is the heat flux out of the top of the convecting region in W m⁻², and A_u is the surface area of the top of the convecting region in m². The subscript u denotes the upper boundary layer; the analogous notation with subscript c applies to the core. M_c is found through the core mass fraction. Just as in the 2D models, we explicitly include a mechanical stagnant lid, sitting atop the upper thermal boundary layer, never participating in convection. ³ Our choice of initial conditions for the governing equations is explained in Section 2.2.5.

Note also that we assume a perfectly spherical planet. For simplicity, and for consistency with our assumption of incompressibility in the 2D models, we treat c_m and other thermodynamic quantities as constant throughout the mantle (i.e., always equal to their reference values at the top of the convecting mantle); in reality these would vary with the adiabatic profile. This assumption would be a greater source of error for more massive planets with higher pressures at the base of the lithosphere. Although Equation (6) simplifies the problem by omitting other heat fluxes like volcanism (see Section 4.4.3), it will suffice in capturing the essential behavior of a cooling convective planet (Jaupart et al. 2015).

2.2.2. Interior Structure

The radius of the planet, R_p , is based on the physically motivated mass–radius relation in Zeng et al. (2016),

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{p}}{{R}_{\oplus }}=(1.07-0.21\,\mathrm{CMF}){\left(\displaystyle \frac{{M}_{p}}{{M}_{\oplus }}\right)}^{1/3.7},\end{eqnarray} \tag{ 7 }$$

while the radius of the core, R_c , is from Zeng & Jacobsen (2017),

$$\begin{eqnarray}&&{R}_{c}={R}_{p}\,{\mathrm{CMF}}^{0.5}.\end{eqnarray} \tag{ 8 }$$

We use a surface gravity g_s consistent with M_p and R_p . Note that Table 2 suggests that the mantle density, ρ_m , is a constant, but Equations (7) and (8) assume that density decreases radially outward such that gravity is constant through the mantle. Our box model can be said to treat ρ_m as a near-surface value, apt for the upper thermal boundary layer typically found at r ≈ 0.99R_p . Note that Equations (7) and (8) entail extrapolating equations of state to pressures beyond their validity range, which could lead to errors in R_p and R_c , compared to more accurate high-pressure equations of state such as in Hakim et al. (2018). Even at 5 M_⊕, however, the radius predicted by Equation (7) is 1.2% smaller than that from Hakim et al. (2018) for an Earth-like core size. This radius error has no effect on rms dynamic topography but decreases ocean basin sizes by 8%. Significant errors in dynamic topography predictions would come with R_p overinflations of more than 20%. In detail, accurate mass–radius relations will require tailoring to specific bulk compositions.

In the parameter study, we vary CMF from 0.0 to 0.4, the quoted range for which Equation (7) is valid. Neglecting any potential silicate mass loss after planet differentiation, oxidation chemistry predicts a theoretical upper CMF of 0.34 (Dyck et al. 2021). We consider values of M_p ranging from 0.1 to 5 M_⊕, corresponding to a Mars-sized body and to an equivalent radius slightly below the accustomed upper limit for rocky planets at 1.6 R_⊕ (Rogers 2015) based on Equation (7) with a CMF of 0.33.

2.2.3. Mantle Rheology

The rheology of rocky mantles is thought to obey an Arrhenius law (Karato & Wu 1993). The Arrhenius functional form yields exceedingly large viscosity contrasts over the cold lithosphere—spawning numerical issues in 2D models that preclude its use there. We exploit the Arrhenius form in the 1D model, but to maintain consistency between our 1D and 2D models, we ignore any pressure dependence and non-Newtonian behavior. We adopt a canonical law for diffusion creep as a function of temperature,

$$\begin{eqnarray}&&\eta (T)={\eta }_{0}\exp \left(\displaystyle \frac{{E}_{a}}{{R}_{b}T}\right),\end{eqnarray} \tag{ 9 }$$

where η is the dynamic viscosity in Pa s, R_b = 8.314 is the gas constant in J mol⁻¹ K⁻¹, E_a is the activation energy in J mol⁻¹, and η₀ is a prefactor with the same units as η. Note that our definition of η₀ does not act as a "reference viscosity" sometimes employed; it just encompasses all pre-exponential terms. In natural systems, the mantle viscosity will also depend on pressure; this caveat is discussed in Section 4.2.

In testing variations of η₀ and E_a , we shall try to capture the uncertainty imparted by unconstrained exoplanet rheologies. Strain rates brought on by the diffusion creep of silicate mantle rock would be strongly affected by both the water content and the bulk mineralogy. For olivine, Karato & Wu (1993) give the canonical wet (water-saturated) and dry (water-free) flow laws: E_a from 240 kJ mol⁻¹ in the former to 300 kJ mol⁻¹ in the latter; water weakens the rock. For the pre-exponential coefficient η₀, the same canonical laws correspond to 1.6 × 10¹¹ and 2.6 × 10¹¹ Pa s, which produces a dry olivine viscosity of ∼10²¹ Pa s at 1600 K.

We also expect to find overall higher viscosities inside planets that have mantles with lower Mg/Si compared to Earth's value of ∼1.3 (Pagano et al. 2015; Spaargaren et al. 2020; Ballmer & Noack 2021). At Mg/Si <1, the upper mantle composition would be dominated by orthopyroxene; at Mg/Si near 2 it would approach pure olivine. Our coarse treatment considers some empirical end members. We have laws for olivine; Zhang et al. (2017) give an Arrhenius flow law for the diffusion creep of enstatite. They find that E_a = 200 kJ mol⁻¹, that wet enstatite is approximately 10 times more viscous than wet olivine at depth, and that virtually dry enstatite is about 100 times more viscous than wet enstatite.

So far, this simple mineralogical paradigm would imply that water-saturated regions of Earth's upper mantle would exhibit the weakest-possible diffusion creep among rocky planets. To be conservative, we set a minimum η₀ of 2.6 × 10¹⁰ Pa s, an order of magnitude weaker than wet olivine (Karato & Wu 1993). The maximum η₀ is set at 5.3 × 10¹³ Pa s, approximating a dry enstatite rheology (Zhang et al. 2017). We test E_a between 200 and 300 kJ mol⁻¹. Both E_a and η₀ are drawn from random uniform distributions. By varying these parameters independently, we are likely overestimating the true uncertainty if they are in fact correlated. Note that we do not self-consistently adapt other bulk properties to account for the unknown mineralogy (an invaluable endeavor, but outside the scope of the current manuscript).

2.2.4. Heat Fluxes

Internal heating.—The radiogenic heat flux at t is

$$\begin{eqnarray}\begin{array}{rcl}{q}_{\mathrm{rad}} & = & \displaystyle \sum _{i=1}^{4}{\chi }_{i}{c}_{i}{h}_{i}\exp \left[\left(\tau -t\right)\frac{\mathrm{ln}2}{{\tau }_{1/2,i}}\right],\\ {\chi }_{i} & = & \left\{\begin{array}{ll}{\chi }_{\mathrm{rad}} & \mathrm{if}\ \ i\geqslant 2\\ 1 & \mathrm{otherwise}\end{array}\right.,\end{array}\end{eqnarray} \tag{ 10 }$$

where we are summing over the heat-producing isotopes ⁴⁰K, ²³⁸U, ²³⁵U, and ²³²Th; c_i is the present-day bulk silicate Earth concentration of the ith isotope in kg kg⁻¹; h_i is the heating contribution in W kg⁻¹; and τ_1/2,i is the half-life in the same units as t. Values for these parameters are taken from Table 1 in O'Neill et al. (2020). Further, for the refractory elements U and Th, we multiply the summand by a common factor χ_rad to reflect potentially extraterrestrial variations in the abundances of these r-process elements. As surveyed in Nimmo et al. (2020), U and Th abundances are conservatively expected to vary across Sun-like stars from between 30% to 300% of the solar value, which—assuming that relative mantle concentrations directly reflect relative stellar abundances (Thiabaud et al. 2015; Hinkel & Unterborn 2018; Putirka & Rarick 2019; Adibekyan et al. 2021)—translates to a range in q_rad of 2.22–14.34 pW kg⁻¹ at 4.5 Gyr, with the baseline value equivalent to 5.36 × 10⁻¹² W kg⁻¹. (We ignore the unconstrained variations in ⁴⁰K, a volatile isotope that in any case contributes less heating with age than refractory U and Th.) Although we do not account for the galactic chemical evolution of U and Th abundances as a function of stellar age (Frank et al. 2014), some of this variation is captured in χ_rad regardless.

Thermal boundary layers.—Across the upper and lower thermal boundary layers, heat fluxes are conductive:

$$\begin{eqnarray}&&{q}_{u,c}={k}_{m}\displaystyle \frac{{\rm{\Delta }}{T}^{u,c}}{{\delta }_{\mathrm{rh}}^{u,c}},\end{eqnarray} \tag{ 11 }$$

where k_m is the mantle thermal conductivity in W m⁻¹ K⁻¹, ΔT^u (respectively ΔT^c ) is the temperature contrast across the upper (lower) boundary layer in K, and ${\delta }_{\mathrm{rh}}^{u}$ (${\delta }_{\mathrm{rh}}^{c}$) is the thickness in m.

The thermal boundary layer thicknesses are controlled by their local Rayleigh numbers:

$$\begin{eqnarray}&&{\delta }_{\mathrm{rh}}^{u,c}=({R}_{\mathrm{lid}}-{R}_{c}){\left(\displaystyle \frac{{\mathrm{Ra}}_{\mathrm{crit}}^{u,c}}{{\mathrm{Ra}}_{\mathrm{rh}}^{u,c}}\right)}^{\beta },\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\mathrm{Ra}}_{\mathrm{rh}}^{u,c}=\displaystyle \frac{\alpha \rho {g}^{u,c}{\rm{\Delta }}{T}^{u,c}{({R}_{\mathrm{lid}}-{R}_{c})}^{3}}{\kappa \eta ({T}^{u,c})},\end{eqnarray} \tag{ 13 }$$

where Ra${}_{\mathrm{rh}}^{u,c}$ is the local Rayleigh number, Ra${}_{\mathrm{crit}}^{u}$ is the critical Rayleigh number for convection, and β is a constant that can be obtained from either experiments or theory. For both thermal boundary layers we take β = 1/3, such that q_u is independent of d; the boundary layers are assumed to be in a state of marginal stability (e.g., Solomatov 1995). The value of β is tied physically to the planet's dominant cooling mechanism, which strongly depends on the tectonic mode (Lenardic 2018; Seales & Lenardic 2020). The choice made here is appropriate for chaotically time-dependent, stagnant lid convection with temperature-dependent viscosity (Solomatov 1995; Solomatov & Moresi 2000). Other fitting choices do not significantly change our results (Thiriet et al. 2019).

For the upper thermal boundary layer, we have ΔT^u = T_m −T_lid, η(T^u ) = η(T_m ), and g^u = g_s , and we fix Ra${}_{\mathrm{crit}}^{u}$ at 450. For the lower layer, this becomes ΔT^c = T_c − T_m , η(T^c ) =η[(T_c + T_m )/2], g^c the gravity at R_c , and after Deschamps & Sotin (2000), Ra_crit,c = 0.28Ra${}_{i}^{0.21}$, with Ra_i the interior Rayleigh number defined for 1D convection in Equation (17). Although Ra_crit,c can be tricky to parameterize, T_c tends to equilibriate with T_m fairly quickly under this setup; hence, q_c ≪ q_u .

Finally, the temperature T_lid at the base of the lid in K (identically, at the top of the convecting region) is obtained for parameterized convection in a similar way to numerical models. The temperature drop between T_m and T_lid is proportional to the so-called viscous temperature scale, ΔT_ν (Davaille & Jaupart 1993):

$$\begin{eqnarray}&&{T}_{\mathrm{lid}}={T}_{m}-{\rm{\Delta }}{T}_{\mathrm{rh}}={T}_{m}-{a}_{\mathrm{rh}}{\rm{\Delta }}{T}_{v},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{T}_{\nu }=\displaystyle \frac{\eta ({T}_{m})}{{\rm{d}}\eta /{\rm{d}}T{| }_{{T}_{m}}}=\displaystyle \frac{{R}_{b}{T}_{m}^{2}}{{E}_{a}}.\end{eqnarray} \tag{ 15 }$$

The coefficient a_rh is empirically determined; we adopt a value of 2.44 for β = 1/3 based on the fits of Thiriet et al. (2019) to 3D spherical convection simulations. The radius R_lid of this temperature coordinate is described in the next section.

2.2.5. Stagnant Lid Thickness and the Final Governing Equation

The lid does not instantly grow or shrink in response to a change in the heat flux coming from the upper thermal boundary layer. Rather, there is a lag in which δ_lid adjusts such that the difference between the flux out of the top of the lid and the flux into the base of the lid is minimized:

$$\begin{eqnarray}&&\frac{d{\delta }_{\mathrm{lid}}}{dt}=\frac{{q}_{\mathrm{lid}}{| }_{{R}_{\mathrm{lid}}}-{q}_{u}}{{\rho }_{m}{c}_{m}({T}_{m}-{T}_{\mathrm{lid}})},\end{eqnarray} \tag{ 16 }$$

where the heat flux profile of the lid, q_lid(r) in W m⁻², is obtained by solving the steady-state conductive heat transfer equation in spherical geometry with boundary conditions (R_lid, T_lid) and (R_p , T_s ) where T_s , is the surface temperature in K, and with internal heating equal to the mantle q_rad (in reality, we might anticipate higher concentrations of lithophiles U, Th, and K in the lid). This steady-state formulation ignores the time dependence of heat conduction in the lid, leading to errors compared to a time-dependent model in the surface heat flux of ≲5 mW m⁻² for a Mars-sized planet. A smaller error is expected for larger planets with thinner lids (Thiriet et al. 2019).

We account for the mass of the convecting region changing with δ_lid by subtracting the lid mass, ${\rho }_{m}4\pi /3({R}_{p}^{3}-{R}_{\mathrm{lid}}^{3})$, from the fixed quantity M_p (1 − CMF). At each time step we also update R_lid = R_p − δ_lid. Thus, Equation (16) presents a third differential equation that must be solved simultaneously with Equation (6). We solve this system of equations using the explicit Runge–Kutta method of order 5. The initial conditions T_m,0, T_c,0, and δ_lid,0 reflect the unknown formation history of the planet—the leftover gravitational energy of accretion and core segregation, and the crystallization of the primordial magma ocean(s). To bypass this uncertainty, we only consider simulations that have converged to a memoryless state. That is, we prime each experiment by running it forward from t = −5 to 0 Gyr and using the solution at 0 Gyr as the initial conditions. Then, Equations (6) and (16) are solved again from t = 0 to τ.

2.2.6. Dynamic Topography

Once we have a solution for the planet's thermal history, we combine these results with Equation (5) to find ${h}_{\ \mathrm{rms}}^{{\prime} }$. Since we have b and Ra_i forming the basis of the topography scaling from 2D experiments, applying Equation (5) to 1D thermal histories requires writing 1D-appropriate analogs of these two variables. An analog of Ra_i is quite straightforward; for parameterized convection this variable is defined a posteriori as

$$\begin{eqnarray}&&{\mathrm{Ra}}_{i}=\displaystyle \frac{{\alpha }_{m}{\rho }_{m}{g}_{s}{\rm{\Delta }}T{\left({R}_{p}-{R}_{c}\right)}^{3}}{{\kappa }_{m}\eta ({T}_{m})},\end{eqnarray} \tag{ 17 }$$

where ΔT = T_s − T_c . This equation is the same as Equation (4) using the dimensional parameters for the mantle in Table 2 and simply letting the interior viscosity η(T_i ) = η(T_m ). For our runs, T_c ≈ T_m . Note also that Ra_i differs from Ra${}_{\mathrm{rh}}^{u}$ (Equation (13)) in that the latter excludes the stagnant lid from its domain.

Meanwhile, b as defined in the exponential viscosity law must be related to Arrhenius law parameters, since in 1D viscosity is modeled with the latter, more realistic law. Moresi & Solomatov (1995) demonstrate such an exponential approximation to an Arrhenius law. The approximation comes from the idea that in the stagnant lid regime it is the local rheological gradient over the upper thermal boundary layer that propels temperature-dependent viscosity convection, rather than the total domain viscosity contrast, Δη (Davaille & Jaupart 1993). One can therefore write $\eta (T)\sim \exp \left[\left({\rm{\Delta }}T/{\rm{\Delta }}{T}_{\nu }\right)T\right]$, where the viscous temperature scale ΔT_ν is rescaled by ΔT to make the temperature prefactor dimensionless. From Equation (15) this implies

$$\begin{eqnarray}&&b=\displaystyle \frac{{\rm{\Delta }}T}{{R}_{b}{T}_{m}^{2}/{E}_{a}}.\end{eqnarray} \tag{ 18 }$$

In 2D applications, setting T_m at the interior temperature just below the upper thermal boundary layer would create a viscosity profile that is most closely aligned to the Arrhenius profile, especially over the key region of the upper thermal boundary layer (Moresi & Solomatov 1995).

Finally, the dimensionless ${h}_{\mathrm{rms}}^{{\prime} }$ resulting from Equations (18), (17), and (5) is scaled by α_m ΔTd (Equation (2)) to get the dimensional h_rms. To clarify, we do consider the whole domain in the dimensionalization, so d = R_p − R_c and again ΔT = T_c −T_s ; the fact that several of these constituents evolve with time means that h_rms is a function of the age of the planet.

These calculations so far have assumed subaerial topography. Water-loaded topography would be higher by a factor of ρ_m /(ρ_m − ρ_w ) ≈ 1.5, where ρ_w is the density of water.

It is worth mentioning at this point that the dependence in several places on T_s —inside the definition of b in particular—means that there is a certain sensitivity of h_rms to this free parameter. For example, all else held constant at the baseline value (Table 2), increasing T_s from 273 to 373 K is associated with a 30% decrease in h_rms. However, because this study is only concerned with temperate planets that have a narrow range in T_s , we do not consider its effect on topography.

We have based our scaling relationship on h_rms (Section 2.1.2), yet it is the peak topography, h_peak, that controls how much water a planet's surface reservoirs can hold at the maximum capacity. Therefore we require the peak topography associated with an rms value in a 3D spherical geometry, given assumptions about topography's distribution.

Appendix A explains the relevant spherical harmonics method in more detail. Suppose we have a loglinear power spectrum, which fiducially describes dynamic topography amplitudes on a sphere. Essentially, for each run of 1D thermal evolution, we transpose the power spectrum vertically such that its frequency-domain rms value matches the spatial-domain rms value expected from the h_rms(Ra_i , b) scaling function. The transposed spectrum is expanded onto a 2D map, h(x, y), which has its own ${h}_{\mathrm{peak}}=\max (h)$. The volumetric ocean basin capacity in cubic meters—the main intended application of our topography modeling—is estimated as

$$\begin{eqnarray}&&\displaystyle \begin{array}{l}{V}_{\mathrm{cap}}=\frac{{\rho }_{m}}{{\rho }_{m}-{\rho }_{w}}\int \left[{h}_{\mathrm{peak}}-h(x,y)\right]\\ \,\,dS=\frac{{\rho }_{m}}{{\rho }_{m}-{\rho }_{w}}4\pi {R}_{p}^{2}{h}_{\mathrm{peak}},\end{array}\end{eqnarray} \tag{ 19 }$$

where the integral is over the surface S and the 2D map is multiplied by the density ratio term to account for water-loaded topography (our purpose here entails that the whole map is underwater, save for the single grid point corresponding to h_peak). The actual basin capacities of Venus, Earth, and Mars defined this way are 3.4, 3.3, and 2.9 Earth oceans, respectively—we expect to find lower values by considering only dynamic topography.

Robust models of dynamic topography power spectra are not available at this time. Instead, for the spectrum needed above, we explore three hypothetical scenarios. The first and simplest model is that all topography behaves like red noise, as per the historical paradigm introduced in Section 1.2 (e.g., Turcotte 1987). The second option is to represent empirical dynamic topography with the observed shape of Venus—although broad regions of Venus's highlands indicate isostatic support, so the resulting spectral distribution should reflect a mix of support mechanisms (e.g., Kiefer et al. 1986; Arkani-Hamed 1996; Simons et al. 1997; Yang et al. 2016); further, Venus may not be a perfectly archetypal stagnant lid planet and be better described instead by a plutonic-squishy lid regime (Lourenço et al. 2020). Option three is to be consistent with the pure dynamic topography we already produced to feed our scaling functions: we derive time-averaged power spectral densities from the numerical topography profiles, to which we fit a generic model.

Although the present study only considers dynamic topography, this same framework could be applied to any kind of topography on a planet as long as we can infer its spectral distribution.

The products of numerically modeled chaotic stagnant lid convection include time-dependent, dimensionless temperature fields and surface dynamic topography profiles (Figure 2). For each case, temporally and horizontally averaged temperature fields are used to calculate ${T}_{i}^{{\prime} }$, Ra_i , and other convective parameters; full outputs can be found in Table B1 in Appendix B. Average ${T}^{{\prime} }$ profiles hardly vary in time; hence, neither does ${T}_{i}^{{\prime} }$ nor the average position of the upper thermal boundary layer's base. Stepping up Ra₁ thins ${\delta }_{\mathrm{rh}}^{{\prime} }$ and lowers the rms height of topography in the regime we explore numerically. Increasing Δη thickens the stagnant lid because high viscosities are reached at lower depths; this is also associated with a slight increase in ${\delta }_{\mathrm{rh}}^{{\prime} }$.

Zoom In Zoom Out Reset image size

3.1.1. Fit to rms Height of Topography

Figure 3 shows the four-parameter linear fit between log(${h}_{\mathrm{rms}}^{{\prime} }$), log(Ra_i ), and $\mathrm{ln}({\rm{\Delta }}\eta )$, using the functional form in Equation (5). Best-fit parameter values and standard deviations are given in Table 3. The residual variance of this fit is ${\sigma }_{\mathrm{res}}^{2}\sim {10}^{-3}$, equal to the sum of squares error divided by the degrees of freedom. Because the fitted data correspond to averages over model time, the standard errors of the mean independent and dependent variables are all small and do not impact the regression.

Zoom In Zoom Out Reset image size

Table 3. Topography Scaling Coefficients and Their Errors Obtained from Fitting a Multiple Linear Regression Model with an Interaction Term to Equation (5)

	A	B	C	D
Best fit	9.581	−0.5818	−1.510	0.07536
Standard deviation	3.298	0.1859	0.4220	0.02379
${\sigma }_{\mathrm{res}}^{2}=1.584\times {10}^{-3}$

Note. The bottom row reports the residual variance, ${\sigma }_{\mathrm{res}}^{2}$, of the fit.

Download table as:  ASCIITypeset image

The key piece of information from this section is that chaotic convection with temperature-dependent viscosity does not lend itself to constant power-law scalings of ${h}_{\mathrm{rms}}^{{\prime} }$ with Ra_i (or Ra₁). The value of Δη is effectively altering the slope of $\mathrm{log}({h}_{\mathrm{rms}}^{{\prime} })$ with $\mathrm{log}({\mathrm{Ra}}_{i})$. Smaller viscosity contrasts of 10⁷ (b = 16) and below are associated with strongly negative slopes. With increasing Δη, the slope grows systematically shallower, until it changes sign between Δη = 10⁸ (b = 18) and Δη = 10⁹ (b = 20). Conversely, the effect of Ra_i on Δη is such that at higher Ra_i above ∼6 × 10⁷, large viscosity contrasts favor high rms topography, while at lower Ra_i below ∼6 × 10⁷, small viscosity contrasts favor high rms topography. At Ra_i ∼ 6 ×10⁷, these slopes "cross over" and the effect of Δη disappears.

Evidently this behavior is governed by a complex, chaotic system; extracting a general mechanistic understanding is compromised by the limited number of runs performed here. The effect of Δη to increase ${h}_{\mathrm{rms}}^{{\prime} }$ may be related to thermal isostatic uplift within the stagnant lid (Kucinskas & Turcotte 1994; Moore & Schubert 1995; Orth & Solomatov 2011). We include thermal isostasy as part of the full dynamic topography. Under a swell, hot, low-density upwelling material extends to shallower depths. To compensate, the cold, dense overlying lithosphere grows thinner, and it is buoyed upward. It can be shown that the maximum amount of thinning is directly proportional to the average lithospheric thickness. Hence, higher-viscosity-contrast convection, with its deeper lid bases, will enable a greater magnitude of thermal thinning. Meanwhile, smaller Ra_i are associated with thicker δ_rh, to which dynamic topography should be proportional (Parsons & Daly 1983). (For a constant Δη, lowering Ra₁ also slightly increases δ_lid and thus the potential for thermal thinning.) We speculate that there is a trade-off whereby the Δη effect dominates when stagnant lids are already thick and when convection is too vigorous to support high topography in its thin thermal boundary layers. Conversely, for lids that are not particularly thick, Ra_i (and δ_rh) becomes more relevant.

A corollary of this is that at the still-higher values of Ra_i expected for realistic rocky planets (up to several orders of magnitude beyond the range amenable to numerics; see discussion in Section 4.4.2), the sensitivity of ${h}_{\mathrm{rms}}^{{\prime} }$ to the viscosity scale becomes quite high indeed. If the absolute viscosity follows an exponential law, $\eta (T)\sim \exp (-{bT})$, high b is associated with low η for the same T, implying low ${h}_{\mathrm{rms}}^{{\prime} }$.

3.2.1. Thermal Evolution

Underlying thermal histories are sampled in Figure 4. Because all test planets are initialized at quasi-equilibriated temperatures and stagnant lid thickness, their evolutionary paths reflect secular cooling alone, which track roughly parallel at around −100 K Gyr⁻¹. Radiogenic heating inevitably declines with age, with surface heat losses lagging behind slightly; the present-day Urey ratios are ∼0.65 depending on planet mass.

Zoom In Zoom Out Reset image size

Interior temperatures and Ra_i increase with M_p as anticipated from simple scaling laws. We expect the heat flux q_u to increase linearly with planet radius for a fixed internal heat generation rate. This implies that ${q}_{u}\propto {M}_{p}^{1/3}$, ignoring compression. We can rewrite Equations (11)–(15) as

$$\begin{eqnarray}&&{\eta }_{u}=\displaystyle \frac{{\rho }_{m}{g}^{u}{\alpha }_{m}{k}_{m}^{3}{a}_{\mathrm{rh}}^{4}{\rm{\Delta }}{T}_{\nu }^{4}}{{\kappa }_{m}{q}_{u}^{3}{\mathrm{Ra}}_{\mathrm{crit}}^{u}},\end{eqnarray} \tag{ 20 }$$

Thus, we have ${\eta }_{u}\propto {M}_{p}^{-1};$ Equation (17) leads to Ra${}_{i}\propto {M}_{p}^{2}$ for approximately the same temperature difference. A five times more massive planet has a 25 times larger Ra_i (see also Stevenson 2003; Kite et al. 2009).

Figure 4 illustrates how uncertainty in the viscosity law parameters E_a and η₀ affects the spread and mean behavior of the dimensional h_rms and its physical constituents over time. Temperature-dependent viscosity exhibits self-regulating behavior: a slight increase in temperature lowers the viscosity, hence more vigorous convection via Equation (1). This leads to more efficient heat loss out of the top of the convecting cell, lowering temperatures in turn. This positive feedback is not visible in a single run (which are already at quasi-steady-state in our case), but we do see the effect at play over the entire ensemble: its range of η_m (t) is always less than an order of magnitude, despite a three-order range in η₀. Meanwhile, T_m is adjusting such that q_u approaches a balanced state for a given q_rad and surface-area-to-volume ratio. Hence, the rheological uncertainty manifests itself in T_m .

We note that these calculated Ra_i values are on average higher for a given M_p than those commonly associated with Venus or Mars. The thermal Rayleigh numbers of real planets require some dexterity to extract, but the few constraints available suggest a value on the order of 10⁶ for Mars (Kiefer 2003; Samuel et al. 2019). Constraints for Venus are even more scarce, but previous work employs Ra at upper mantle temperatures on the order of 10⁷–10⁸ (Huang et al. 2013; King 2018). This discrepancy is partly explained by the more viscous mantles we permit in this exoplanet study. Further caveats to our Ra_i estimates are discussed in Section 4.4.3.

The dimensional h_rms reflects a trade-off between b, Ra_i , and the dimensionalization factor dΔT α_m through Equations (2) and (5). Extrapolating Figure 3 would imply that, in the b-Ra_i regime of the 1D models, high h_rms is favored with high b and low Ra_i . Thus, deep, hot, weak mantles are doubly inhibited from having any remarkable topography. It is clear from Figure 4 that deeper mantles are not enough to make up for lost ${h}_{\mathrm{rms}}^{{\prime} }$.

Ultimately, the thermal state plays a main role in limiting the amplitude of dynamic topography. Hotter mantles necessitate lower viscosities, more vigorous convection, and thinner thermal boundary layers. Within these thinner boundary layers, there may be less scope for density variations related to thermal expansion. If we know some property of a planet to have a strong effect on its interior temperatures, then we might expect it to also impact its dynamic topography.

3.2.2. Dynamic Topography as a Function of Bulk Exoplanetary Properties

We now test the topographic reaction to planet age, mass, CMF, and radioisotope budget (Figure 5). We find h_rms to decrease with M_p and χ_rad and increase with age and CMF. Assuming that the x-axes in this figure cover the limits within which we expect to find most rocky exoplanets, then it is plausible that the resulting y range marks the variability of pure dynamic topography that nature could manifest, if our scaling relationship indeed applies. The fact that h_rms drops by the largest absolute amounts over M_p and χ_rad reflects the geodynamic significance of these parameters, as well as the spread over which we would expect to find rocky planets. The senses of change of h_rms with M_p and χ_rad are predictable from their known effects on T_m . That is, hotter interiors are expected for massive, U- and Th-rich planets, hence lower h_rms. Uncertainty in h_rms predictions is tied to uncertainty around the underlying thermal histories—yet another clue to the immeasurable usefulness of characterizing this uncertainty more rigorously (e.g., Seales & Lenardic 2020).

Zoom In Zoom Out Reset image size

The raw values of h_rms predicted by our scaling relationship are on the order of hundreds of meters, while the hottest planets can exhibit mere tens of meters of dynamic topography. In fact, due to inherent self-regulation, it is difficult to achieve significantly higher topographies in our 1D model while keeping to Earth-like values of the free parameters. This result may seem very low when compared to the heights of typical topographic features seen across the solar system. However, a fair comparison requires isolating an rms height of just the dynamic component of topography; this is not model independent, as we will discuss (Section 4.3.4).

We have tested three fiducial spectral models to find a relationship between the rms and peak value of dynamic topography. The theoretical red-noise model, the empirical Venus model, and the numerical dynamic topography model all produce an h_peak that is, on average, some constant scalar multiple of h_rms. For both numerical dynamic topography and the total Venus topography, h_peak ≈ 3.5h_rms, and for red-noise topography, h_peak ≈ 3.9h_rms. (For a pink noise structure similar to Earth's observed dynamic topography, h_peak ≈ 4.0h_rms.)

We use our h_peak estimations to derive the ocean basin volume capacity V_cap as a function of planet mass (Figure 6). This quantity represents the smallest volume of surface liquid water that would entirely inundate a planet. The actual land fraction requires knowing the ocean mass. We leave sea level as an unknown quantity and simply consider fiducial surface water budget scenarios. Specifically, we treat the amount of surface water as a constant mass fraction of M_p . This parameterization brackets the planet's total water budget with its volatile partitioning between the interior and exterior—in reality the amount of water stored in the mantle would affect the planet's thermal evolution through its rheology (and melting history, which is not modeled).

Zoom In Zoom Out Reset image size

As the basin volume capacity changes with M_p , so too does the water volume corresponding to this mass fraction (we assume a density of 1000 kg m⁻³; salt water is slightly denser). Figure 6 can be read as follows: for a given surface water budget, the planet mass where this contour intersects the basin capacity gives the most massive planet that could sustain land with dynamic topography alone. For example, a 1 M_⊕, 4.5 Gyr old planet endowed with solar U and Th could hold about 0.3 Earth oceans on its surface. The internal heating rate has a strong influence on V_cap.

Figure 6 compares different assumptions about the spectral distribution of topography, which would affect the relationship between the peak and rms topography. The dynamic topography and Venus models overlap identically, and the red-noise spectrum results in only slightly larger V_cap, seemingly because they are very similar in the low-degree regions where most of their power is concentrated. The basin volume corresponding to an infinitesimally small but nonzero land area is insensitive to the distribution of topography at high frequencies.

We can formulate these results in terms of a simple scaling analysis. Equation (19) can be written as ${M}_{\mathrm{cap}}\,=4\pi {R}_{p}^{2}{\rho }_{w}{\rho }_{m}/({\rho }_{m}-{\rho }_{w}){h}_{\mathrm{peak}}$, where M_cap is the ocean basin capacity in kg. For Earth's ocean mass (1.4 × 10²¹ kg), this means that a peak topography h_peak less than 2.7 km leads to a waterworld. If h_peak were independent of planet mass, we would expect ${M}_{\mathrm{cap}}\propto {M}_{p}^{2/3}$ owing to the increase in surface area alone (the mass–radius relation in Equation (7) gives a slightly shallower power owing to compression). However, we have h_peak strongly decreasing with increasing mass. For dry olivine and solar U and Th abundances, ${h}_{\mathrm{peak}}\propto {h}_{\mathrm{rms}}\propto {M}_{p}^{-0.5}$. From Equation (19), ${M}_{\mathrm{cap}}\propto {R}_{p}^{2}{h}_{\mathrm{peak}}$, so ${M}_{\mathrm{cap}}\propto {M}_{p}^{0.04}$ using Equation (7). Warmer, less viscous interiors decrease this exponent, so the most massive rocky planets have the smallest basin capacities even though they have the largest surface areas. If the pressure of a topographic load is balanced only by a constant compressive strength of the crust rock, we have h_peak ∝ g⁻¹, and the resulting proportionality ${M}_{\mathrm{cap}}\propto {M}_{p}^{0.08}$ is also quite flat (though the overall basin capacity would be higher). We are being conservative about how likely planets are to have dry land by considering only dynamic topography.

It is important to emphasize that the basin capacities shown in Figure 6, based on dynamic topography alone, are likely underestimating the true value. The observed topographies of Venus, Earth, and Mars produce basin capacities of 3.4, 3.3, and 2.9 Earth oceans, respectively, whereas the model produces basin capacities of <1 Earth ocean. The peak and rms elevations of our terrestrial planets are much higher than those predicted by the dynamic topography scaling here. Other mechanisms contribute to supporting higher topography on planets. Also, our model may underpredict dynamic topography for a given planet mass, as we will discuss in the next section.

Figure 6 suggests that reasonable changes to the spectral distribution of topography have no strong effect on how peak dynamic topography scales with planet mass, and hence on the volume of water that could be contained below this highest point. Our concern with topography's characteristic harmonic structures might thus seem somewhat tangential to (or in the worst case, distracting from) the basic problem that this study purports to address. However, these details would become more of a concern if the field can mature—and especially if we hope, someday, to use informed topography distributions as a boundary condition in exoplanet climate models (e.g., Turbet et al. 2016; Rushby et al. 2019). For example, the volume calculated in Equation (19) represents the amount of water that would flood a planet exactly, leaving just an island with infinitesimally small area. Yet in principle one could also calculate the maximum basin size associated with any arbitrary land fraction. These intermediate land fractions may be much more sensitive to spectral complexities, such as wide plains or anisotropic mountain ranges.

The initial questions here have justified simplified harmonic structures of topography as such. Specifically, we have presumed a loglinear model of the power spectral density, which is to say that the variance of elevation is a power-law function of the horizontal distance scale, and that this relationship is constant over the whole planet (as proposed in, e.g., Turcotte 1987). Contemporary workers now know the behavior to be much more nuanced. Local estimates of topography's spectral slope can appear notably inconstant—the surface roughness is heterogeneous—but these differences are entwined by further power laws of other statistical moments, out to virtually infinite order, all culminating neatly in a mathematical model with three scale-invariant parameters (e.g., Pelletier 1999; Gagnon et al. 2006; Lovejoy & Schertzer 2007; Ali Saberi 2013; Liucci & Melelli 2017; Rak et al. 2018; Landais et al. 2019a; Keylock et al. 2020). Landais et al. (2019b) have demonstrated the use of such a descriptive model for synthesizing surface relief of arbitrary rocky planets. Thus, the framework exists for representing full global topography layouts to a high degree of statistical realism and with few parameters. The hitch is that these parameters are empirical on a case-by-case basis: the gain in descriptive accuracy may not translate to predictive power for distant exoplanets. At present there is no theory tying the pattern to the (geophysical) process. If this gap could be bridged with more work based on Earth and solar system bodies, then these realistic mathematical models could be applied, and higher-order insight about the topographies of exoplanets might not necessarily be a fantasy.

Any deterministic prediction of h_rms will be hindered by the unknown mantle rheology. Increasing the activation energy of viscosity from 240 to 300 kJ mol⁻¹ will double h_rms for an Earth-mass planet, all else being equal. This uncertainty propagation is built into our model via the scaling functional form in Equation (5). E_a enters this equation twice, in both b and Ra_i (via η_m ). Particularly in the high-Ra_i regime, small changes in the viscosity contrast parameter b create large changes in ${h}_{\mathrm{rms}}^{{\prime} }$ (Figure 3).

We have attempted to capture some of the rheological uncertainty by varying E_a and η₀, the free parameters in the Arrhenius viscosity law (Equation (9)). However, we cannot claim that our results are propagating nature's true variability. First, the underlying covariance of these parameters is not known. The prior range employed by our study covers only pure olivine and pure orthopyroxene, and roughly so at that. Spaargaren et al. (2020) also parameterize the mineralogical control on viscosity with an extra prefactor that increases over three orders of magnitude, calibrated between ferropericlase-rich (high Mg/Si) and stishovite-rich (low Mg/Si) lower mantle compositions (Ballmer et al. 2017; Xu et al. 2017). Relating the rheological parameters to the lower or upper mantle composition in a realistic way requires not only a complex thermodynamic model predicting these mineral compositions but also a data set of strain rates from experiments and ab initio mineral physics. The actual strain rate of an olivine-orthopyroxene aggregate is certainly not a simple combination of diffusion creep flow laws. Further, in practice, real mantle viscosities will be strongly sensitive to their water content, unlikely to ever be known for a given exoplanet.

The second reason why we are not capturing the true variation is that our fixed rheological model ignores structural uncertainty by design. We have only considered diffusion creep with no pressure dependence, but the creep mechanism depends on shear stress and is not known a priori. Including pressure dependence in the parameterization (with adiabatic profiles from an interior structure model, for example) would lead to higher viscosities and sluggish flow in the lower mantle. Importantly, and in particular for more massive planets, this fact could render the viscosity self-regulation less efficient (Stamenković et al. 2012), meaning that internal temperatures for evolved planets become much more sensitive to initial temperature conditions, and the resulting h_rms scatters more widely (overall, retaining a hotter mantle at older ages will reduce h_rms). Uncertainty would grow severer still if one allowed for complex rheological features such as a low-viscosity asthenosphere (Bodur & Rey 2019), which manifests in smaller-scale dynamic topography on Earth (Hoggard et al. 2016). Finally, technically, the lithosphere itself obeys a distinct viscoelastic rheology, and coupling these dynamics to a convection model would also modify its topography amplitudes (Patočka et al. 2017)—we have ignored elastic effects in this attempt (Section 4.3.3).

All this rheological uncertainty is worth discussing because dynamic topography is apparently sensitive to both viscosity's absolute value and how it changes over the boundary layers (Hager et al. 1989). Low viscosities imply higher temperatures and low convective stresses. For the isoviscous case, the association of low viscosity with low topography can be seen clearly in Table 2 of Lees et al. (2020), from which we get a numerical scaling of ${h}_{\mathrm{rms}}^{{\prime} }\propto {\eta }^{-0.6}$, with interior temperature and lithospheric thickness fixed. If we have two isoviscous layers with a stiffer top layer (i.e., approximating a cool viscous lithosphere), then there is an analytical solution for the surface normal stress induced by a spherical density anomaly at some depth (Equation (34) in Morgan 1965). In this solution, the effect of relative viscosity is strongest when density anomalies are nearer the surface.

In determining a scaling relationship for the rms and peak amplitudes of dynamic topography from numerical convection, we have assumed that details of our methodology can produce generalizable results. This section discusses some important caveats.

4.3.1. Low-order Power

4.3.2. Geometry and Heating Mode Effects

4.3.3. Filtering in the Lithosphere

4.3.4. Paucity of Ground Truths

4.4.1. Sensitivity to Functional Form

4.4.2. Extrapolation across Rayleigh Numbers

4.4.3. Accuracy of Interior Rayleigh Number Estimates

4.4.4. Model Validity at High Planet Mass

Agol et al. (2021) give preliminary constraints on the surface water content of the TRAPPIST-1 planetary system, for different values of the CMF and assuming that all water exists as a condensed surface layer. Although the problem is degenerate, planets e–g appear consistent with water layers deeper than Earth's, on the order of at least 0.1% of the planet mass. Other independent estimates have produced similar results (Acuña et al. 2021). This water budget would place TRAPPIST-1e to g well above the upper water mass limit for maintaining land with dynamic topography. Note, however, that the high rates of tidal heating inferred for some of these planets (Barr et al. 2018) would reduce dynamic topography beyond what is modeled here.

As we have previously emphasized, however, the true limit to elevation differences on a planet will be higher than that suggested by purely dynamic topography. To estimate a planet's total scope for land, we can calculate the minimum value of h_peak required for an instance of land on a planet with a given radius and surface water content. We find that any instance of land on TRAPPIST-1e would require a peak topographic amplitude of ∼40 km (a minimum rms topography of ∼10 km), given 0.3 wt.% surface water (Agol et al. 2021 estimate for a CMF of 0.25). Then, one could compare this minimum to a rough estimate of the overall maximum elevation.

In Section 1 we motivated a crustal strength limit: for a surface load of ρ gh, somewhere in the crust below, at a depth of about 1/4 times the load width, a minimum stress difference Y of 1/2 to 1/3 ρ gh is sustained (Jeffreys 1929). This result assumes a flat earth model of elastic stress distributions and holds for various load configurations of horizontal scale less than a few hundred kilometers. Melosh (2011) illustrates that the force balance given by

$$\begin{eqnarray}&&Y\approx 0.5{\rho }_{c}{gh},\end{eqnarray} \tag{ 21 }$$

with a crust density ρ_c = 2700 kg m⁻³ and Y set at an effective crustal strength on the order of 100 MPa, will roughly reproduce the maximum elevations of Venus, Earth, and Mars (Figure 7). While this estimate is certainly an oversimplification, a more rigorous effort will naturally become very complicated, due not least to the difficulty in predicting, from planetary bulk properties, a value of Y corresponding to the maximum h.

Zoom In Zoom Out Reset image size

In typical crustal strength models, the strength increases with depth (lithostatic pressure) according to the rock's resistance to frictional sliding in the relatively cool, shallow part—the brittle regime—until viscosity is low enough to favor ductile deformation instead, and strength starts to decrease with depth (temperature). Thus, the strongest part of the crust is near this brittle−ductile transition. However, the resulting strength maxima of ∼500 MPa or more for Earth-like conditions (Katayama 2021) would imply ∼40 km of peak topography using Equation (21); it is a limit not necessarily reached in practice. Further complicating the application of Equation (21), crustal strength profiles are strongly sensitive to the temperature profile and porosity of the crust—generating these profiles for arbitrary exoplanets must attend to assumptions on these facets (Byrne et al. 2021)—and surface gravity has a nonlinear effect on brittle strength through its influence on porosity and fracture density (Heap et al. 2017). For example, doubling the thermal gradient will approximately halve the maximum Y—and thus h—in a dry case, and including hydrostatic pore fluid pressure shows a similar decrease (Katayama 2021).

A parallel approach to estimating maximum elevation differences from crustal concerns comes from isostasy. The height of a topographic feature above a plain is h_A = (t_R − t_avg)(ρ_m − ρ_c )/ρ_c , where t_R is the thickness of the crust below the feature and t_avg is the average crustal thickness of the plain. For a basaltic crust (the primary crust formed from an Earth-like bulk composition), the maximum value of t_R is set by the phase transition from basalt to denser and unstable eclogite: the crust cannot be much thicker than the depth of this transition. This fact limits the peak isostatic h_A to about 15 km for a Venus-like case (Jull & Arkani-Hamed 1995). However, the depth of this phase transition depends sensitively on the crust thermal structure, and estimating h_A in practice requires knowing t_avg.

Finally, the height limits of volcanoes in particular must follow tighter rules. Magma will only rise to the top of a vent—and contribute to a growing pile of lava—so long as the vertical pressure gradient across the system is positive. Castruccio et al. (2017) write this limit as ${h}_{\max }=({\rm{\Delta }}\rho /{\rho }_{m})H+{\rm{\Delta }}{P}_{i}/({\rho }_{m}g)$, where Δρ is the density contrast between the crust and the magma, H is the depth from the surface to the magma chamber, and Δ_Pi is the critical overpressure to trigger an eruption (the pressure that would crack the magma chamber roof, related to the tensile strength of the crust and tellurically on the order of ∼20 MPa). Although a narrow range of H can be argued for on Earth, related to the magma water content and crustal rheology (Huber et al. 2019), this concept has not yet been expanded to comparative planetology.

In light of the above complexities, it is difficult to find a middle ground between the oversimplification of Equation (21) using a universal crustal strength estimate and a careful case-by-case application (Barton & Shen 2018). We will employ the former for the present purpose of comparing peak dynamic topography to peak total topography. We consider Y = 100 MPa, ostensibly representing the compressive strength of granite—the difference in compressive strength between an average granite and average basalt seems to be smaller than the spread seen across individual basalt samples in various laboratory conditions (Heap et al. 2017)—but also include scalings for half and double this strength value.

Figure 7 plots the containable ocean mass fraction scalings corresponding to both this crust strength limit and the dynamic topography limits calculated previously. For a given scaling relationship, points above the line would be waterworlds. We see that planet e may have coexisting land and water if its crust could withstand around 200 MPa of normal stress. Although these strengths can be achieved on Earth, it is not immediately obvious that they would be available at the right loci. Note that it is very difficult in practice to put a lower limit on these water budgets. Nevertheless, according to Agol et al. (2021), TRAPPIST-1e through g could easily be wet enough that estimating their land propensities may seem moot. However, our growing catalog of planets may soon present a case study closer to the waterworld−land world transition.

Another takeaway from Figure 7 is that for the most massive rocky planets amplitudes of dynamic topography in the most favorable case seem to approach the overall limit. Scalings for different internal heating scenarios have different slopes because, as surface heat fluxes increase with the surface-area-to-volume ratio, larger planets are penalized such that any extra radiogenic heat would escape less easily. Thus, more internal heating per unit volume in more massive planets will have a more drastic effect on topography.

At the moment, it is not guaranteed that constraints on these or any rocky exoplanet water budgets could be tightened much in the future. With current Bayesian inference methods, uncertainties on retrieved water mass fractions may be capped at around σ ≈ 10 wt.%, independent of the observational uncertainty on the planet mass itself (Otegi et al. 2020). Meanwhile, topography can avert waterworlds only for water mass fractions of ≲1 wt.%. Therefore, with respect to predictions about a given exoplanet, any topographic contribution to land coverage could be washed out by the uncertainty on the inferred water budget.

In Figure 5, we predicted how dynamic topography might vary as a function of several properties broadly deemed observable. None of these properties will be perfectly known, or even necessarily constrained well enough such that they are not the dominant source of uncertainty, but we will leave a more detailed assessment of this uncertainty to future work.

In any case, an obvious fact emerging from our scaling law application is that there is a pivotal future role to be filled for any constraints on rocky planet compositions. This study provides yet another example of how higher-order properties of planetary interiors govern their surface character. Namely, mantle viscosities, radiogenic heating rates, and core mass fractions all relate to planetary ratios of certain major elements: viscosities decrease with Mg/Si, radiogenic heating rates increase with U/Si and Th/Si, and core mass fractions increase with Fe/O. Exoplanet compositional parameters are not completely inaccessible because refractory element ratios are expected to generally preserve themselves between a star and its planets (Thiabaud et al. 2015; Hinkel & Unterborn 2018; Putirka & Rarick 2019; Adibekyan et al. 2021). Although pilot work is surely needed, this useful fact means that element abundances from stellar spectra offer a promising constraint on planetary interior dynamics. Additionally, measurements of the same element ratios in polluted white dwarf spectra could inform the underlying natural distributions of bulk rocky planet composition across nearby star systems (Bonsor et al. 2021).

Observables for exoplanetary topography itself would be buried quite deep. McTier & Kipping (2018) proposed that extreme topographic features could induce scatter in an exoplanet's transit photometry, but the associated signal would not be detectable with realistic photometric precision. Proposed next-generation direct imaging missions might be capable of enough precision for the exo-cartography of small planets—solving the inverse problem of 2D albedo distributions from time-resolved light curves—which might discriminate between land and ocean surfaces (Cowan & Fujii 2018; Farr et al. 2018; Lustig-Yaeger et al. 2018; Aizawa et al. 2020; Kawahara 2020). Interpretations of the data may remain highly model dependent and burdened by cloud removal, however (Paradise et al. 2021; Teinturier et al. 2022). Ocean fractions might also be discerned from near-infrared polarimetric observations (Takahashi et al. 2021). A land fraction between zero and unity would necessitate some surface roughness, leading to an upper limit on the water budget given some inferences about topographic propensity.

We choose our Case 4 simulation (Table 1) from which to extract a scalable model spectrum of the surface dynamic topography, since its temporal distribution of ${h}_{\mathrm{rms}}^{{\prime} }$ is the narrowest. A type-2 orthonormalized discrete cosine transform of this profile produces a Fourier representation,

$$\begin{eqnarray}\begin{array}{rcl}{f}_{p} & = & 2\gamma \displaystyle \sum _{n=0}^{N-1}{h}_{n}^{{\prime} }\cos \left(\displaystyle \frac{\pi p(2n+1)}{2N}\right),\\ \gamma & = & \left\{\begin{array}{ll}\sqrt{\tfrac{1}{4N}}, & \mathrm{if}\ p=0\\ \sqrt{\tfrac{1}{2N}}, & \mathrm{otherwise},\end{array}\right.\end{array}\end{eqnarray} \tag{ A1 }$$

from which we can find a 1D power spectral density,

$$\begin{eqnarray}&&{\phi }_{0}^{1{\rm{D}}}=2{\rm{\Delta }}{x}^{{\prime} }{\left({f}_{p}\right)}^{2},\end{eqnarray} \tag{ A2 }$$

as a function of dimensionless wavenumber,

$$\begin{eqnarray}&&{k}^{{\prime} }=\displaystyle \frac{\pi }{{L}^{{\prime} }}p,\end{eqnarray} \tag{ A3 }$$

where ${h}_{n}^{{\prime} }$ is the height of dynamic topography at sample point n, N is the number of sample points in the spatial profile (fixed by the mesh size), p = [0,...,N − 1], ${L}^{{\prime} }=8$ is the dimensionless box width, and ${\rm{\Delta }}{x}^{{\prime} }={L}^{{\prime} }/N$. We calculate ${\phi }_{0}^{1{\rm{D}}}$ at every model time step and use the average for our baseline spectrum. This spectrum has an rms amplitude ${h}_{\mathrm{rms},\ 0}^{{\prime} }$.

There is an upper wavenumber limit, ${k}_{\max }^{{\prime} }$, at around the equivalent wavelength of the upper thermal boundary layer thickness, δ_rh, where features narrower than this are not meaningful for the dynamic topography. We also observe all spectra roughly rolling off to a constant value at wavenumbers below around twice the convection cell depth, so we set ${k}_{\min }^{{\prime} }=2d$. In log-log space, ${\phi }_{0}^{1{\rm{D}}}$ is approximately linear from ${k}_{\min }^{{\prime} }$ to ${k}_{\max }^{{\prime} }$. Therefore, we approximate the power spectra by two line segments. We fit a constant slope between ${k}_{\min }^{{\prime} }$ and ${k}_{\max }^{{\prime} }$ and assign a value of ${\phi }_{0}^{1D}({k}_{\min }^{{\prime} })$ wherever ${k}^{{\prime} }\lt {k}_{\min }^{{\prime} }$. This fit is done to the average power spectral density over all time steps for the given simulation. We interpolate this fitted function such that it has a discrete value at each integer spherical harmonic degree l, where $l={k}^{{\prime} }{R}_{p}^{{\prime} }$ − 0.5, from l = 1 to the nearest degree to ${k}_{\max }^{{\prime} }$. That is, we do not scale ${k}_{\max }^{{\prime} }$. While realistically ${k}_{\max }^{{\prime} }$ would increase with Ra₁, the effect on ${h}_{\mathrm{rms}}^{{\prime} }$ is small (less than one part in a thousand) because these high wavenumber bands hold such little relative power. For this generic spectrum we assume a dimensionless planet radius ${R}_{p}^{{\prime} }=2$ (a core radius fraction of 0.5 for a dimensionless mantle depth of 1; varying ${R}_{p}^{{\prime} }$ has negligible effects on the results).

Figure A1 shows the 1D power spectral densities ${\phi }_{h}^{\mathrm{PSD}}$ of dynamic topography computed from our 2D numerical modeling experiments, normalized as a percentage of the total power. Between ${k}_{\min }^{{\prime} }$ and ${k}_{\max }^{{\prime} }$, the loglinear slopes of the topography spectra are roughly similar within the noise for all Ra₁, Δη cases. Due to our limited number of 2D runs, however, we cannot really make a compelling case for this statement, and we would not back our interim result outside of its intended, rather inconsequential usage here. For example, we might expect more vigorous, higher-Ra convection to exhibit more smaller-scale drips from the upper thermal boundary layer, leading to slightly more topographic power at high wavenumbers—although the total power would be virtually unaffected by these high-frequency features. Note also that because the spatial domain topography is 1D, data paucity will always entail a certain amount of noise, compared to a 2D grid of topography from a 3D convection simulation.

Zoom In Zoom Out Reset image size

Also in Figure A1 is the observed topography spectrum of Venus from Wieczorek (2015). On Venus, elastic and compositional sources of topography are superimposed on dynamic topography. Venus's spectrum thus provides an empirical modification of the pure dynamic topography. As a third and final spectral model, we have the theoretical red-noise spectrum given by the power law ${\phi }_{h}^{\mathrm{PSD}}\propto {k}^{-2}$ and a roll-off wavenumber the same as the numerical spectrum. Compared to the numerical dynamic topography, Venusian topography and red-noise topography both have a shallower slope and retain more power at higher wavenumbers—as expected from the high-frequency nature of topography created by impact cratering and volcanism. The Venus spectrum additionally shows a peak at degree l = 3. Note that these (normalized) spectra represent different geophysical and geomorphologic processes and are therefore not expected to have the same absolute rms value.

We use the pyshtools.SHCoeffs.from_random() function to obtain a set of spherical harmonic coefficients consistent with ${\phi }_{0}^{1D}$ (Wieczorek & Meschede 2018). This function requires a power per l (dimensional units m²), so we apply a conversion from ${\phi }_{0}^{1D}$ (dimensional units m² m). First, we find the effective 2D power spectral density assuming radial symmetry, ${\phi }_{\mathrm{iso}}^{2D}$ (dimensional units m² m²), which would correspond to our 1D spectrum:

$$\begin{eqnarray}&&{\phi }_{\mathrm{iso}}^{2D}=\displaystyle \frac{1}{{k}^{{\prime} }}{\phi }_{0}^{1D}.\end{eqnarray} \tag{ A4 }$$

The power per l is

$$\begin{eqnarray}&&{S}_{l}=\displaystyle \frac{{\phi }_{\mathrm{iso}}^{2D}\left(2l+1\right)}{4\pi {R}_{p}^{{\prime} 2}}.\end{eqnarray} \tag{ A5 }$$

With these normalizations, the total power per coefficient,

$$\begin{eqnarray}&&{S}_{{lm}}=\displaystyle \frac{{S}_{l}}{2l+1},\end{eqnarray} \tag{ A6 }$$

is proportional to ${\phi }_{\mathrm{iso}}^{2D}$. In converting our spectra into 2D equivalents, we are presupposing that 2D Cartesian and 3D spherical models result in approximately similar topography power spectra with consistent ${h}_{\mathrm{rms}}^{{\prime} }$. Using the output from Lees et al. (2020), we have verified that constant-viscosity convection in Cartesian geometry indeed produces similar spectra between 2D and 3D, but the assumption remains a caveat until dedicated 3D spherical realizations can test it. Nevertheless, we already know that it is incorrect to try fitting a scaling function to 2D numerical h_peak directly—this quantity is certainly sensitive to details of the model setup, as we have mentioned in Section 2.1.2.

If we are seeking a spatial map of a hypothetical spectrum other than ${\phi }_{0}^{1D}$ (i.e., different rms value), we take advantage of the fact that numerical dynamic topography spectra will appear to have roughly consistent slopes between ${k}_{\min }^{{\prime} }$ and ${k}_{\max }^{{\prime} }$, and hence scale S_l appropriately,

$$\begin{eqnarray}&&\bar{{S}_{l}}={S}_{l}{\left(\displaystyle \frac{{h}_{\mathrm{rms},1}^{{\prime} }}{{h}_{\mathrm{rms},0}^{{\prime} }}\right)}^{2},\end{eqnarray} \tag{ A7 }$$

where ${h}_{\mathrm{rms},1}^{{\prime} }$ refers to the desired rms of the new spectrum.

We can now obtain our set of coefficients via pyshtools: random spherical harmonic coefficients are generated from a normal distribution with unit variance, subject to the strong assumption of no correlation between wavenumbers.

Then, we again use pyshtools to expand the random spherical harmonic coefficients onto a Gauss–Legendre quadrature grid. At this stage we can dimensionalize the spatial domain topography with Equation (2), given the results of the parameterized convection model. A sample elevation map is shown in Figure A2. Because the randomly generated spherical harmonic coefficients are not unique for a given power spectrum, we reduce the noise by generating 500 sets of coefficients and taking the average of the resulting peak elevation values.

Zoom In Zoom Out Reset image size