Half-life of Silicate Earth heating

Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 14 January 2019

If we wanted to represent the decay of radiogenic heat production from several major radionuclides (238U, 232Th, 40K) by one simple exponential, what would be the representative half-life?

Constraints:

  1. best fits the actual sum of individual exponential decays
  2. gives the actual present-day heat production rate (in 1-parameter fit)
In [1]:
import numpy as np
from math import log
import matplotlib.pyplot as plt
In [2]:
plt.rcParams['figure.figsize'] = (12,8)
In [3]:
Mearth = 5.97218e24 # mass of Earth in kg (Chambat et al 2010)
Mcore = 1.93e24 # mass of core in kg (Yoder 1995)
Mbse = Mearth-Mcore # mass of BSE in kg

nuclides = ["238U", "235U", "232Th", "40K"]
halflives = np.array([4.4683, 0.70348, 14.1, 1.262]) # in Gyr
abundances = np.array([20.0e-9, 20.0e-9, 75.4e-9, 280e-6]) # in kg-element per kg-rock
powerspec = np.array([94.247, 4.046, 26.180, 0.003387]) * 1e-6 *1e-12 # in TW per kg-element
H0 = abundances*powerspec*Mbse # H0[] in TW in BSE
H0sum = sum(H0) # H0 in TW in BSE
lambdas = log(2.)/halflives
In [4]:
print ("Th/U mass ratio = {:.3f}".format(abundances[2]/abundances[0]))

print ("K/U mass ratio = {:.0f}".format(abundances[3]/abundances[0]))

heat_production_BSE = H0sum # in TW
print ("CHECK: Heat production in BSE = {:.2f} TW".format(heat_production_BSE))

Th/U mass ratio = 3.770
K/U mass ratio = 14000
CHECK: Heat production in BSE = 19.76 TW
In [5]:
age = np.linspace(0, 4.568, 128) # even time grid over the age of the Earth
In [6]:
# heat production vs. age: sum of exponentials
Hage = np.zeros(age.size)
for i in range(halflives.size):
    Hage = Hage + H0[i]*np.exp(lambdas[i]*age)
In [7]:
from scipy.optimize import curve_fit

1-parameter fit: only decay constant

In [8]:
# 1 parameter fit – only decay constant

def func(a, lam):
    return H0sum*np.exp(lam*a)

p_opt, p_cov = curve_fit(func, age, Hage)

print (p_opt)
print (p_cov)
lam_opt = p_opt[0]
lam_err = np.sqrt(p_cov[0][0])

print ("\ndecay constant fit = ({:.4f} ± {:.4f}) Gyr-1\n".format(lam_opt, lam_err))

halflife_fit = log(2.)/lam_opt
print ("half-life fit = {:.4f} Gyr".format(halflife_fit))

[ 0.32932529]
[[  5.50184199e-06]]

decay constant fit = (0.3293 ± 0.0023) Gyr-1

half-life fit = 2.1047 Gyr

2-parameter fit: decay constant and H0

In [9]:
# 2 parameter fit – decay constant and H0

def func2(a, lam, h0):
    return h0*np.exp(lam*a)

p_opt, p_cov = curve_fit(func2, age, Hage)
print (p_opt)
print (p_cov)

lam2_opt = p_opt[0]
lam2_err = np.sqrt(p_cov[0][0])
H0_opt = p_opt[1]
H0_err = np.sqrt(p_cov[1][1])

print ("\ndecay constant fit = ({:.4f} ± {:.4f}) Gyr-1".format(lam2_opt, lam2_err))
print ("present-day H0 = ({:.4f} ± {:.4f}) TW\n".format(H0_opt, H0_err))

halflife2_fit = log(2.)/lam2_opt
print ("half-life fit = {:.4f} Gyr".format(halflife2_fit))

[  0.40150962  15.16906412]
[[  2.35913110e-05  -1.23735667e-03]
 [ -1.23735667e-03   7.03165979e-02]]

decay constant fit = (0.4015 ± 0.0049) Gyr-1
present-day H0 = (15.1691 ± 0.2652) TW

half-life fit = 1.7264 Gyr
In [10]:
Hfit = H0sum*np.exp(lam_opt*age)
Hfit2 = H0_opt*np.exp(lam2_opt*age)

plt.plot(age, Hage, label='actual')
plt.plot(age, Hfit, label='1-parameter fit')
plt.plot(age, Hfit2, "--", label='2-parameter fit')
plt.xlabel("Age in Gyr")
plt.ylabel("Heat production in TW")
plt.legend()
plt.show()