Abstract
The aim of this paper is to model the level of volcanic activity on Earthlike 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 onedimensional Laplace equation. We then calculate the temperature field of the planet by numerically solving a spherically symmetric boundary value problem. Finally, we relate discrete timesteps in our simulation to realtime via an empiricallyinformed mapping. Our results indicate that Earth will remain volcanically active for a total of $\sim 6.0$ billion years since its formation, while Earthlike planets of $0.5$ and $2.0$ Earth masses will be active for $\sim 5.5$ and $\sim 7.7$ billion years respectively. Our model incorrectly predicts that Earthlike planets above $2.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 timeevolution of volcanic activity on Earthlike 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, lowvelocity 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,
$\frac{\partial T}{\partial t}=\nabla\cdot\left(\alpha\nabla T\right)$  (1) 
where $T$ 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 $k$, specific heat capacity $C_{p}$, and the density $\rho$ so that
$\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 steadystate equilibrium solutions of the heat equation to the volcanic activity of a planet.
Region  Radius (Earth Radii)  Density (g $\cdot$ cm${}^{3}$ )  Conductivity(W/mK) 

Inner Core  0.192  13.0885  8.8381$r^{2}$  80 
Outer Core  0.546  12.5815  1.2638$r$  3.6426$r^{2}$  5.5281$r^{3}$  10 
Lower Mantle  0.895  7.9565  6.4761$r$ + 5.5283$r^{2}$ 3.0807$r^{3}$  5 
Transition Zone  0.937  11.2494  8.0298$r$  5 
Upper Mantle  0.965  7.1089  3.8045$r$  5 
LVZ  0.996  2.691 + 0.6924$r$  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=0$, to that temperature and the surface to the effective surface temperature of Earth ($255$ K) as calculated via the StefanBoltzmann equation ((2)). Since the planets are assumed to be Earthlike, 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 ($C_{p}\sim 1$). A summary is given in Table 1.
By treating planets purely as conductively heated bodies, we avoid the laborious computation of selfgravitational 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.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 timeevolution of the temperature field within the simulated planet, we first considered the one dimensional heat equation:
$\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
$\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 $N$ concentrically stacked shells. The innermost and outermost shells, $S_{0}$ and $S_{N}$, are at fixed boundary temperatures. For each shell in the interior region, $S_{i}$, its new temperature is computed via Equation 5
$T(S_{i})=\frac{w_{i1}T(S_{i1})+w_{i+1}T(S_{i+1})}{w_{i1}+w_{i+1}}$  (5) 
Where the weights, $w_{i}$, are given by
$w_{i}=4\pi\alpha_{i}r_{i}^{2}$  (6) 
Starting at the $1$st shell ($i=1$), the new temperature was recursively computed for each $i$th 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 reintroduce 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. solutiontime).
It is assumed in our simulation that all planets begin in a volcanic state, and become nonvolcanic when a thermodynamic steadystate is reached. The equilibrium time, $t_{eq}$, at which a planet becomes nonvolcanic, 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\%$ between the $(i1)$th and $i$th iteration, then $t_{eq}$ was set to $t_{i}$.
3 Results
The timeevolved temperature field for a planet simulated at $1.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.5$ billion years, must correspond to a number of iterations less than $t_{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 uniformtemperature 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 $t_{eq}$ is considerably less than that of the Earth. Figure 4 demonstrates that the Moon cools rapidly in comparison to Earth, thus reaching a steadystate 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.5$ billion years, must be greater than $t_{eq}$ for the Moon.
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 $t_{eq}$ falls between the values obtained for Earth and the Moon. If it is given that Mars is currently near its $t_{eq}$ in realtime, then the simulated $t_{eq}$ may be set to $4.5$ billion years ((6)). Knowing that our simulations start at the formation of the planet, and that a linear timescale is being used, we now infer a mapping between realtime and simulationtime:
$t_{ratio}=\frac{6771\text{ steps}}{4.5\text{ billion years}}$  (7) 
Using this result, it is now possible to predict the timeevolved volcanic activity of a planet of arbitrary mass by simulating the evolution of it’s temperature field, and then mapping simulationtime to realtime. The results of this mapping are shown in Figure 6.
Figure 6 shows that our simulation correctly predicts that the mass of a planet is positively correlated with the time until thermal steadystate is reached, which corresponds to the end of volcanic activity on the planet. Our simulation also correctly predicts that, at the present realtime ($4.5$ billion years since planet formation), the Moon and Mercury are volcanically inactive while Venus and Earth are volcanically active. Using the mapping between realtime and simulationtime, the simulation predicts that Earth will continue its volcanic activity until $6$ billion years of age.
When simulating planets beyond the threshold of 2 Earth masses, it was found that the temperature field does not reach steadystate. 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 Earthlike planets were performed only up to 2.6 Earth masses, suggesting that simulation of terrestrial planets above this mass is nontrivial ((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 Earthlike planets was studied by equating volcanic inactivity to thermal steadystate. The timeevolution of planetary temperaturefields was solved using a novel method, in which a Laplace equation was solved iteratively and a mapping between realtime and simulationtime was found by comparison to realworld data.
It was found that Earth is still in its nonsteadystate 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 nonterrestrial behaviour due to continual selfinduced internal temperature increase, in the absence of convection.
Future work in this area could include the nonlinear scaling of density and conductivity profiles and the inclusion of a simple model of convective effects to determine how large Earthlike planets could be simulated to exhibit terrestrial behaviour.
References
 The preliminary reference earth model. Phys. Earth Plan. Int 25, pp. 297–356. Cited by: 1, 2.
 The physics of atmospheres. Cambridge University Press. Cited by: 2.
 Geodynamics and rate of volcanism on massive earthlike planets. Astrophys.J. 700, pp. 1732–1749. Cited by: 3.
 Highresolution simulations of the final assembly of earthlike planets 2: water delivery and planetary habitability. Astrobiology Special Issue on M Stars 7. Cited by: 2, 3.
 Electrical and thermal conductivities of feni 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: \urlhttp://www.sciencedaily.comÂ/releases/2001/11/011109075016.htm Cited by: 3.
Additional Assets
 moon_heatmap.png 101 KB
 numerical_solution.png 12.5 KB
 mass_vs_time.png 27 KB
 mass35_3d.png 415 KB
 earth_3d.png 399 KB
 mars_heatmap.png 174 KB
 earth_heatmap.png 246 KB
Showing 1 Reviews

1
Note: This paper was written for a 48hr undergraduate competition (http://www.uphysicsc.com/), in which the research, coding and writing were done within that time frame.
(Stars below are arbitrary.)
License
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.
Wow! This must be a record of some sort!
Wow! This must be a record of some sort!
You must try it after get interrupt your main software audiodg. I know it always remain some record for what options you choose.