# How to calculate the loss of primordial heat for a given age and mass of a planet?

I´m working on a computer program which generates realistic planets for worldbuilding purposes. I´ve run into an issue with calculating my planets heat loss. While radiogenic and tidal heat generation and loss were easy enough to pin down, calculating primordial heat seems way harder. On Earth, we currently lose 12 - 30 TW of primordial heat. Jupiter on the other hand is even today 40 K hotter than it should be due to its primordial heat. What I'm looking for is a reasonably accurate formula giving me the current heat loss in TW for a given age and mass of the planet.

Assume that a planet is simply a pure rocky sphere without an atmosphere. If it isn't too complicated the modification of the formula for a water or pure hydrogen sphere would be nice.

• For a hard-science estimate, that's a tough one. Heat loss is depending on many factors, including nuances of planetary composition. For example, I bet the figures would be substantially different for Earth, "dry Earth" and "water Earth". – Alexander Jun 6 at 19:32
• A certain Lord Kelvin attempted the same, without taking into account the radioactive decay en.wikipedia.org/wiki/… – L.Dutch - Reinstate Monica Jun 6 at 20:06
• @Alexander I know that, but since I'm not writing a physics simulation, just a program working on a "close enough" basis, I gave the baseline of the planet in the last paragraph. Atmosphere and oceans shall simply be ignored. – TheDyingOfLight Jun 6 at 20:21
• @L.Dutch I'm aware of Lord Kelvins story and I already got the other heat sources (tidal, solar and radiogenic) pinned down with reasonable accuracy. – TheDyingOfLight Jun 6 at 20:23
• Other heat sources is not the only mistake Lord Kelvin did. Convection heat transfer works very differently inside the planet depending on temperature and composition. Just radiation of heat through ocean vs continental crust is 50% different. – Alexander Jun 6 at 20:38

Like most things in science, the answer to "how do I model x" is "it depends how accurate you want to be." The simplest thing would be to pick some number for what you think the surface temperature of the planet should be, and then use the Stefan-Boltzmann law to find out how much the surface radiates. Of course, the problem is that we likely don't know the planet's surface temperature. If you want your model to only be for soon after the planet forms, we can do pretty well by assuming the planet is a sphere of uniform temperature caused by all of the gravitational binding energy.

If it supposed to apply for a time long after the formation of the planet, looking at surface temperatures for other planets in our solar system for inspiration might be a good starting point. However, you have to keep in mind that you'll probably want to look at temps for the dark side of the planet, because trying to solve this for planets where one side is in the sun makes it significantly harder as the symmetry of the problem is destroyed. So, for the rest of the answer I'm going to assume the planet in question is far enough away from its star that there isn't such more light coming from the direction of it's star than there is from any other direction in space.

What follows is a more accurate, but much more complicated model. Be warned: this is a fairly sophisticated mathematical problem and requires mathematical experience up to at least partial differential equations. I'm not sure what your mathematical experience is, so if you expand on that in the comments I can add a section more tailored to you.

With that out of the way, I think L.Dutch has an answer that's on the right track-- this is just attempting to flesh it out a bit more. First off, to model how the heat moves through Earth, we use the heat equation:

$$\frac{\partial u}{\partial t} = \alpha \nabla^2u + \frac{q}{c_p \rho}$$

where $$u$$ is the heat density, $$c_p$$ the specific heat capacity, $$\rho$$ the density, and q is the heat generated\lost in the material due to outside sources. Note that the temperature is related to these quantities by $$T=\frac{u}{c_p \rho}$$.

Now, technically this is all you need, along with some boundary values. In practice, it might be unclear what some of these values should correspond to, so I'll walk through that. First, it will help simplify stuff if we assume spherical symmetry. This isn't a requirement, and if we want to model a planet close to it's star or be super duper accurate we can't do it. But, since we assumed it wasn't, and since this is already complicated enough, this allows us to make the substitution

$$\nabla^2 \rightarrow \frac{1}{r^2}\frac{\partial}{\partial r}(r^2 \frac{\partial}{\partial r})$$

and the assumption that $$u = u(r, t)$$ (ie there is no $$\theta$$ or $$\phi$$ dependence). Now, to avoid boundary condition problems, we will make the spatial domain of our problem all of $$\mathbb{R}^3$$. But since this means the problem includes space, which can't conduct heat, we must make sure that $$\alpha(r)=0$$ when $$r$$ is greater than the radius of our planet.

Next, we have to determine the form for $$q$$. This will either make complete sense or none at all and require a ton of explanation, so I'll just leave the result which is that

$$q=-4\pi R_p^2 \sigma \epsilon (\frac{u}{c_p \rho})^4\delta(r-R_p) + \mathbf{1}_{r \in [0, R_p]}D_0r^2e^{-t/\lambda}$$

where the left term accounts for radiative cooling of the planet and the right for radioactive generation of heat within the planet as per L.Dutch's answer. All the variable names are the standard ones from the Stefan-Boltzmann law, $$R_p$$ is the radius of the planet, and $$\mathbf{1}$$ is the indicator function.

Finally all that's left is to specify boundary conditions, which should just specify $$u(r,0)$$ and that $$u$$ is bounded. The exact form of these will be that the temperature of space (ie $$r>R_p$$) will be the temperature of the cosmic microwave background, and that inside it will be some distribution of the initial binding energy of the planet (uniform is easiest and probably pretty accurate). I think these conditions should be enough to solve the equation, but it's possible that the non-linear $$u^4$$ term messes stuff up and the equation is no longer parabolic, in which case I don't have enough math knowledge to help.

After all that, you "simply" have to solve this monstrosity of an equation. A numerical solver is probably easiest, but you could attempt an analytical solution. Unfortunately, the non-linear $$u^4$$ term takes away the strongest tool in our arsenal (Green's functions) so I don't have any clue how you'd proceed, or if an analytical solution is even possible (it likely isn't).

As one final note, I should say that I've made a few simplifying assumptions. Aside from the assumption of spherical symmetry, I've also implicitly assumed that there is only one radioactive species and that the density of rock is constant throughout the planet, among others. However, if you understand any of this and desire a more accurate approach, I should have hopefully given you enough to work with to extend it to a more accurate model.

Like I said before, this is a pretty complicated problem so don't be worried if you don't understand this approach-- I just wanted to leave a record for how this could be solved in an accurate manner. My apologies if it is also rather terse-- I don't have time right now to add more explanations. If I have more time later, I might edit some in.

• +1 for including your assumptions – awsirkis Jun 9 at 16:39
• Thanks for the great answer and sorry for the late reaction. I´m a German high school student who just finished school and while I´ve always been decent at math your formula is somewhat beyond me. Would you mind explaining the variables of the final equation in detail and elaborate on what you mean by a numerical solver? – TheDyingOfLight Jun 16 at 17:31
• How do I fit Tidal heating in there? This post is my basis for calculating it and I assume that I can simply add it to the formula giving q after I´ve adjusted the units. I that correct? – TheDyingOfLight Jun 16 at 17:35
• @TheDyingOfLight No problem! I'll add a more thorough description sometime if I have more time, but I used all the standard variable names so you should be able to get more info about what the parameters mean by looking at the wikipedia links I provided and L. Dutch's answer. As for a numerical solver, basically what that means is that you use a computer to turn the continuous problem involving derivatives into a discrete one using differences, which the computer can then solve by performing a bunch of arithmetic. Unfortunately, this can get a little difficult for PDE's so finding the right... – el duderino Jun 18 at 1:27

## Heat is fungible.

If a planet was once very hot, and has been radiating heat, and has been generating heat by nuclear fission, then the following are indistinguishable, and both somewhat misleading:

• The planet is staying warm because as the primordial heat radiates away it's replaced by radiogenic heat.
• The planet's primordial heat is being preserved while it radiates away radiogenic heat.

## Primordial heat is (almost) non-renewable.

Basically the idea here is that primordial heat comes from gravitational potential energy.

To a first approximation, we could just use the gravitational binding energy. The binding energy is usually defined as the energy required to rip an object apart, but in the case of gravity we can equivalently say that it's the energy released when the object clumps together.

You may also need to account for the (small amount of) heat the materials had prior to clumping together, the non-uniform density of the planet, and how the average density and the "density curve" (out from the center) change over time. Celestial impacts and sustained meteor bombardments may also be thought of as adding "primordial" heat.

## How are you actually modeling this stuff?

Specifically, the way you've worded the question suggests that you're accounting for heat loss on a per-source basis. If you have a good-enough model for heat loss (by approximately-black-body radiation for example), then use that for all heat loss as a single out-flow of energy.

You'll have a few in-flows (which will not be evenly distributed through the body of the planet): nuclear fission, solar radiation, tidal friction, and "primordial". That said; for a rocky middle-age planet, no "primordial" energy will be generated at all.
So long as the planet isn't changing mass or size, and the distribution of mass within the planet isn't changing either, no gravitational potential energy is being converted into heat.

• Heat conductivity is a major issue. A convective mantle makes for a warmer surface, greatly increasing the heat loss. Also the core density increases upon cooling, which sets free additional potential energy as heat. – Karl Jun 7 at 0:44
• @Karl Indeed probably pretty hard to predict this too not to mention that if you want it to be accurate over geological timescales then you probably really need to start modelling volcanism too sure a minor eruption might not make any real difference but those flood basalt eruptions that bring literally millions of cubic kilometers of hot magma directly to the surface are likely to be less trivial. And we are not exactly great at predicting such geological events as yet so the available formulas might leave something to be desired here. – MttJocy Jun 7 at 12:05
• The stochasticity of volcanism doesn't seem like a big problem; no eruption will be a large fraction of the planet's mass, so on the timescales in question you can just use the average rate of volcanism. Probably it would be right to treat it as just another form of convection. – ShapeOfMatter Jun 7 at 13:35
• That said, I haven't tried here to go into how to model heat loss. On the one hand that's hard and I'm not qualified, on the other hand it sounds like the OP already has some system they're comfortable with. (L.Dutch's answer is correct that if you can somehow know the surface temperature, and make reasonable approximations about it's albedo and emissivity, then calculating the loss is easy.) – ShapeOfMatter Jun 7 at 13:41

The energy radiated per unit surface by a black body at surface temperature $$T(t)$$ at time t can be calculated by Stefan-Boltzmann law

$$P(t)=\epsilon\sigma T^4(t)$$

The temperature T depends by the amount of energy $$H(t)$$ stored in the planet, this being given by two components: the gravitational binding energy and the radioactive decay.

For the gravitational binding energy:

A gravitational binding energy is the minimum energy that must be added to a system for the system to cease being in a gravitationally bound state. A gravitationally bound system has a lower (i.e., more negative) gravitational potential energy than the sum of its parts—this is what keeps the system aggregated in accordance with the minimum total potential energy principle.

For a spherical mass of uniform density, the gravitational binding energy U is given by the formula $$U=3GM^2 \over 5R$$ where G is the gravitational constant, M is the mass of the sphere, and R is its radius

For the radioactive decay, assuming a decay constant $$\lambda$$ and an initial emission rate $$D_0$$:

$$D(t)=D_0e^{-t/\lambda}$$ (if you have more than a single radioactive species, that becomes a summation over the various decay rates and initial concentrations)

This so far allows you to estimate the amount of heat available in the mass of the planet. The main problem is, how to correlate the surface temperature to that amount of thermal energy still present below the surface.

The empiric way to estimate it is to measure the thermal flow at different depths, and try to fit those points with some heat transfer relation (remember the Kola deep drill?).

• I think you've miss-transcribed the Stefan-Boltzmann law? There shouldn't be a time component, and it certainly shouldn't be increasing with the forth power of time. – ShapeOfMatter Jun 7 at 11:37
• @ShapeOfMatter, forth power of the absolute temperature, and since it is a black body there is no emissivity – L.Dutch - Reinstate Monica Jun 7 at 11:44
• @L.Dutch. NahUh. Although planets and stars are neither in thermal equilibrium with their surroundings nor perfect black bodies, black-body radiation is used as a first approximation for the energy they emit.[6] They are called 'black bodies" because at temps below ~1000°K they do not emit significant visible light. And yes time: The radiant emittance {\displaystyle j^{\star }}￼ has dimensions of energy flux (energy per time per area), and the SI units of measure are joules per second per square metre, or equivalently, watts per square metre. – pHred Jun 7 at 12:31
• I'm not sure I follow what @pHred is saying, but they're right that the Power $P$ has a time factor in its units. That said, if you know the temperate $T$, you don't care how much time has elapsed. I think the formula on Wikipedia is correct (and if it's not let's take the fight over there). $P = A \sigma T^4$ (there's no $t$ at all) – ShapeOfMatter Jun 7 at 13:31
• @ShapeOfMatter L.Dutch is right, he's just explicitly making clear in his answer that the temperature is a function of time, which it will be because the planet will slowly cool off. To actually find the rate that this happens you'd use energy conservation to write a differential equation and then solve it. – el duderino Jun 7 at 14:53