Assessing the Accuracy of 2‐D Planetary Evolution Models Against the 3‐D Sphere

Regardless of the steady increase of computing power during the last decades, numerical models in a 3D spherical shell are only used in specific setups to investigate the thermochemical convection in planetary interiors, while 2D geometries are typically favored in most exploratory studies involving a broad range of parameters. The 2D cylindrical and the more recent 2D spherical annulus geometries are predominantly used in this context, but the extent to how well they reproduce the 3D spherical shell results in comparison to each other and in which setup has not yet been extensively investigated. Here we performed a thorough and systematic study in order to assess which 2D geometry reproduces best the 3D spherical shell. In a first set of models, we investigated the effects of the geometry on thermal convection in steady‐state setups while varying a broad range of parameters. Additional thermal evolution models of three terrestrial bodies, namely Mercury, the Moon, and Mars, which have different interior structures, were used to compare the 2D and 3D geometries. Our investigations show that the 2D spherical annulus geometry provides results closer to models in a 3D spherical shell compared to the 2D cylindrical geometry. Our study indicates where acceptable differences can be expected when using a 2D instead of a 3D geometry and where to be cautious when interpreting the results.


Introduction
Geodynamic modeling is a powerful approach to investigate the dynamics of the mantle and lithosphere of terrestrial planets and to explore the evolution of their interior that is not directly observable.Such models vary in their complexity and often employ different geometrical approaches to investigate physical processes such as secular cooling, mantle melting, and the generation of a magnetic field.When applying these models to interpret specific observations of the Earth and other planets, care must be taken in particular for the choice of geometry (see Noack & Tosi, 2012, for an overview of geometries), as this may significantly impact quantities such as the mantle temperature, the convection velocity, and the heat flux of the simulations.
The role of two-dimensional studies in the field of thermochemical mantle convection modeling is still predominant despite an ever-increasing computing power.Although the formulation of 3D grids has seen improvements in previous years with the Yin-Yang grid (Kageyama & Sato, 2004), the spiral grid (Hüttig & Stemmer, 2008a), and the adaptive 3D grid (Kronbichler et al., 2012), among others, simulations with a full spherical shell geometry remain highly expensive in terms of computational power, hence making them inappropriate to study broad ranges of parameters or conduct large exploratory studies.As an alternative, geometrical analogues to the 3D spherical shell have been extensively used, namely the 2D spherical axi-symmetric (van Keken & Yuen, 1995) and the more popular 2D cylindrical annulus (Jarvis, 1993), which is typically referred to as the 2D cylindrical geometry or 2D cylinder.The 2D axi-symmetric geometry has been used in earlier studies of mantle convection (e.g., Jarvis et al., 1995;van Keken & Yuen, 1995), but in addition to the artificial boundaries formed by the poles which trap down-and up-wellings, an asymmetry between the polar and the equatorial regions exists (van Keken, 2001).The cylindrical geometry on the other hand, while resolving the problems of the artificial boundaries at the poles imposed by the axi-symmetric geometry, still exhibits an important drawback.The ratio of the two surfaces (the planetary surface and the core surface) is different in the 2D cylindrical geometry compared to the 3D spherical shell.This leads to a mismatch in heat flux values between these geometries, as the heat flux of the core mantle boundary (CMB) is underestimated and the surface heat flux is overestimated in the 2D cylindrical geometry compared to a 3D spherical shell with the same ratio between the core and planet radius (i.e., radius ratio).
In order to mitigate this problem, van Keken (2001) introduced a scaling of the 2D cylindrical geometry such that the ratio of outer and inner areas of the cylinder matches the ratio obtained for the spherical shell.This scaling, however, while correcting the surface ratio discrepancy of the cylinder, still uses the volume elements of a cylinder.Additionally, this scaling creates an artificially smaller core, which in turn modifies the convection pattern in the mantle, leading for example, to a crowding of the structures near the CMB, a behavior that would be less present in a 3D spherical shell using the original, non-scaled, radii.
To overcome this major drawback of the cylindrical geometry, another 2D geometry called "spherical annulus" has been proposed by Hernlund and Tackley (2008).This geometry effectively uses a second degree of curvature by considering the same volume elements as the 3D geometry.Since no scaling is necessary for this geometry, it keeps the same radius ratio as the 3D spherical shell.
In the study of Hernlund and Tackley (2008), the 2D spherical annulus showed promising results to approximate the 3D spherical shell geometry, with mean temperature and Nusselt number well reproduced for steady-state thermal convection calculations.While these results are highly valuable, they correspond only to the case of an Earth-like radius ratio and only consider thermal convection simulations in the Boussinesq approximation.
More recently, Guerrero et al. (2018) performed a more extensive study using the 2D spherical annulus for stagnant lid convection models and compared the temperature distribution between the 2D spherical annulus and the 3D spherical shell.They showed that the 2D spherical annulus systematically overestimates the mantle temperature, especially in low radius ratio cases.However, an extensive study investigating the ability of the 2D spherical annulus to reproduce results obtained in a 3D spherical shell and a systematic comparison with the 2D cylinder for various setups has never been conducted so far.
In this study, we present simulations of thermal convection in the 2D spherical annulus and compare the results to the 2D cylinder and the 3D spherical shell.First, we focus on simple steady-state convection models using the Boussinesq approximation.We vary the Rayleigh number, the radius ratio, and the heating mode for isoviscous cases and run additional temperature-dependent viscosity models to determine which of the two 2D geometries (i.e., cylinder or spherical annulus) is able to best reproduce the 3D spherical shell results.The set of equations that was used for this comparison is described in Section 2.1, the grid geometries are displayed in Section 2.2, and a description of the cases investigated here is available in Section 2.3.A detailed analysis of the results is presented in Section 2.4.
In a second step, we run thermal evolution models with the same geometries in three separate scenarios with Moon-like, Mars-like, and Mercury-like structures to investigate how well the 2D spherical annulus reproduces the results of the 3D spherical shell geometry.The three planetary bodies were chosen since they cover a wide range of interior structures with small and large core to planetary radius ratio and they are all thought to have been in a stagnant lid regime over their entire thermal history, which makes them comparable in terms of their tectonic regime (Breuer & Moore, 2015).In Section 3.1 we list the equations used for the thermal evolution models.A description of the employed parameters and setup of the models is given in Section 3.2.The results are described in Section 3.3, while additional post-processing is presented in Section 3.4.A discussion of the steady-state and thermal evolution models is presented in Section 4, followed by conclusions in Section 5.

Steady-State Mantle Convection
In a first set of calculations we focus on the comparison between steady-state models in 2D and 3D geometries.For this purpose we test a large number of parameter combinations for isoviscous as well as temperaturedependent viscosity scenarios.

Mathematical Model
We investigate the spatial and temporal evolution of mantle flow by solving the conservation equations of mass, momentum, and energy.Here, the conservation equations are scaled using the mantle thickness D as length scale, the temperature drop across the mantle ΔT as temperature scale, and D 2 κ with κ being the thermal diffusivity as time scale.A Table listing the scaling factors is available in Supporting Information S1 (Table S1).In the remainder of this study, all parameters are non-dimensional.By assuming a Newtonian rheology, an infinite Prandtl number, and considering the Boussinesq approximation, the non-dimensional conservation equations read (Schubert et al., 2001;van Zelst et al., 2022): In Equations 1-3, u is the velocity vector, η is the viscosity, T is the temperature, e r is the unit vector in radial direction, P is the dynamic pressure, and t is the time.
The parameter Ra denotes the thermal Rayleigh number, a non-dimensional number, which controls the vigor of the convection in the mantle.H is the internal heating rate of the mantle that is given by Ra , where Ra Q denotes the Rayleigh number associated with internal heating.The Rayleigh numbers Ra and Ra Q read: where ρ ref is the reference density, g ref is the reference gravitational acceleration, α ref is the reference thermal expansivity, κ ref is the reference thermal diffusivity, η ref is the reference viscosity, k ref is the reference thermal conductivity, and H dim is the internal heating rate per unit mass (W/kg).
For the steady-state models, we use a constant or temperature-dependent viscosity that follows the Frank-Kamenetskii approximation (Frank-Kamenetskii, 1969), which is a linearized form of the Arrhenius law: The parameter Δη T is the logarithm of the viscosity contrast due to temperature and T ref is the reference temperature at which a non-dimensional viscosity equal to 1 is attained.For the thermal evolution simulations presented further in this study, we use another parametrization of the viscosity (Equation 15), which is discussed more in-depth in Section 3.1.

Grid Geometries
We use the numerical code GAIA (Hüttig & Stemmer, 2008a, 2008b;Hüttig et al., 2013) to model the mantle convection in the interior of rocky planets.GAIA solves the conservation equations (Equations 1-3) in their dimensionless form in 2D and 3D geometries.For the 2D geometry, we use both the classical cylindrical geometry (van Keken, 2001) and the spherical annulus geometry following the approach of Hernlund and Tackley (2008).
In the 2D cylindrical geometry, the volume of the grid cells is typically formulated using the equations for a 2D cylinder.In comparison to a volume element in a 3D spherical shell (dV S ), a volume element in the 2D cylinder (dV Cyl ) is written as (Anton et al., 2012): Geochemistry, Geophysics, Geosystems where r is the radial distance, ϕ is the azimuthal angle, z is the height, θ S is the polar angle of the spherical coordinates, θ Cyl is the azimuthal angle of the cylindrical coordinates, and the superscripts S and Cyl indicate the spherical and the cylindrical geometry, respectively.
While the 2D cylindrical geometry is the intersection of a plane through a cylinder, the 2D spherical annulus is a geometry based on the intersection of a spherical shell with a plane passing through its center.The subsequent coordinate system used in the 2D spherical annulus results in the reduction of the traditional spherical system to the equatorial plane.A volume element in the spherical annulus geometry (dV SA ) is the integration of the classic spherical volume element dV S on θ S (Besserer, 2012): The total volume of a spherical shell V S tot of inner radius r i and outer radius r o can then be written as follows: Thus, the 2D spherical annulus has a second degree of curvature, and uses an effective three dimensional formulation for its element of volume, as shown in Figure 1.
The 2D cylindrical geometry is scaled according to the scaling introduced by van Keken (2001), where the inner and outer radii of the 2D cylindrical grid (i.e., the core radius and the planetary radius, respectively) are changed such that the ratio between the surface area of core and mantle is the same as the ratio obtained in a 3D spherical shell geometry.The equations used to correct the inner and outer radii of the 2D cylinder are the following: where r ic and r oc are the inner and the outer radius of the 2D cylinder, respectively.The inner and outer radii of the spherical shell are denoted by r is and r os , respectively.In the following, we will refer to this type of geometry that considers the scaling of the inner and outer radii as the "scaled cylinder geometry."

Case Definition
In the first part of this study, we performed steady-state simulations in order to investigate the effects of the ratio of the inner to outer radius and of heating modes on the results obtained with the 2D scaled cylinder, 2D spherical annulus, and 3D spherical shell geometry.Our aim is to compare 2D and 3D geometries and determine for which scenarios the 2D spherical annulus gives closer results to the 3D spherical shell compared to the 2D scaled cylinder.To this end, we use models heated from below (purely bottom-heated), from within (purely internally heated), and from both below and within (mixed-heated).We use an initial random perturbation of the temperature field with an amplitude of 5%, vary the Rayleigh number Ra of our simulations from 10 4 up to 10 8 , and the radius ratio f from 0.2 to 0.8 for our isoviscous setup.
For the 2D geometries, we use between 1.1 × 10 4 and 6.7 × 10 4 grid points for low Rayleigh number simulations and between 4.8 × 10 4 and 4.1 × 10 5 for simulations with a Rayleigh number higher than 10 6 .For the 3D geometry, we use between 2.04 × 10 6 and 2.94 × 10 6 grid points.A more in depth description of each grid and its associated lateral and radial resolution is available in Table S2 in Supporting Information S1.A short comparison of our results to the ones of Hernlund and Tackley (2008) for isoviscous steady-state cases is presented in Section S3 in Supporting Information S1, and shows that the 2D spherical annulus implemented in GAIA can accurately reproduce the diagnostic output quantities provided by Hernlund and Tackley (2008).
For our models we use free-slip boundary conditions.The temperature of the upper boundary T surf is set to zero, while the one of the lower boundary is set to one for the bottom-heated and mixed-heated cases.For the purely internally heated cases we use a zero heat flux at the core-mantle boundary (CMB).
The simulations are run until a statistical steady-state is reached.Then, output quantities such as the average temperature, root-mean-square velocity, and top temperature gradient are computed using an arithmetic average of the values obtained over the last 10% of the simulation.While for purely steady-state models this is the same as taking the last output, for quasi steady-state time-dependent or periodic models this ensures to retrieve representative average values.The top temperature gradient dT/dr is the temperature gradient at the top of the domain, calculated between the last two shells of the grid.
An additional, more complex set of simulations includes a temperature-dependent viscosity, which leads to the formation of a stagnant lid at the top of the convecting domain.These simulations use radius ratios representative for a Moon-like ( f = 0.2), Mars-like ( f = 0.5), and Mercury-like ( f = 0.8) interior.We use here thermal and radiogenic Ra numbers with values similar to those expected for planetary mantles (Laneuville et al., 2013;Plesa et al., 2015;Tosi, Grott, et al., 2013), that is, Ra = 5 × 10 6 and Ra Q = 5 × 10 7 (see values of Ra and Ra Q in Table 1).We use the Frank-Kamenetskii parametrization for the viscosity (Equation 5) and set a viscosity contrast Δη T to 10 8 at a reference temperature T ref of 0.5 to ensure that we are in a stagnant lid convection regime.In this setup, the total number of nodes used for the 2D grids lies between 4.8 × 10 4 and 2.8 × 10 5 , accounting for a radial resolution of 100 shells and a lateral resolution between 472 and 2,828 points per shell.For the 3D simulations we use a radial resolution of 70 shells and a lateral resolution of 40,962 points per shell, being equivalent to a total number of nodes of 2.9 × 10 6 .
In total, we run 180 isoviscous simulations and 36 temperature-dependent viscosity cases.All parameters and the grid resolutions for these steady-state simulations are listed in Table S2 in Supporting Information S1.

Isoviscous Convection
For the first set, consisting of isoviscous simulations, we provide a thorough and systematic comparison between the 3D spherical grid and the two 2D geometries, namely the spherical annulus and the scaled cylindrical geometry.
We provide a summary of the comparison in Figure 2, where the analysis of the 180 simulations has been divided into three subplots, one for each heating mode.Each subplot contains two rows, one for the 2D spherical annulus (first row) and another one for the 2D scaled cylindrical geometry (second row).The colors and the values in percent indicate the relative error to the 3D spherical shell results.The computation details of the relative error can be found in Supporting Information S1 (Section S4) along with tables containing the values for each simulation in CSV format.Figure 2 shows that the mean domain temperature of the 3D geometry is more accurately reproduced by the 2D spherical annulus geometry than by the cylindrical one for every heating mode, up to values of Ra = 10 7 .This difference is even more noticeable in the low radius ratio setups (i.e., f = 0.2 and 0.4).We also see for the purely internal heated cases and high Rayleigh number cases (i.e., Ra = 10 7 and Ra = 10 8 ) that this difference is even more dramatic with an overestimation of the temperature in the cylindrical geometry by up to ∼20% for a small radius ratio (i.e., f = 0.2), while it reduces to only a few percent in the 2D spherical annulus geometry.
In the case of the root mean square velocity, we see a general increase of the accuracy with the 2D spherical annulus geometry.This increase is less pronounced than what is observed for the mean temperature.We see a slightly better match for the purely bottom-heated and mixed-heated setups, where simulations with Ra <10 7 provide results that are slightly underestimated by the 2D scaled cylinder and slightly overestimated by the 2D spherical annulus.However, the cases with Ra ≥ 10 7 show a net improvement, with an underestimation decreasing for cases with an intermediate radius ratio (0.4 and 0.6).This improvement is yet less visible in the case of the internal heated cases where only the 0.2 and 0.8 radius-ratios show marginally better results.For the purely internally heated setup, however, almost no improvement is visible, with the notable exception of the simulation with a Ra = 10 8 and a radius ratio of 0.2.This case, showing an overestimation of approximately 40% in the case of the 2D scaled cylinder, becomes overestimated by the annulus only by ∼20% in the case of the 2D spherical annulus.A closer look at the simulations shows that the internally heated simulations in 2D tend to require more convective strength in order to transport the same amount of heat compared to a 3D simulation, as already shown by Hernlund and Tackley (2008).This is clearly seen in the comparison of v rms at Ra = 10 4 , where we observe that the steady-state root mean square velocity in 2D geometries is much greater than in the 3D spherical shell geometry, leading to an overestimation of 99%.
In comparison to the 2D scaled cylinder, the 2D spherical annulus gives less conclusive results for the top temperature gradient, as it does not provide better results.We can see an improvement in case of small radius ratios (i.e., f = 0.2) with either purely internal heating or mixed heating, or in the case of high Rayleigh numbers (i.e., Ra ≥ 10 7 ) for purely internally heated cases.
Otherwise, for the majority of the setups and heating modes, the 2D spherical annulus does not show significant improvement in reproducing the 3D spherical shell values compared to the scaled cylindrical geometry.

Stagnant Lid Convection
In the second part of the steady-state simulations, we modeled stagnant lid convection while varying the heating mode and the radius ratio in order to study how accurately the 2D geometries can reproduce the results obtained with a 3D spherical shell scenario (see Table S2 in Supporting Information S1).While we compared the 2D scaled cylinder and the 2D spherical annulus with the 3D spherical shell in the isoviscous cases, we included simulations in a 2D non-scaled cylinder geometry for the temperaturedependent models to test the effects of this scaling on the results.In this type of setup, where the viscosity is strongly temperature-dependent, convection is expected to operate in a stagnant lid regime.This rigid layer that forms at the top of the domain restricts convective heat transfer to the deep interior.
The results are summarized in Figure 3, where we report the relative error of the 2D cylindrical, 2D scaled cylindrical, and 2D spherical annulus geometry compared to the 3D spherical shell.We inspect the mean temperature, the root mean square velocity, the top temperature gradient, and the stagnant lid thickness.For each simulation we use a Ra of 5 × 10 6 and Ra Q of 5 × 10 7 .We also show a more detailed comparison of these geometries for three specific setups, namely a bottom-heated case with a radius ratio of f = 0.2, a mixedheated case with a radius ratio of f = 0.2 and an internally heated case with a radius ratio of f = 0.5.The case of a purely bottom-heated and small radius ratio setup was selected to better investigate the large difference of temperature between 2D and 3D as reported by Guerrero et al. (2018).The purely internally heated case with a medium radius ratio enables us to better distinguish the difference between downwellings in 3D and 2D geometries.Lastly, the mixed-heated case with a small radius ratio setup shows that the hotter temperature anomaly occurring in the purely bottom-heated 2D geometries tends to be mitigated by the addition of internal heating.We illustrate the differences in the temperature distribution between the geometries in these three cases in Figures 4-6, while in Figure 7 we compare the profiles of temperature, radial velocity and root mean square velocity among the geometries.
The purely bottom-heated setup in Figure 3 indicates that the 2D spherical annulus provides only a marginally better approximation of the temperature in the 3D spherical shell compared to the scaled and non-scaled cylinder, while not showing any improvement in the other diagnostic values.The highest discrepancy is observed for the radius ratio f = 0.2, with >50% of overestimation on average for the 2D geometries.For the 2D spherical annulus, the mean temperature matches best the one of the 3D spherical shell when the radius ratio increases, an effect also seen for the v rms .The top temperature gradient does not follow this pattern, however, while the error in the stagnant lid thickness is the lowest for f = 0.5.The largest difference to the 3D spherical shell for the low radius Averaged relative error to the 3D geometry with three different heating modes, that is, purely bottom-heated (top row), purely internally heated (middle row), and mixed-heated (bottom row) for steady-state isoviscous convection simulations.For each heating mode the first row shows the relative errors for the 2D spherical annulus and the second row represents the 2D scaled cylinder.Blue colors indicate an underestimation of the results obtained in the 3D geometry whereas a red color indicates an overestimation.The thermal and internal Rayleigh numbers vary from 10 4 to 10 8 and the radius ratio from 0.2 to 0.8.The radius ratio indicated on the plot is the one corresponding to reference 3D spherical shell simulations.In the case of the 2D cylinder, the radius ratio is scaled, thus a 0.2 radius ratio becomes 0.04 following the scaling formula from van Keken (2001) (Equation 12).The columns from left to right, show the relative error of the mean temperature, of the v rms , and of the top temperature gradient.ratio setups ( f = 0.2) in the purely bottom-heated setup that we found in our models is in agreement with Guerrero et al. (2018).The study by Guerrero et al. (2018) considered stagnant lid convection with only basal heating, and observed that for small radius ratio setups, the 2D spherical annulus would systematically show a hotter temperature when compared to the 3D spherical shell.This behavior is indeed confirmed by our simulations, where we can observe that the different geometry of plumes between the 2D spherical annulus and the 3D spherical shell create considerable differences in temperature for the low radius ratio, bottom-heated simulations.This is shown in Figures 4 and 7, where the 3D model presents the highest maximum temperature as well as a higher excess temperature (i.e., the residual temperature in Figure 4).On the other hand, its average domain temperature and its root mean square velocity are significantly lower than in the 2D geometries (almost by one order of magnitude in the case of the v rms ).The 3D spherical shell cases depict a maximum radial velocity larger by one order of magnitude compared to the 2D simulations.This clearly shows that the upwellings in the case of the 3D spherical shell are hotter and faster than their counterparts in 2D, and could be explained by the fact that plumes in a 2D Figure 3. Average relative error to the 3D geometry with three different heating modes (i.e., purely bottom-heated, purely internally heated, and mixed-heated) for steady-state convection simulations with a temperature-dependent viscosity.Blue colors indicate an underestimation of the results obtained in the 3D geometry whereas a red color indicates an overestimation.We compare the mean temperature of the domain, the root mean square velocity, the temperature gradient at the top of the domain, and the stagnant lid thickness.Each subplot shows the relative error in a given heating mode.Every column represents a different geometry, namely the 2D nonscaled cylinder, the 2D scaled cylinder, and the 2D spherical annulus.The thermal Rayleigh number is 5 × 10 6 , the internal Rayleigh number is 5 × 10 7 , and the radius ratios are 0.2, 0.5, and 0.8.geometry are rather sheet-like than column-like as they are in a 3D geometry, leading in turn to a larger mean temperature in the 2D geometries compared to the 3D spherical shell.As explained by Guerrero et al. (2018), this could be the consequence of a topological difference between flows in a 3D spherical shell and a 2D spherical annulus.The difference in the number of degree of freedom between 2D and 3D geometries leads to more sheetlike plumes in 2D which are colder and slower, but spreading in a larger part of the domain, compared to plumes in 3D which are faster, hotter, more localized, and columnar.Moreover, this difference is strongly accentuated for low radius ratios, where the number of plumes is small and such geometrical effects significantly affect the results.
For purely internally heated scenarios, we observe that the annulus fares remarkably better than the cylindrical geometries for a low radius ratio.In contrast to the bottom-heated cases, the average relative error of the mean temperature and the v rms increase as the radius ratio increases.The highest radius ratio (i.e., f = 0.8) shows the highest discrepancy between 2D and 3D geometries, with more than 50% and 20% of overestimation for the v rms and the top temperature gradient, respectively, while the stagnant lid thickness is underestimated on average by 15%.The snapshots in Figure 5 show that the flow is largely dominated by the thermal instabilities in the cold thermal boundary layer.We observe a larger number of downwellings in the case of a 3D spherical shell model for  d-i) representing the temperature distribution at different depths in a temperature-dependent viscosity case with purely basal heating for a radius ratio of f = 0.2, with Ra = 5 × 10 6 .First row: Temperature snapshots of three simulations using a 3D spherical shell (left), a 2D spherical annulus (middle), and a 2D scaled cylinder (right).An isotherm is also displayed in the case of the 3D spherical shell at T = 0.6 to show the distribution of upwellings.The white circles indicate the non-dimensional depth of 0.1, 0.5, and 0.9 used for the histograms in the second and third rows.Second row: Histograms showing the temperature distribution at different depths, namely at a non-dimensional depth of 0.1 (d), 0.5 (e), and 0.9 (f).Third row shows histograms of temperature variations relative to the average temperature at a given non-dimensional depth of 0.1 (g), 0.5 (h), and 0.9 (i).The histograms show the number of points at a given depth, which attain a specific temperature value, as percentage of the total number of points at that depth.The bin size is 0.01 in all the panels.a given internal heating rate when compared to the 2D geometries, which leads to a lower average domain temperature and root mean square velocity (Figure 7).Moreover, the minimum radial velocity and the minimum domain temperature (see Figures 7c and 7f) indicate colder and slower downwellings in 3D compared to 2D.
For cases heated both from below and from within (i.e., mixed-heated cases), the 2D spherical annulus shows the best match to the 3D spherical shell results for the mean temperature, the top temperature gradient, and the stagnant lid thickness when the radius ratio is the lowest (i.e., f = 0.2).For the highest radius ratio considered here (i.e., f = 0.8), the results obtained with the 2D spherical annulus and the 2D cylindrical geometries (both scaled and non-scaled) show similar errors, with the exception of v rms that seems to be best reproduced by the 2D scaled cylinder geometry.
When combining basal heating with internal heating, the discrepancy previously arising from the low radius ratio and purely basal heated simulations seems to be mitigated by the addition of internal heating.The difference in the temperature distribution between the 2D spherical annulus and the 3D spherical shell tends to disappear.Figures 6,7b,and 7e show that the difficulty of the 2D spherical annulus to reproduce the 3D spherical shell temperature distribution for low radius-ratio is reduced when combining basal and internal heating, and the difference in temperature and velocity of the hot plumes between the 3D and the 2D geometries is reduced compared to the purely basal heated setup.These 3D plumes will be on average faster, but not hotter than the plumes observed in 2D.
In summary, for the low radius-ratio stagnant lid cases presented here, the mean temperature is overall better reproduced by the 2D spherical annulus, showing a net improvement in the purely internally heated and mixed- heated cases, while only exhibiting a marginal amelioration in the purely bottom-heated cases.Additionally, in the internally heated and intermediate to low radius ratio simulations ( f ≤ 0.5), the 2D spherical annulus can also reproduce the 3D velocities, the top temperature gradient, and the stagnant lid thickness better than the 2D cylindrical geometry.It is also interesting to note that in the majority of the cases the 2D scaled cylinder geometry gives closer results to the 3D spherical shell when compared to the non-scaled cylinder.This confirms the results of van Keken (2001), indicating that the scaling of the core and mantle surface areas should always be used, when employing a cylindrical geometry.

Thermal Evolution Models
In the following section we compare the 2D and 3D geometries in a more complex setup by modeling the thermal evolution in a stagnant lid regime.Compared to the previously discussed steady-state calculations, these models illustrate which of the 2D geometries (2D cylinder, 2D scaled cylinder, and 2D spherical annulus) can best reproduce the 3D spherical shell results when mantle and core cooling are considered.

Mathematical Model
For the thermal evolution models we use the Extended Boussinesq Approximation (Schubert et al., 2001) to account for adiabatic heating and cooling.The energy equation (Equation 3) becomes:

DT
where u r is the radial component of the velocity vector, T surf is the surface temperature, α is the thermal expansivity, and Φ is the viscous dissipation given by Φ = τ : ε/2, where τ is the deviatoric stress tensor and ε the strain rate tensor.The dissipation number Di is defined as follows: where c p is the mantle specific heat capacity.
In the thermal evolution models, we consider a temperature and pressure dependent viscosity that follows the Arrhenius law for diffusion creep.The non-dimensional equation for the viscosity reads (Roberts & Zhong, 2006): where E and V are the activation energy and the activation volume, respectively (Hirth & Kohlstedt, 2003;Karato & Wu, 1993).T ref and z ref are the reference temperature and depth, respectively, at which the reference viscosity is attained.The non-dimensional T ref and z ref values correspond to a dimensional reference temperature of 1,600 K and a dimensional reference pressure of 3 GPa, respectively.
The temperature of the lower boundary T CMB evolves following a 1-D energy balance, assuming a core with constant density and heat capacity (Stevenson et al., 1983): where c c is the specific heat capacity of the core, ρ c is the core density, V c is the core volume, q c is the heat flux at the CMB, and A c is the core surface area.Here we do not consider core crystallization.
As appropriate for thermal evolution models, we take into account the decay of heat producing elements (i.e., 238 U, 235 U, 232 Th, and 40 K).The amount of radioactive heat sources is calculated using the concentrations listed in Table 1.
In these models we consider the effect of a 50 km laterally homogeneous crust.This crust is enriched in radiogenic elements while the mantle is depleted according to the following mass balance: where M is the mass and Q is the heating rate.This setup also considers the blanketing effect of the crust by using a lower thermal conductivity k in the crust compared to the mantle.

Case Definition
We test three scenarios for the thermal evolution of a Mars-like, a Moon-like, and a Mercury-like planet.We model the entire evolution of these planets to determine the variations over time of several key output quantities.These three planetary bodies were chosen because of their different interior structures.In the Mars-like case, the ratio between core and planetary radius is f = 0.544.The Moon-and Mercury-like scenarios represent two endmembers in terms of their radius ratios, with f = 0.224 and f = 0.828, respectively.
In the case of the 2D grids, we use a radial resolution of ∼10 km for Mars-and Moon-like cases, and a ∼5 km radial resolution for the thin mantle of Mercury.In the case of the 3D grids, the radial resolution lies between 22 and 9 km, while the lateral resolution is between 59 and 30 km at the surface (40,962 points per shell).More details about the resolution of each grid are available in Table S2 in Supporting Information S1.
In the 2D cylindrical geometry, we use both the classical 2D cylinder and the scaling of van Keken ( 2001), as done previously in the stagnant lid steady-state simulations.It is worth to note that for the peculiar Moon-like interior structure, the scaling leads to a radius ratio of f = 0.0503 (compared to a radius ration of f = 0.2241 in the nonscaled geometry).Since this is an extreme case, we test whether the scaling of the 2D cylinder geometry is appropriate for such interior structures in a thermal evolution scenario.
It is important to note that for each planet we use different Rayleigh numbers (i.e., Ra and Ra Q ).These are calculated self-consistently using the mantle thickness, internal heat sources, and temperature difference across the mantle for each planet.The parameters are listed in Table 1.While for Mars-and Moon-like simulations we use a reference viscosity of 10 21 Pa s (Laneuville et al., 2013;Plesa et al., 2015), for the Mercury-like case we perform additional tests with a reference viscosity that is lower by two orders of magnitude (i.e., 10 19 Pa s).Since for a reference viscosity of 10 21 Pa s the thin Mercurian mantle typically falls into a conductive state after a few Gyr of evolution (Tosi et al., 2015), by using a lower reference viscosity, we test additional scenarios, in which convection can be sustained over most of the evolution.
The 50 km primordial crust, which was used in all thermal evolution scenarios, is enriched by a factor of 2 in radiogenic elements compared to the bulk heat sources of the planet and has a two times lower thermal conductivity than the mantle (see Table 1).While the crustal enrichment factor used here is representative for the Mercurian crust (Tosi, Grott, et al., 2013), for the Mars-and Moon-like scenarios it allows to increase the complexity of the simulation without having a thermal evolution dominated by the crustal heat sources.A second set of thermal evolution models neglecting the effects of a primordial crust is listed in Section S7 in Supporting Information S1.

Results
In the following, we present the results obtained for the thermal history of Mars, the Moon, and Mercury for all geometries.Similar to the steady-state calculations, we show in Figure 8 the difference of the 2D non-scaled cylinder, 2D scaled cylinder, and 2D spherical annulus to the 3D spherical shell geometry after 4.5 Gyr of evolution (i.e., at present day).In addition to the error shown in percent, we also list the difference between the dimensional values for the mean temperature, CMB temperature, root mean square velocity, stagnant lid thickness, as well as surface and CMB heat fluxes (Figure 8).We examine the results presented in Figure 8 simultaneously with the temporal trends illustrated in Figures 9-11 for each planetary body.
Figure 9 shows the evolution of the output quantities of interest for all four geometries (2D non-scaled cylinder, 2D scaled cylinder, 2D spherical annulus, and 3D spherical shell) for a Mars-like structure.Our results show clearly that in the case of a Mars-like structure, the values obtained during the entire thermal evolution with a 3D geometry are more closely reproduced by the 2D spherical annulus than by the 2D cylinder, irrespective of whether the latter is scaled or not.It is to note that the 2D spherical annulus geometry is reproducing especially well (with an error smaller than 2%) the mean and CMB temperatures (see Figure 8 for the actual present-day error), while the velocities are systematically overestimated by the 2D geometries (27.7% for the 2D spherical annulus and 29.6% for the 2D scaled cylinder).The velocity discrepancy between the 2D and 3D geometries directly affects the calculation of the stagnant lid thickness.The stagnant lid thickness obtained in the 3D spherical shell is underestimated by 6.2% for the 2D spherical annulus and 8.1% for the 2D scaled cylinder.When trying to reproduce the heat fluxes obtained by the 3D spherical shell model at present day, the 2D spherical annulus is somewhat better than the 2D scaled cylinder, with an approximated underestimation of the CMB heat flux by 22% compared to 28% for the 2D scaled cylinder, while the surface heat flux will be overestimated by 6% and 8% for the 2D spherical annulus and the 2D scaled cylinder, respectively.While we focus in this part mostly on present-day values, when examining the entire thermal evolution of the planet we observe in the early and middle stages (between 1 and 3 Gyr) a relative error even larger than the one observed at present day, as seen on Figures 9c and 9e.The Mars-like setup is the most challenging setup to reproduce for the 2D spherical annulus, and although exhibiting relatively low errors, it still shows the highest discrepancies among the thermal evolution cases with different interior structures.The overestimation of the mean temperature, mean velocity, surface heat flow, and the underestimation of the stagnant lid thickness can be linked directly to the results obtained for the steady-state mixed-heated, temperature-dependent viscosity case with a radius ratio of 0.5, which shows the same type of relative error (compare Figure 3).
In the case of a Moon-like setup (Figure 10), the difference between the 2D spherical annulus and the 3D spherical shell is significantly lower than for the Mars-like case in terms of the mean temperature, the stagnant lid thickness and the surface heat flow.However, the differences between the 2D spherical annulus and the 2D cylinder (whether scaled or not) are typically larger than in the Mars-like case; as can be seen on Figures 10a-10c and 10f.The temperature evolution is very well reproduced by the 2D spherical annulus (less than 2% of error compared to more than 12% for the cylindrical geometries).Similarly, the surface heat flux and the stagnant-lid thickness show a good match between the 2D spherical annulus and the 3D geometry (see Table S6 in Supporting Information S1).Concerning the cylindrical geometries, the effects of the scaling are plainly visible on the overall temperatures and heat fluxes.As van Keken (2001) showed, the 2D non-scaled cylinder tends to overestimate the relative importance of the CMB radius compared to planetary radius and requires a scaling of the radii.However, even when the scaling is applied the results are still largely different compared to a 3D spherical shell geometry.A better approximation of the 3D results is obtained by the 2D spherical annulus, where such scaling is not needed.We see that in the non-scaled cylindrical geometry the CMB heat flux stays at around 1.86 mW m 2 even at present day, meaning that the core is actively heated by the mantle, although in the other geometries the core is already cooling (see Figures 10b and 10d).This behavior was also seen for Mars with the non-scaled cylindrical geometry.Similar to what has been seen previously for the steady-state low aspect-ratio cases with temperaturedependent viscosity and mixed heating (see Figure 3), the scaled and non-scaled cylindrical geometries show large disagreements with 3D results in all studied metrics.Nevertheless, even for the 2D spherical annulus, the mean velocity and the CMB heat flux present the largest errors among the investigated quantities.
In Figure 11, we show the results of the thermal evolution for a Mercury-like interior structure.Here we used two sets of simulations in order to illustrate the case of a weakly convecting mantle with a reference viscosity of η ref = 10 21 Pa s resulting in a Rayleigh number of Ra = 3.49 × 10 4 and the case of a more vigorously convecting mantle with a reference viscosity lowered by two orders of magnitude, thus increasing the Rayleigh number to Ra = 3.49 × 10 6 .Here, too, the global trend that has already been observed for Mars and the Moon is evident.As shown in Figure 8, the 2D spherical annulus geometry reproduces the 3D spherical shell results with an approximate error of less than 1%, with a notable exception for the velocities, which are highly overestimated (more than 90% of relative error).The very high relative error of the v rms is explained by the present-day state of the Mercurian mantle.In our simulations, a Mercury-like planet falls into a quasi-conductive state after a couple of Gyr of evolution irrespective of the geometry (Figure 11), which in turn gives very low absolute v rms values.
Yet the absolute difference of velocity between the geometries is very small (less than 1 × 10 3 cm/year).Despite this high relative error of the velocity, the stagnant lid thickness is, however, quite well reproduced by the 2D geometries, giving a maximum relative error of 19 km (or an underestimation of the 3D spherical shell results by 6.5%).

Melting and Mechanical Thickness
In addition to the previously analyzed variables such as mantle temperature, surface heat flux and v rms , Figure 12 displays additional post-processing parameters, that is, the degree of partial melting and the mechanical lithosphere thickness.The latter have been calculated in order to better characterize the differences in thermal evolution between the 3D spherical shell, the 2D spherical annulus, and the 2D scaled cylinder geometries.The fraction of molten mantle and the thickness of the mechanical lithosphere as a function of time are displayed along with temperature profiles at different times during the evolution.For the calculation of the mechanical lithosphere thickness, we follow the approach of Grott and Breuer (2008) and use a dry rheology for the crust and the mantle as a first approximation.The calculations involving partial melting of the mantle are highly simplified and do not include the effects of latent heat or mantle depletion.We compute a volumetrically averaged degree of melting by comparing the local mantle temperature with the solidus and liquidus of Takahashi (1990) and assuming a linear melt relation between solidus and liquidus.We then compare the volume of melted mantle to the total volume of the mantle.While the quantities presented in Figure 12 are based on simple post-processing of the thermal evolution results, they are meant to provide first order implications for the thermal evolution modeling with 2D and 3D geometries (for additional information concerning the post processing, see Sections S8 and S9 in Supporting Information S1).
When reproducing 3D simulations, the 2D spherical annulus is well suited to compute melting and mechanical lithosphere thickness.The 2D scaled cylinder shows a stronger deviation with respect to the 3D results, especially for the Moon and Mercury cases.The mechanical thickness for a Moon-like interior structure will be underestimated on average by 5% for the 2D spherical annulus and 20% by the 2D scaled cylinder, while the amount of melting will be overestimated by 10% at present day in the case of the 2D spherical annulus and by an average of 35% by the 2D scaled cylinder.In our simulations the temperature profiles of the Moon are predominantly higher than the solidus and thus show a high degree of partial melting.While these profiles stay fairly similar until 1-1.5 Gyr (Figure 12d), the 2D scaled cylinder will show early on a large overestimation of the partial melt production compared to the 3D spherical shell geometry.For the Mercury-like setup, the 3D mechanical thickness is well reproduced by the 2D geometries, with the 2D spherical annulus having close to 0% error at present day and the 2D scaled cylinder having an approximate underestimation of 5%.Our simulations show a rather high degree of partial melting in the early phases of Mercury's evolution, with a peak at approximately 0.3-0.4Gyr (Figure 12c), which is well reproduced by the 2D spherical annulus with an overestimation of 7% while the 2D scaled cylinder fares worse with an overestimation of around 20%.In the case of a Mars-like structure (Figures 12b and 12e) the differences between 2D and 3D are larger, in particular for partial melting, which the 2D spherical annulus geometry will overestimate by a maximum of 30% while the 2D scaled cylinder by more than 35% at around 1.5 Gyr.The mechanical lithosphere thickness will be underestimated on average by 10% in the 2D spherical annulus case and by 14% in the 2D scaled cylinder.As seen previously in Figure 8, the 2D spherical annulus has the most difficulty in reproducing the results for a 3D Mars-like setup.Nevertheless, the approximation of the 3D diagnostics is still better for the 2D spherical annulus than for the 2D scaled cylindrical geometry.

Discussion
Our results confirm that the 2D spherical annulus can reproduce the 3D spherical shell geometry better than the 2D cylindrical geometry, consistent with the previous study of Hernlund and Tackley (2008).Using simulations of increasing complexity, our systematic study shows for the first time in great detail the difference in applying a 2D geometry instead of a more realistic 3D spherical shell domain when modeling thermal convection in planetary mantles.
With the 2D cylindrical geometry, whether scaled or not, the results differ substantially from the 3D spherical geometry results in both steady-state and thermal evolution scenarios.The necessity of choosing between a scaled and a non-scaled 2D cylinder in modeling geodynamic processes inevitably results in a trade-off between a more accurate representation of the deep interior structures in the case of the non-scaled 2D cylinder, which is especially important when studying thermochemical features (Kameyama, 2022;Nakagawa & Tackley, 2004;Stegman et al., 2003;Yu et al., 2019), and a correct representation of the heat fluxes as well as root mean square velocity of the mantle (Deschamps et al., 2010;Mulyukova et al., 2015) in case of the scaled 2D cylinder.To circumvent these inaccuracies, the systematic use of the 2D spherical annulus geometry in reproducing thermochemical convection in 3D spherical shell is strongly recommended.
In the case of steady-state simulations, we demonstrated that the 2D spherical annulus has the largest error in the high radius ratio scenarios (i.e., f = 0.6 and 0.8).Thus, the efficiency of the 2D spherical annulus in reducing the error to the 3D spherical shell (compared to the results of the 2D scaled cylinder) is most visible in the case of a low radius ratio configuration (i.e., f = 0.2).In general, large discrepancies between 2D and 3D geometries occur for the mean temperature in the case of bottom-heated and temperature-dependent viscosity setups, as also seen by Guerrero et al. (2018).
The purely internally heated cases show the largest difference between the 2D and 3D geometries as also reported by Hernlund and Tackley (2008).Simulations using mixed heating (i.e., bottom and internal heating), on the other hand, show the smallest errors between the 2D spherical annulus and the 3D spherical shell, as the difference in the temperature distribution between these two geometries tends to decrease (see Figure 6).This fact becomes important when modeling the thermal evolution of terrestrial planets, since the silicate mantles of planets always exhibit heating caused both by the presence of radiogenic elements in the mantle and by core cooling.The smaller error between the 2D spherical annulus and the 3D spherical shell observed in mixed-heated cases makes the 2D spherical annulus an acceptable alternative to model more complex scenarios (Figure 8), for which a 3D geometry can be too expensive.The trend of the relative error in the steady-state stagnant lid simulations with mixed heating is also observed in the case of thermal evolution models: that is, the errors in the surface heat flux and stagnant lid thickness decrease with decreasing radius ratio (cf. the errors obtained for the Moon and Mars in Figure 8).However in the case of Mercury, while the steady-state simulations would predict the largest errors, the low Rayleigh number for the Mercurian mantle and the transition to a conductive state during the thermal evolution strongly reduce the discrepancy between 2D and 3D geometries.
The results also indicate that 2D geometries generally tend to underestimate the mechanical thickness of the lithosphere, while overestimating the amount of melting during the evolution.Again, the 2D spherical annulus gives systematically better results and is better suited to replicate melting and mechanical thickness obtained for a 3D spherical shell geometry compared to the 2D scaled cylinder geometry.The underestimation of the mechanical lithosphere thickness by the 2D geometries is explained primarily by warmer mantle temperatures.The hotter temperature profiles (as seen on Figure 12) in the case of the 2D spherical annulus and especially in that of the 2D scaled cylinder, reduces the thickness of the stagnant lid and thus the depth of transition between brittle and ductile deformation.The 2D spherical annulus also approximates better what is used as a proxy to represent the amount of melt in the thermal evolution cases (Figures 12d-12f).Since the mantle temperature obtained in a 3D spherical shell geometry is better approximated by the 2D spherical annulus than by the 2D cylindrical geometry, also the melting history is more accurately reproduced by the 2D spherical annulus.Moreover, as shown in Figure 1, the 2D spherical annulus grid is constructed in a way that the total sum of all volumes matches the one of a 3D spherical shell.Therefore, the calculation of a melt volume in the annulus is much closer to the 3D spherical shell than in the 2D scaled cylinder geometry.The largest error reduction occurs for the 2D spherical annulus in the Moon-like case, showing again that the improvement brought by the 2D spherical annulus geometry will be the highest in a small radius-ratio setup, whereas for an intermediate radius-ratio setup (Mars-like), the annulus shows the least error reduction from 2D scaled cylinder to 2D spherical annulus.
Although the 2D spherical annulus geometry can better reproduce the 3D spherical shell results compared to the 2D cylindrical geometry, care must be taken especially when modeling partial melting.The higher temperature and thus the higher amount of partial melting could lead to an overestimation of the crustal thickness in the 2D geometries as compared to the 3D spherical shell geometry.This would have further implications for the thermochemical evolution.A thicker crust may lead to higher interior temperatures and thus crustal production rate due to a pronounced "crustal blanketing" (e.g., Schumacher & Breuer, 2006), in which the reduced thermal conductivity of the crust will prevent efficient cooling of the mantle.However, a higher crustal production rate also leads to a stronger depletion of the mantle in radioactive elements promoting a stronger cooling of the interior.Moreover, partial melting processes can strongly affect interior outgassing and associated mantle depletion in volatile elements.Which effect dominates depends on the values of the thermal conductivity in the crust and the enrichment factor of the radioactive element from the mantle into the crust.In any case, care should be taken when 2D geometries (both spherical annulus and cylinder) are employed to study partial melting and subsequent crust production or degassing.This is particularly true for planets with an intermediate radius ratio, like Mars or Venus, for which the error in the amount of partial melt is the largest.These processes and the differences between 2D and 3D geometries for such scenarios need to be quantified in future studies.

Conclusions
We provided a systematic comparison between different geometries in order to determine how accurately 2D geometries reproduce 3D results.To do so, we investigated (scaled and non-scaled) 2D cylinder, 2D spherical annulus, and 3D spherical shell geometries in a series of scenarios.We started with isoviscous steady-state models, included the effects of a temperature-dependent viscosity, and finally tested the different geometries for thermal evolution scenarios.Our main findings are the following: 1.While it is obvious that a 3D spherical shell geometry is preferable to a 2D geometry, this is not always feasible due to the high computational cost.When applying models of varying complexities, we demonstrated that the 2D spherical annulus geometry is able to reproduce the 3D spherical shell models much better than the 2D cylinder, especially for the low radius ratio setups.The latter is also clearly seen when modeling the thermal evolution of the Moon.2. For steady-state scenarios, our models show that the 2D geometries will mostly overestimate the mean temperature compared to the 3D spherical shell, a result largely explained by the geometry of mantle plumes (i.e., sheet-like in 2D vs. columnar-like in 3D).This discrepancy decreases with an increasing Rayleigh number but is more accentuated for low-radius ratio cases, a result already observed by Guerrero et al. (2018).
The difference in temperature between the 2D and 3D geometries decreases for mixed-heated cases (i.e., heated both from below and from within).This is especially true in the case of the 2D spherical annulus, since it uses the same cell volumes as a 3D spherical shell.3.For thermal evolution models, we found that for medium ratios between inner and outer radius (e.g., Mars-like structure), the differences in the results for the 2D and 3D geometries are larger than for extreme radius ratios.
In contrast to the temperature-dependent steady-state cases, where the difference in surface heat flux and stagnant lid thickness between 2D geometries and 3D geometries is largest for high radius ratios, the difference obtained for Mercury-like evolution parameters is minimal.This is due to the low Rayleigh number of Mercury, which leads to the transition to a conductive state during its thermal history.4. When investigating melting processes with 2D geometries in thermal evolution setups, the results for small (Moon) or large (Mercury) radius ratios can be reproduced well with 2D spherical annulus.However, caution is required for medium radius ratios (e.g., Mars and Venus): Although the 2D spherical annulus approximates the results of a 3D spherical shell geometry better than the 2D cylinder, even in this case, crustal production could be overestimated by up to 30% compared to a 3D spherical shell simulation, which may lead to a different thermal history of the interior.This should be taken into account when interpreting the results of the models.
Future studies need to test the accuracy of the 2D spherical annulus in reproducing the 3D spherical shell geometry in even more complex scenarios considering variable thermal conductivity and expansivity (Tosi, Yuen, et al., 2013), chemical buoyancy (Nakagawa et al., 2010), as well as partial melting of the mantle and its influence on thermal evolution.

Figure 1 .
Figure 1.(left) 3D representation of a single cell in the 2D spherical annulus geometry.(middle) The spherical annulus geometry with six radial shells.The filled cells show the individual cells of the grid.Since the volumes and areas are calculated according to this shape, the sum of all volumes matches the volume of a 3D spherical shell.This shape creates accurate heat conduction profiles of a 3D spherical shell while maintaining the typical 2D connectivity.(right) Cylindrical projection with an arbitrary height of the same grid.

Figure 2 .
Figure2.Averaged relative error to the 3D geometry with three different heating modes, that is, purely bottom-heated (top row), purely internally heated (middle row), and mixed-heated (bottom row) for steady-state isoviscous convection simulations.For each heating mode the first row shows the relative errors for the 2D spherical annulus and the second row represents the 2D scaled cylinder.Blue colors indicate an underestimation of the results obtained in the 3D geometry whereas a red color indicates an overestimation.The thermal and internal Rayleigh numbers vary from 10 4 to 10 8 and the radius ratio from 0.2 to 0.8.The radius ratio indicated on the plot is the one corresponding to reference 3D spherical shell simulations.In the case of the 2D cylinder, the radius ratio is scaled, thus a 0.2 radius ratio becomes 0.04 following the scaling formula from van Keken (2001) (Equation12).The columns from left to right, show the relative error of the mean temperature, of the v rms , and of the top temperature gradient.

Figure 4 .
Figure 4. Temperature snapshots (a-c) and histograms (d-i) representing the temperature distribution at different depths in a temperature-dependent viscosity case with purely basal heating for a radius ratio of f = 0.2, with Ra = 5 × 10 6 .First row: Temperature snapshots of three simulations using a 3D spherical shell (left), a 2D spherical annulus (middle), and a 2D scaled cylinder (right).An isotherm is also displayed in the case of the 3D spherical shell at T = 0.6 to show the distribution of upwellings.The white circles indicate the non-dimensional depth of 0.1, 0.5, and 0.9 used for the histograms in the second and third rows.Second row: Histograms showing the temperature distribution at different depths, namely at a non-dimensional depth of 0.1 (d), 0.5 (e), and 0.9 (f).Third row shows histograms of temperature variations relative to the average temperature at a given non-dimensional depth of 0.1 (g), 0.5 (h), and 0.9 (i).The histograms show the number of points at a given depth, which attain a specific temperature value, as percentage of the total number of points at that depth.The bin size is 0.01 in all the panels.

Figure 5 .
Figure 5. Same as Figure 4 but for a temperature-dependent viscosity case in a purely internally heated scenario with a radius ratio of f = 0.5 and Ra Q = 5 × 10 7 .

Figure 6 .
Figure 6.Same as Figures 4 and 5 but for a temperature-dependent viscosity case in a mixed-heated scenario with a radius ratio of f = 0.2, Ra = 5 × 10 6 , and Ra Q = 5 × 10 7 .

Figure 7 .
Figure 7. Profile comparison between the 3D spherical geometry (solid lines), the 2D annulus geometry (dashed lines) and the 2D scaled cylinder (dotted lines).First row (panels a-c): Temperature profiles.Second row (panels d-f): Radial velocity profiles.Third row (panels g-i): Root mean square velocity profiles.Each column corresponds to one of the three specific cases shown in Figures 4-6 (left to right), namely the purely bottom-heated case with f = 0.2, the mixed-heated case with f = 0.2, and the purely internally heated case with f = 0.5.The minimum and maximum temperature profiles show the minimum and maximum values attained at each depth, respectively.The positive values of the radial velocity indicate the velocity of mantle plumes, while the negative values indicate the velocity of cold downwellings (in panel d, the minimum radial velocity in the 3D spherical shell geometry is small, but non-zero).

Figure 8 .
Figure 8. Relative error to the 3D in the case of thermal evolution simulations.In each panel we vary the planet on the y axis and the geometry on the x axis.We investigate the relative error for: the mean temperature (a), the core mantle boundary (CMB) temperature (b), the root mean square velocity (c), the stagnant lid thickness (d), the surface heat flux (e), and the CMB heat flux (f).Blue colors indicate an underestimation of the results obtained in the 3D spherical shell geometry whereas red colors indicate an overestimation.2D Cyl stands for non-scaled 2D cylindrical geometry, 2D Sca Cyl is the 2D cylindrical geometry with the scaling by van Keken (2001), and 2D Ann stands for 2D spherical annulus geometry.Each panel shows the relative error in % and the absolute error in dimensional units, when compared to the 3D spherical shell geometry.

Figure 9 .
Figure 9.Time series of a Mars-like case with an initial crust of 50 km and with four different geometries (2D non-scaled cylindrical, 2D scaled cylindrical, 2D spherical annulus, and 3D spherical shell).The panels show the mean temperature (a), the core-mantle boundary (CMB) temperature (b), the surface heat flux (c), the CMB heat flux (d), the average root mean square velocity of the domain (e), and the stagnant lid thickness (f) from 4.5 Gyr ago to present day, respectively.The shaded areas show the min.-max.variations during the evolution.The 2D non-scaled cylinder has been added to show the effect of the scaling introduced by van Keken (2001).

Figure 10 .
Figure 10.Time series of a Moon-like case with an initial crust of 50 km and with four different geometries (2D non-scaled cylindrical, 2D cylindrical, 2D spherical annulus, and 3D spherical shell).The panels show the mean temperature (a), the core-mantle boundary (CMB) temperature (b), the surface heat flux (c), the CMB heat flux (d), the averaged root mean square velocity of the domain (e), and the lid thickness (f) from 4.5 Gyr ago to present day, respectively.The shaded areas show the min.-max.variations during the evolution.The 2D non-scaled cylinder has been added to show the effect of the scaling introduced by van Keken (2001).

Figure 11 .
Figure 11.Time series of a Mercury-like case with an initial crust of 50 km and with four different geometries (2D non-scaled cylindrical, 2D cylindrical, 2D spherical annulus, and 3D spherical shell) and two different reference viscosities.The panels show the mean temperature (a), the core-mantle boundary (CMB) temperature (b), the surface heat flux (c), the CMB heat flux (d), the average root mean square velocity of the domain (e), and the lid thickness (f) from 4.5 Gyr ago to present day, respectively.The dotted lines represent simulations with a reference viscosity of η ref = 10 19 Pa s, while the solid lines represent the cases with a reference viscosity of η ref = 10 21 Pa s.The maximum and minimum of the output quantities are not displayed here, since these variations are negligible.The 2D non-scaled cylinder has been added to show the effect of the scaling introduced by van Keken (2001).

Figure 12 .
Figure 12.Time series and temperature profiles of the three scenario investigated (from left to right; the Moon, Mars, and Mercury).The first row shows the mechanical thickness of the lithosphere (in km) and the second row shows the fraction of molten mantle (in %) as a function of time.Solid line represents the 3D spherical geometry, the dashed line corresponds to the 2D spherical annulus and the dotted line to the 2D scaled cylinder.The following panels show the solidus (orange line), as well as average and maximum temperature profiles (black and red lines, respectively) at different times as indicated in each panel and illustrated by vertical dotted red lines on panels (d)-(f).The solidus is calculated following the work from Takahashi (1990).

Table 1
Parameters for Thermal Evolution Calculations for Mars, Moon, and Mercury That the non-dimensional radii are scaled according to Equation 12 for the scaled cylinder geometry.