#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Wed Apr 11 2023 @author: dbyrne @email:djbyrne20@gmail.com This file contains the python functions used in the manuscript "Primordial and 244Pu-derived xenon missing from Earth's convecting mantle" By D.V. Bekaert et al., submitted to Nature in April 2023 """ import numpy as np import matplotlib.pyplot as plt def Nres(eta): return k(eta)/(Mmantle*0.9) def k(eta): return (Qp_input/eta)*(np.exp(eta*TE*1e6)-1) def decay(i0, t, lam): return i0*np.exp(-lam*t) def Q(t, eta, Qp, f=1): #processed mantle flux in g/yr Qexp = Qp*np.exp(eta*(TE-t)*1e6) TE_yr = TE*1e6 m = (2*Qp/(TE_yr**2)) * ( ((np.exp(eta*TE_yr)-1)/eta) - TE_yr) Qlin = m*(TE-t)*1e6 + Qp return f*Qexp + (1-f)*Qlin def gen_degas_profile(t, convecting_mantle_fraction, eta, Qp, f=1): #generate a fractional degassing profile return Q(t, eta, Qp, f)*1e6/(Mmantle*convecting_mantle_fraction) def regas_Xe130(degas_profile, mod_Xe130, f_subXe): #calculate initial Xe130 concentration in the mantle in mol/g initial_Xe130 = mod_Xe130*(1-f_subXe)/np.prod((1-degas_profile)) return initial_Xe130 def prop_CC(Xe130): return Xe130*Mmantle/(ME*CC_Xe130) def sigmoid(z, sharpness=1, inflection=0, skew=0): sigmoid_curve = 1/(1+np.exp(-sharpness*(z-inflection)/(z**skew))) return sigmoid_curve def gen_crustalgrowth_profile(t, sharpness=1, inflection=0, skew=0): sigmoid_profile = sigmoid(t, sharpness, inflection, skew) sigmoid_profile[0] = 0 growth_profile = sigmoid_profile/sigmoid_profile[-1] return growth_profile def U_evolution(t, crustalgrowth_profile): #Returns U evolution within the mantle, crust in mol/g initial_BSE_U = decay(U238_BSE, -TE, lam238)*MBSE fextracted = crustalgrowth_profile BSE_U = decay(initial_BSE_U, t, lam238) fextracted = crustalgrowth_profile*U_fcrust mantle_U = BSE_U*(1-fextracted)/Mmantle crustal_U = BSE_U*(fextracted)/Mcrust return mantle_U def Pu_evolution(t, crustalgrowth_profile): #Returns Pu evolution within the mantle, crust in mol/g initial_BSE_U = decay(U238_BSE, -TE, lam238)*MBSE initial_BSE_Pu = initial_BSE_U * CC_Pu_U fextracted = crustalgrowth_profile BSE_Pu = decay(initial_BSE_Pu, t, lam244) fextracted = crustalgrowth_profile*U_fcrust mantle_Pu = BSE_Pu*(1-fextracted)/Mmantle crustal_Pu = BSE_Pu*(fextracted)/Mcrust return mantle_Pu def I_evolution(t, crustalgrowth_profile, I_ppb): initial_BSE_I129 = MBSE * I_ppb * I_ratio * 1e-9 / 127 fextracted = crustalgrowth_profile*0.95 BSE_I129 = decay(initial_BSE_I129, t, lam129) mantle_I129 = BSE_I129*(1-fextracted)/(Mmantle) crustal_I = BSE_I129*(fextracted)/Mcrust return mantle_I129 def Xe136_ingrowth(t, crustalgrowth_profile, element="both"): U = U_evolution(t, crustalgrowth_profile) Pu = Pu_evolution(t, crustalgrowth_profile) ingrowth_from_U = (U-decay(U, 1, lam238))*Ysf_U_Xe136 ingrowth_from_Pu = (Pu-decay(Pu, 1, lam244))*Ysf_Pu_Xe136 output = np.zeros(np.size(t)) if element == "U": output += ingrowth_from_U elif element == "Pu": output += ingrowth_from_Pu elif element == "both": output += ingrowth_from_Pu output += ingrowth_from_U return output def Xe129_ingrowth(t, crustalgrowth_profile, I_ppb): I = I_evolution(t, crustalgrowth_profile, I_ppb) ingrowth_from_I = (I-decay(I, 1, lam129)) return ingrowth_from_I def Xe_evolution(t, initial_Xe130, degas_profile, crustalgrowth_profile, I_abundance, misfit_evaluation=False): initial_Xe136 = initial_Xe130 * CC_Xe136_Xe130 #chondritic 136/130 initial_Xe129 = initial_Xe130 * CC_Xe129_Xe130 #chondritic 129/130 fiss_Xe136_Pu = Xe136_ingrowth(t, crustalgrowth_profile, element='Pu') fiss_Xe136_U = Xe136_ingrowth(t, crustalgrowth_profile, element='U') rad_Xe129_I = Xe129_ingrowth(t, crustalgrowth_profile, I_abundance) Xe130 = np.zeros(np.size(t)) Xe136 = np.zeros(np.size(t)) Xe136_Pu = np.zeros(np.size(t)) Xe136_U = np.zeros(np.size(t)) Xe129 = np.zeros(np.size(t)) Xe129_I = np.zeros(np.size(t)) Xe130[0] = initial_Xe130 Xe136[0] = initial_Xe136 Xe129[0] = initial_Xe129 for i in range(1, len(t)): f_degas = 1 - degas_profile[i] Xe130[i] = Xe130[i-1] * f_degas Xe136[i] = (Xe136[i-1] + fiss_Xe136_Pu[i] + fiss_Xe136_U[i]) * f_degas Xe136_Pu[i] = (Xe136_Pu[i-1] + fiss_Xe136_Pu[i]) * f_degas Xe136_U[i] = (Xe136_U[i-1] + fiss_Xe136_U[i]) * f_degas Xe129[i] = (Xe129[i-1] + rad_Xe129_I[i]) * f_degas Xe129_I[i] = (Xe129_I[i-1] + rad_Xe129_I[i]) * f_degas if misfit_evaluation: return Xe130[-1], Xe136[-1], Xe129[-1], Xe136_Pu[-1], Xe136_U[-1], Xe129_I[-1] else: return Xe130, Xe136, Xe129, Xe136_Pu, Xe136_U, Xe129_I def add_giantimpact(degas_profile, start, end, magnitude): modified_degas_profile = degas_profile modified_degas_profile[start:end] = magnitude return modified_degas_profile def model_misfit(model_output): Xe130, Xe136, Xe129, Xe136_Pu, Xe136_U, Xe129_I = model_output Pu_TF = Xe136_Pu/(Xe136_Pu+Xe136_U) I_Pu = Xe129_I/Xe136_Pu misfit = 0 misfit += ((Xe136/Xe130) - X_Xe136_Xe130)**2 / X_Xe136_Xe130 misfit += ((Xe129/Xe130) - X_Xe129_Xe130)**2 / X_Xe129_Xe130 misfit += (Pu_TF - X_Pu_TF)**2 / X_Pu_TF misfit += (I_Pu - X_I_Pu)**2 / X_I_Pu return misfit def print_model(model_output): Xe130, Xe136, Xe129, Xe136_Pu, Xe136_U, Xe129_I = model_output Pu_TF = Xe136_Pu/(Xe136_Pu+Xe136_U) I_Pu = Xe129_I/Xe136_Pu fig = plt.figure(figsize=(8, 10)) ax1 = fig.add_subplot(411) ax2 = fig.add_subplot(412) ax3 = fig.add_subplot(413) ax4 = fig.add_subplot(414) ax1.plot(time, Xe136/Xe130, 'g-') ax2.plot(time, Xe129/Xe130, 'm-') ax3.plot(time, Pu_TF, 'b-') ax4.plot(time, I_Pu, 'r-') ax1.plot([0, TE], [X_Xe136_Xe130, X_Xe136_Xe130], 'g--') ax2.plot([0, TE], [X_Xe129_Xe130, X_Xe129_Xe130], 'm--') ax3.plot([0, TE], [X_Pu_TF, X_Pu_TF], 'b--') ax4.plot([0, TE], [X_I_Pu, X_I_Pu], 'r--') ax1.set_xticks([]) ax2.set_xticks([]) ax3.set_xticks([]) ax1.set_ylabel(r"$^{136}Xe/^{130}Xe$") ax2.set_ylabel(r"$^{129}Xe/^{130}Xe$") ax3.set_ylabel(r"$Pu/TF$") ax4.set_ylabel(r"$I/Pu$") ax4.set_xlabel(r"time (Myr)") plt.show() def run_model(eta, I, modXe, fexp, impact=False, t_impact=0, mag_impact=0.05): degas = gen_degas_profile(time, CMF_input, eta, Qp_input, f=fexp) if impact: degas = add_giantimpact(degas, t_impact, t_impact+10, mag_impact) Xe_init = regas_Xe130(degas, modXe, f_subXe_input) return Xe_evolution(time, Xe_init, degas, cr, I, misfit_evaluation=True) #Setting up base constants time = np.linspace(0, 4567, 4568) TE = 4567 #Age of the Earth / Myr ME = 5.98e27 #Mass of Earth in g Mmantle = 4e27 #mass of mantle in g Mcrust = 2.2e25 #mass of continental crust in g MBSE = Mmantle+Mcrust Qp_input = 6.1e17 #modern mantle processing rate in g/yr eta_input = 8.1e-10 #mantle processing rate change with time CMF_input = 0.9 #fraction of convecting mantle mod_Xe130_input = 1.79e-18 #modern mantle Xe130 abundance in mol/g mod_Xe136_input = mod_Xe130_input * 2.44 mod_Xe129_input = mod_Xe130_input * 7.41 f_subXe_input = 0.8 #fraction of subducted Xe130 in modern mantle mod_Xe130_corrected = mod_Xe130_input*(1-f_subXe_input) mod_Xe136_corrected = mod_Xe136_input -\ (mod_Xe130_input*f_subXe_input*2.1761) mod_Xe129_corrected = mod_Xe129_input -\ (mod_Xe130_input*f_subXe_input*6.4963) X_Xe129_Xe130 = mod_Xe129_corrected/mod_Xe130_corrected X_Xe136_Xe130 = mod_Xe136_corrected/mod_Xe130_corrected #CC abundances CC_Xe130 = 5.38e-14 # mol/g CC_Xe129_Xe130 = 6.533 CC_Xe136_Xe130 = 1.988 CC_Pu_U = 0.0068 I_ratio = 1.1e-4 # initial ratio of 129/127 I #modern composition U_BSE = 0.021/(1e6*238.029) #molg-1 U238_BSE = U_BSE*0.99274 U235_BSE = U_BSE*0.0072 U_crust = 1.3/(1e6*238.029) U238_crust = U_crust*0.99274 U235_crust = U_crust*0.0072 U_fcrust = U_crust*Mcrust/(U_BSE*MBSE) #decay constants lam238 = np.log(2)/4468 lam235 = np.log(2)/704 lam129 = np.log(2)/15.7 lam244 = np.log(2)/82 Ysf_U_Xe136 = 6.83E-12/lam238 Ysf_Pu_Xe136 = 5.96E-07/lam244 #Change following constants to desired target Pu/TF and I/Pu #Etna values are 0.03 and 80; MORB are 0.26 and 11 X_Pu_TF = 0.03 X_I_Pu = 80 #%% Example model run #Set model parameters eta = 8e-10 #parameter controlling total degassing I = 10 #BSE I abundance in units of ppb modXe = 3e-18 #modern Xe130 mantle abundance in mol/g fexp=0.5 #shape of degassing curve between exponential (1) and linear(0) #generate crustal growth profile cr = gen_crustalgrowth_profile(time, sharpness=0.03, inflection=1250, skew=0.35) #generate degassing profile dg = gen_degas_profile(time, CMF_input, eta, Qp_input, f=fexp) #calculate initial Xe130 abundance Xe_init = regas_Xe130(dg, modXe, f_subXe_input) #run Xe evolution forward model model = Xe_evolution(time, Xe_init, dg, cr, I) #print mantle Xe isotope evolution with time print_model(model)