Computational Modeling of Volcanism on Earth-like Planets


The aim of this paper is to model the level of volcanic activity on Earth-like planets over time. We consider the level of volcanic activity on a planet to be a function of the planet’s thermodynamic state. A planet is considered to become volcanically inactive once it has reached a thermodynamic steady state. We first model the heat flow in the planet’s interior by the heat equation, which we reduce to a one-dimensional Laplace equation. We then calculate the temperature field of the planet by numerically solving a spherically symmetric boundary value problem. Finally, we relate discrete time-steps in our simulation to real-time via an empirically-informed mapping. Our results indicate that Earth will remain volcanically active for a total of 6.0\sim 6.0 billion years since its formation, while Earth-like planets of 0.50.5 and 2.02.0 Earth masses will be active for 5.5\sim 5.5 and 7.7\sim 7.7 billion years respectively. Our model incorrectly predicts that Earth-like planets above 2.02.0 Earth masses continually increase their internal temperatures if their conductivity and density profiles are assumed to be identical to that of Earth, which suggests some limitations of the model.

1 Introduction

Planetary volcanism brings together many different fields of science. It is well known that terrestrial planets have evidence of volcanism. However, there is no current model for the time-evolution of volcanic activity on Earth-like planets. To make accurate predictions of the Earth’s internal temperature field, we first require an understanding of Earth’s interior, much of which eludes direct measurement.

Conventionally, the Earth is thought to be composed of different layers; the major ones considered here are the inner core, outer core, lower mantle, transition zone, upper mantle, low-velocity zone (LVZ) and crust. Each layer possesses a density profile as given by the PREM model presented by Anderson and Dziewonski, and a conductivity profile given by Stacey and Anderson ((1; 5)). These results are summarized in Table 1 .

The physical nature of heat flow in a material is typically modeled by the Heat Equation,

Tt=(αT)\frac{\partial T}{\partial t}=\nabla\cdot\left(\alpha\nabla T\right) (1)

where TT is the temperature and α\alpha is the thermal diffusivity constant. In practice, it is difficult to measure this diffusivity and it is usually expressed in terms of the conductivity kk, specific heat capacity CpC_{p}, and the density ρ\rho so that

α=kCpρ\alpha=\frac{k}{C_{p}\rho} (2)

Using these equations, we model heat dissipation throughout the Earth’s interior. Our aim in this paper is to relate the time needed to reach steady-state equilibrium solutions of the heat equation to the volcanic activity of a planet.

Table 1: Major regions of Earth’s interior and corresponding density and thermal conductivity values. rr is the radius in Earth radii.
Region Radius (Earth Radii) Density (g \cdot cm-3{}^{-3} ) Conductivity(W/mK)
Inner Core 0.192 13.0885 - 8.8381r2r^{2} 80
Outer Core 0.546 12.5815 - 1.2638rr - 3.6426r2r^{2} - 5.5281r3r^{3} 10
Lower Mantle 0.895 7.9565 - 6.4761rr + 5.5283r2r^{2} -3.0807r3r^{3} 5
Transition Zone 0.937 11.2494 - 8.0298rr 5
Upper Mantle 0.965 7.1089 - 3.8045rr 5
LVZ 0.996 2.691 + 0.6924rr 5
Crust 1.000 2.75 4.5

2 Methods

In order to simplify the problem a number of assumptions must be made. Current theories on planet formation suggest that terrestrial planets begin as molten rock and are hence volcanic ((4)). From this, we assumed that every planet begins as a uniform sphere with a constant temperature equal to it’s current core temperature. We then fixed the innermost core, r=0r=0, to that temperature and the surface to the effective surface temperature of Earth (255255 K) as calculated via the Stefan-Boltzmann equation ((2)). Since the planets are assumed to be Earth-like, we imposed these same boundary conditions on all planets.

We then set the density profile of each planet in accordance with the commonly accepted PREM model ((1)). The thermal diffusivity was calculated from Equation 2 where a specific heat capacity was assumed to be similar throughout the Earth (Cp1C_{p}\sim 1). A summary is given in Table 1.

By treating planets purely as conductively heated bodies, we avoid the laborious computation of self-gravitational effects, material compression, radioactive heating, convective heat transfer and radiative heat flux. In our simulation, the net result of these effects manifest as fixed boundary temperatures at the core and surface of the planet.

Table 1 shows the density and conductivity of Earth. The density was modeled as a piecewise function, with regions corresponding to the major layers of Earth: the inner core, outer core, inner mantle, transition zone, outer mantle, low velocity zone, and crust. The density was verified by simulating Earth and summing the density function over each spherical shell to obtain a total mass of 1.01.0 Earth masses. This summation technique was used to compute total mass for each planet simulated. Conductivity was also modeled as a piecewise function, with linear regions based on values from the literature ((5)). The primary purpose of the conductivity function was to modulate the behaviour of the temperature field within the interior of the planet.

In order to determine the time-evolution of the temperature field within the simulated planet, we first considered the one dimensional heat equation:

Tt=r(α(r)Tr)\frac{\partial T}{\partial t}=\frac{\partial}{\partial r}\left(\alpha(r)\frac{% \partial T}{\partial r}\right) (3)

By setting the time dependence to zero, we transform the problem into a spherically symmetric Laplace equation where, in general, α\alpha may also depend on radius. Thus

r(α(r)Tr)=0\frac{\partial}{\partial r}\left(\alpha(r)\frac{\partial T}{\partial r}\right)=0 (4)

By applying a numerical method, we iteratively solve the Laplace equation within the interior of the simulated planet. Figure 1 shows a schematic of the setup for our numerical solution to the boundary value problem. The planet is composed of NN concentrically stacked shells. The innermost and outermost shells, S0S_{0} and SNS_{N}, are at fixed boundary temperatures. For each shell in the interior region, SiS_{i}, its new temperature is computed via Equation 5

Figure 1: Schematic diagram of iterative procedure used to model the time-evolution of solutions to the Laplace Equation.
T(Si)=wi-1T(Si-1)+wi+1T(Si+1)wi-1+wi+1T(S_{i})=\frac{w_{i-1}T(S_{i-1})+w_{i+1}T(S_{i+1})}{w_{i-1}+w_{i+1}} (5)

Where the weights, wiw_{i}, are given by

wi=4παiri2w_{i}=4\pi\alpha_{i}r_{i}^{2} (6)

Starting at the 11st shell (i=1i=1), the new temperature was recursively computed for each iith shell using a weighted average of the nearest neighbour shells. The weighted average for each neighbour shell was determined by multiplying the surface area of the neighbour shell by its thermal diffusivity.

Successive iterations of our method to solve the boundary value problem can be interpreted as the evolution of time, and so we re-introduce time dependence into the simulation via this route. In other words, there exists some bijection between real time and the number of iterations in our numerical solution (i.e. solution-time).

It is assumed in our simulation that all planets begin in a volcanic state, and become non-volcanic when a thermodynamic steady-state is reached. The equilibrium time, teqt_{eq}, at which a planet becomes non-volcanic, is found by measuring the relative root mean squared (RMS) change in the temperature field between one iteration and the next. When the relative RMS change was less than 1%1\% between the (i-1)(i-1)th and iith iteration, then teqt_{eq} was set to tit_{i}.

3 Results

The time-evolved temperature field for a planet simulated at 1.01.0 Earth masses is presented in Figure 2. Since Earth is still volcanically active, it can be inferred from Figure 2 that the current age of the Earth, 4.54.5 billion years, must correspond to a number of iterations less than teqt_{eq}. It is interesting to note that the simulation predicts an eventual cooling of the mantle and core, which is consistent with physical theory (as any planet will eventually radiate all its heat away).

Figure 3 shows a 3D plot of the temperature field, evolving through time. Transient behaviour at the beginning of the simulation shows that the outer regions of the initial uniform-temperature mass quickly cool, while the core material experiences much slower cooling. The temperature distribution gains a distinct piecewise appearance as it equilibrates; this is a result of the thermal conductivity profile that has been applied to the planet.

The Moon is a body much smaller than Earth, with a mass roughly one hundredth of Earth’s. As expected, its teqt_{eq} is considerably less than that of the Earth. Figure 4 demonstrates that the Moon cools rapidly in comparison to Earth, thus reaching a steady-state at a much earlier time.

With the knowledge that the Moon has no current volcanic activity, it can be concluded that the current age of the Moon, which is again roughly 4.54.5 billion years, must be greater than teqt_{eq} for the Moon.

Figure 2: Time-evolved temperature contour plot for Earth. Horizontal lines show boundaries between piecewise density regions, and vertical line is the equilibrium time.
Figure 3: A three-dimensional representation of the contour map.
Figure 4: Time-evolved temperature contour plot for the Moon. Horizontal lines show boundaries between piecewise density regions, and vertical line is the equilibrium time.

At present, it is unknown whether Mars has any active volcanoes ((6)). The simulated data for Mars is shown in Figure 5, and it can be seen that the teqt_{eq} falls between the values obtained for Earth and the Moon. If it is given that Mars is currently near its teqt_{eq} in real-time, then the simulated teqt_{eq} may be set to 4.54.5 billion years ((6)). Knowing that our simulations start at the formation of the planet, and that a linear time-scale is being used, we now infer a mapping between real-time and simulation-time:

Figure 5: Time-evolved temperature contour plot for Mars. Horizontal lines show boundaries between piecewise density regions, and vertical line is the equilibrium time.
tratio=6771 steps4.5 billion yearst_{ratio}=\frac{6771\text{ steps}}{4.5\text{ billion years}} (7)

Using this result, it is now possible to predict the time-evolved volcanic activity of a planet of arbitrary mass by simulating the evolution of it’s temperature field, and then mapping simulation-time to real-time. The results of this mapping are shown in Figure 6.

Figure 6: Time taken to reach thermodynamic steady-state as a function of planetary mass. Red line shows linear fit.

Figure 6 shows that our simulation correctly predicts that the mass of a planet is positively correlated with the time until thermal steady-state is reached, which corresponds to the end of volcanic activity on the planet. Our simulation also correctly predicts that, at the present real-time (4.54.5 billion years since planet formation), the Moon and Mercury are volcanically inactive while Venus and Earth are volcanically active. Using the mapping between real-time and simulation-time, the simulation predicts that Earth will continue its volcanic activity until 66 billion years of age.

Figure 7: Three-dimensional plot of a planet at 2.162.16 Earth masses showing a continually increasing temperature field.

When simulating planets beyond the threshold of 2 Earth masses, it was found that the temperature field does not reach steady-state. As shown in Figure 7, the core of the planet radiates heat at such a rate that the interior temperature of the planet increases without bound. We attribute this effect to our assumption that all planets would have a radial density and conductivity profile similar to that of Earth’s, scaled linearly to the radius of the simulated planet. Our results suggest that scaling Earth’s conductivity and density profiles to larger radii is too simplistic of a method for properly modeling larger planets. It is also worth noting that in the recent work by Raymond et al., simulations of Earth-like planets were performed only up to 2.6 Earth masses, suggesting that simulation of terrestrial planets above this mass is non-trivial ((4)). Furthermore, we suspect that a lack of convective effects in our model contributes to this inconsistency, as convection would facilitate the dissipation of heat. The model of Kite et al. takes these effects explicitly into account, yielding consistent results for volcanism on planets with masses up to 25 times that of the Earth ((3)). The advantage of our model, in comparison, lies in its simplicity.

4 Conclusion

In this study, the volcanic activity of Earth-like planets was studied by equating volcanic inactivity to thermal steady-state. The time-evolution of planetary temperature-fields was solved using a novel method, in which a Laplace equation was solved iteratively and a mapping between real-time and simulation-time was found by comparison to real-world data.

It was found that Earth is still in its non-steady-state phase, and will continue to exhibit volcanic activity until it reaches an age of 6 billion years. Additionally, it was found that planets with masses greater than 2.0 Earth masses having a density and conductivity profile that is linearly proportional to Earth’s density profile would exhibit non-terrestrial behaviour due to continual self-induced internal temperature increase, in the absence of convection.

Future work in this area could include the non-linear scaling of density and conductivity profiles and the inclusion of a simple model of convective effects to determine how large Earth-like planets could be simulated to exhibit terrestrial behaviour.


  • D. L. Anderson and A. M. Dziewonski (1981) The preliminary reference earth model. Phys. Earth Plan. Int 25, pp. 297–356. Cited by: 1, 2.
  • J. T. Houghton (1986) The physics of atmospheres. Cambridge University Press. Cited by: 2.
  • E. S. Kite, M. Manga and E. Gaidos (2009) Geodynamics and rate of volcanism on massive earth-like planets. Astrophys.J. 700, pp. 1732–1749. Cited by: 3.
  • S. N. Raymond, T. Quinn and J. I. Lunine (2006) High-resolution simulations of the final assembly of earth-like planets 2: water delivery and planetary habitability. Astrobiology Special Issue on M Stars 7. Cited by: 2, 3.
  • F. D. Stacey and O. L. Anderson (2001) Electrical and thermal conductivities of fe-ni- si alloy under core conditions. Physics of The Earth and Planetary Interiors 124, pp. 153–162. Cited by: 1, 2.
  • [6] Volcanoes still active on mars? new evidence for ongoing volcanism and water release. Note: \url­/releases/2001/11/011109075016.htm Cited by: 3.

Additional Assets

Showing 1 Reviews


This article and its reviews are distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and redistribution in any medium, provided that the original author and source are credited.