# Computational Modeling of Volcanism on Earth-like Planets

### Abstract

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 $\sim 6.0$ billion years since its formation, while Earth-like 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 Earth-like 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 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,

 $\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 steady-state 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 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 ($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 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.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:

 $\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_{i-1}T(S_{i-1})+w_{i+1}T(S_{i+1})}{w_{i-1}+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 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, $t_{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\%$ between the $(i-1)$th and $i$th iteration, then $t_{eq}$ was set to $t_{i}$.

### 3 Results

The time-evolved 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 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 $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 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.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 real-time, 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 time-scale is being used, we now infer a mapping between real-time and simulation-time:

 $t_{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 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.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 $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 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.

## References

• 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: \urlhttp://www.sciencedaily.comÂ­/releases/2001/11/011109075016.htm Cited by: 3.

• 101 KB
• 12.5 KB
• 27 KB
• 415 KB
• 399 KB
• 174 KB
• 246 KB

### Showing 1 Reviews

• Edward Kim
Originality of work
Quality of writing
Quality of figures
Confidence in paper
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.)

This review has 2 comments. Click to view.
• Joshua Nicholson

Wow! This must be a record of some sort!

• Joshua Nicholson

Wow! This must be a record of some sort!

• ankit saini

You must try it after get interrupt your main software audiodg. I know it always remain some record for what options you choose.