#gradual accretion and differentiation code used in Dodds et al., 2020, JGR Planets. #requires Python 3 to run #Input parameters are R0, Rf, t0, dt_acc. import numpy as np import time as TIME import random import matplotlib.pyplot as plt w_max=1 w=-1 datafile_path="LONGDUR.dat" while w=w_max: break print(w+1,'/',w_max) R0_int=random.randrange(60,600,5) R=(random.randrange(R0_int+5,605,5))*1E3 #final, compacted radius phi_0=0.25 Rf=((1-phi_0)**(-1/3))*R #final uncompacted radius exp_data[2,0]=R0 exp_data[3,0]=Rf exp_data[29,0]=R bigg=6.67E-11# rhom=3000 rho_0=(1-phi_0)*rhom cpm=800 cpc=850 diffm=9E-7 rho_fe=7800 eta_fe=1E-2 xcount_max=10 check2=0 #for reynolds number t_spin=10*60*60 mag_diff=1.3 j_undiff=-10 #radioactive decay of 26al h0=0.355 c0=7E-7 thalf=0.717E6*365*24*3600 #iron melting properties Ts_fe=1234# K, eutectic melt point Tl_fe0=1810# K, pure iron melt point cb=31#initial wt% S exp_data[4,0]=cb ceu=32 #eutectic wt% S Lfe=270E3 #JKg-1, iron latent heat Tl_fe=Tl_fe0#-(18*cb) Tl_m=Tl_fe0-(18*cb) Tl_fe=Tl_m exp_data[5,0]=Tl_m m=-18 chi_crit=cb/ceu rho_l_fe_0=6980 m_mol_S=32 m_mol_Fe=56 v0_fe=0.4 #from CM chondrites, measured volume fraction of iron x0_fe=(rhom*v0_fe/rho_fe)*((1-(v0_fe*(1-(rhom/rho_fe))))**(-1)) #compaction constants A=3.8E-5 Ea_comp=60*4184 #J/mol Tsquish=710 #K, currently a constant, will need to change Rgas=8.3145 b=1E-3 #grain size - av chondrule size #conductivity constants - Krause et al, 2011, LPSC kb=diffm*rhom*cpm#Wm-1K-1 phi_1=0.08 phi_2=0.17 k_fe=30 t0_int=random.randrange(5E4,2E6,1E4) t0=t0_int*365*24*60*60 dt_acc_int=random.randrange(5E2,(4.5E6-t0_int),1E4) dt_acc=dt_acc_int*365*24*60*60 exp_data[0,0]=t0/(1E6*365*24*3600) exp_data[1,0]=dt_acc/(1E6*365*24*3600) print('R0',R0_int,'R',round(R/1E3),'t0',round(exp_data[0,0],3),'dt_acc',round(exp_data[1,0],3)) dr_max=700 dt=2E10 #s, need to check this for possible set of parameters if exp_data[0,0]<0.5: dr_max=700 dt=2E10 if R0_int<150: dr_max=700 dt=2E10 dr_min=((1-phi_0)**(1/3))*dr_max i_max=int(R/dr_min) i_i=int(R0/dr_max) T0=200#K initial temperature of material t_end=150E6*365*24*60*60 t_end2=12E6*365*24*60*60 j_max=int((t_end-t0)/dt) j_end=int((t_end2-t0)/dt) j_acc=int((t0+dt_acc)/dt) #end of accretion time=np.empty(j_max) spikes=np.empty((j_max,4)) spikes.fill(0) radius=np.empty((i_max,j_max)) rho=np.empty((i_max,j_max)) rho_old=np.empty((i_max,j_max)) r=np.empty(i_max) rho_old2=np.empty(i_max) g=np.empty((i_max,j_max)) P=np.empty((i_max,j_max)) P[:,0]=0 #set initial pressure to zero dP=np.empty(i_max) phi=np.empty((i_max,j_max))# phi.fill(phi_0) T=np.empty((i_max,j_max)) dr=np.empty((i_max,j_max)) dr.fill(dr_max) k=np.empty((i_max,j_max)) k.fill(0) Rp_c=np.empty(j_max) Rp_c[0]=R0 chi=np.empty((i_max,j_max)) chi.fill(0) chi_fe=np.empty((i_max,j_max)) chi_fe.fill(0) critra=np.empty((i_max,j_max)) rayleigh=np.empty((i_max,j_max)) frac_S=np.empty((i_max,j_max)) frac_S.fill(0) frac_S_core=np.empty(j_max) frac_S_core.fill(0) Tpeak=np.empty(j_max) TC=np.empty(j_max) TC.fill(0) TCMB=np.empty(j_max) dT_core=np.empty((i_max,j_max)) dT_core.fill(0) unstable=np.empty((i_max,j_max)) unstable.fill(0) radius_convecting=np.empty(j_max) radius_convecting.fill(0) mag_rey=np.empty(j_max) mag_rey.fill(0) F1=np.empty(j_max) F1.fill(0) F2=np.empty(j_max) F2.fill(0) buoyancy_flux=np.empty(j_max) mag_rey_ave=np.empty(j_max) mag_rey_ave.fill(0) dTcdt=np.empty(j_max) dTcdt.fill(0) dcon_x=np.empty(j_max) dcon_x.fill(0) dcon=np.empty(j_max) time[0]=t0/(1E6*365*24*60*60) #initial radial set up for i in range(0,i_i,1): radius[i,0]=(dr_max*i)+dr_max T[:i_i,0]=T0 phi[:i_i,0]=phi_0 rho[:i_i,0]=rho_0 eta=np.empty((i_max,j_max)) eta[:,0]=4E30 eta.fill(4E30) eta_cmb_old=4E30 chi_sil=np.empty((i_max,j_max)) T_sil_s=1400 T_sil_l=1800 T_sil_crit=1600 a_sil=25 eta_0=1E21 #Pas, reference viscosity eta_1=1E14 #Pas, second reference viscosity eta_sil_l=100 #Pas, viscosity of silicate melts eta_fe_l=1E-2 #Pas a_exp=9.2E-5 #thermal expansivity of core - different value? a_sil_exp=4E-5 #thermal expansivity of silicates E_vis=3E5 #activation energy for silicate viscosity profile Tvis=1400 gamma=E_vis/(Rgas*Tvis*Tvis) convecting=np.empty((i_max,j_max)) TM=np.empty(j_max) surfflux=np.empty(j_max) rho_sil=np.empty((i_max,j_max)) rho_sil[:i_i,0]=rho_0 chi_fe_tot=np.empty((i_max,j_max)) chi_fe_tot.fill(x0_fe) m_l_fe=np.empty(j_max) m_l_fe.fill(0) rcore=np.empty(j_max) rcore.fill(0) f1_ave=np.empty(j_max) f1_ave.fill(0) #new viscosity profile a0=64.17 a1=-1/29 a2=-5 a3=-1625 a4=15 E1_ave=0 time_ave=np.empty(j_max) time_ave.fill(0) x=-1 diffc=k_fe/(rho_l_fe_0*cpc) dt_crit=(dr_min*dr_min)/diffc if dt_crit<=dt: print("CFL not met initially") else: print("CFL met initially") j=0 Rp_u_old=R0 dRp_old=0 i_s0=i_i i_s1=0 dRp_u_sum=0 #discriminant dRp_u=0 i_con=-1 b_int=-1 i_core=-1 i_con_final=-1 j_con=-1 j_fullcon=-1 j_corefull=-1 j_fullsize=-1 j_congrow=-1 jfreeze=-1 j_end=-1 j_cond=-1 j_unstable=-1 j_ave=-1 Rb=0 Rp_u=R0 i_old_2=-1 while j0: i_s1=i_s0+int(dRp_u_sum//dr_max) dRp_u_sum=mod_dis if i_s1>=i_max-1: i_s1=i_max-1 #for i in range(0,i_max-1,1): for i in range(i_s0,i_s1+1,1): #adding new material radius[i,j-1]=radius[i_s0-1,j-1]+(dr_max*(i-(i_s0-1))) rho[i,j-1]=rho_0 rho_sil[i,j-1]=rho_0 phi[i,j-1]=phi_0 T[i,j-1]=T0 #gravity profile for i in range(0,i_s1+1,1): r[i]=radius[i,j-1] rho_old2[i]=rho[i,j-1] if i==0: g[i,j]=(bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max) else: g[i,j]=(((r[i-1]/r[i])**2)*g[i-1,j])+((bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max)) for i in range((-i_max+i_s1),-(i_max+1),-1): if i==(-i_max+i_s1): P[i,j]=0 #surface pressure =0 else: dP[i]=(rho[i,j-1]*g[i,j]*dr_max) P[i,j]=P[i+1,j]+dP[i] #compaction for i in range(0,i_s1+1,1): phi_old=phi[i,j-1] #old value of phi g_phi=((1-phi_0)/(1-phi_old))**(2/3) f_phi=0.5*((3/(np.pi*((2*(g_phi**(1/2))*(3-g_phi))-3)))**(1/3)) P_grain=P[i,j]*((np.pi)/(2*(3**0.5)*((4*(3**0.5)*((1-phi_old)**(2/3))*f_phi*f_phi)-1))) dln_phi=A*(P_grain**(2/3))*(b**(-3))*np.exp(-Ea_comp/(Rgas*T[i,j-1]))*dt phi_old=phi[i,j-1] ln_phi_old=np.log(1-phi_old) ln_phi=ln_phi_old+dln_phi phi[i,j]=1-np.exp(ln_phi) if phi[i,j]<=0: phi[i,j]=0 rho_sil[i,j]=rhom*(1-phi[i,j]) k[i,j]=kb*((np.exp((-4*phi[i,j])/phi_1)+np.exp(-4.4-((4*phi[i,j])/phi_2)))**(1/4)) #Krause et al, 2011, LPSC expression k[i,j]=kb*np.exp(100*phi[i,j]*(-0.1246)) rho[i,j]=((chi_fe_tot[i,j-1]/rho_fe)+((1-chi_fe_tot[i,j-1])/rho_sil[i,j]))**(-1) for i in range(0,i_s1+1,1): #over entire body including surface node - need to get planetary radius if i==0: radius[i,j]=((rho[i,j-1]/rho[i,j])**(1/3))*radius[i,j-1] dr[i,j]=(radius[i,j]) else: Vold=(4/3)*np.pi*((radius[i,j-1]**3)-(radius[i-1,j-1]**3)) Vbelow=(4/3)*np.pi*(radius[i-1,j]**3) radius[i,j]=((3/(4*np.pi))*(((rho[i,j-1]/rho[i,j])*Vold)+Vbelow))**(1/3) dr[i,j]=(radius[i,j]-radius[i-1,j]) #compacted radius of asteroid Rp_c[j]=radius[i_s1-1,j] T[i_s1,j]=T0 #set surface temp eta[i_s1,j]=4E30 for i in range(0,i_s1,1): #thermal evolution bit h=h0*c0*np.exp((-1)*np.log(2)*t/thalf) if i==0: dTdr=0 d2Tdr2=(T[i+1,j-1]-T[i,j-1])/((radius[i+1,j-1])**2) dkdr=(k[i+1,j]-k[i,j])/(radius[i+1,j]-radius[i,j]) else: dTdr=(T[i+1,j-1]-T[i-1,j-1])/(radius[i+1,j-1]-radius[i-1,j-1]) d2Tdr2=(1/dr[i,j-1])*((T[i+1,j-1]-T[i,j-1])/(radius[i+1,j-1]-radius[i,j-1])-(T[i,j-1]-T[i-1,j-1])/(radius[i,j-1]-radius[i-1,j-1])) dkdr=(k[i+1,j]-k[i-1,j])/(radius[i+1,j]-radius[i-1,j]) E_cond=((dkdr*dTdr)+(k[i,j]*((((1/radius[i,j-1])*dTdr)+d2Tdr2)))) #contribution from conduction if T[i,j-1]=0.9999: chi[i,j]=1 if T[i,j]>Tl_fe: chi[i,j]=1 chi_fe[i,j]=x0_fe*chi[i,j] if T[i,j]>=Tl_fe0: T[i,j]=Tl_fe0 #viscosity stuff #silicate melt fraction - assuming linear melting - will need to modify the specific heat capacity - should probably add into cfl check if T[i,j]ceu: frac_S[i,j]=ceu #looking for convection con_check=0 #i_con_old=i_con i_con=-1 i=i_s1 convecting[i_s1,j]=-5 eta[i_s1,j]=1E30 while i_con<=0: i=i-1 critra[i,j]=20.9*((gamma*(T[i,j]-1800))**(4)) rayleigh[i,j]=(g[i,j]*a_sil_exp*rho[i,j]*rho[i,j]*(T[0,j]-T[i,j])*(radius[i,j]**3)*cpm)/(eta[i,j]*k[i,j]) if rayleigh[i,j]>critra[i,j]: i_con=i if i_con>=i_s1: i_con=i_s1-1 break elif i<=1: break #mass of liquid iron present in body if i_con>=0: i=-1 m_l_fe[j]=0 mcoretot=0 mscore=0 dTs=0 mcontot=0 T_m_mixed=0 while i<=i_con: i=i+1 if i==0: vnode=(4/3)*np.pi*(radius[i,j]**3) else: vnode=(4/3)*np.pi*((radius[i,j]**3)-(radius[i-1,j]**3)) mnode=rho[i,j]*vnode mfenode=chi_fe[i,j]*mnode msnode=frac_S[i,j]*mnode*chi_fe[i,j] dTsnode=mfenode*T[i,j] mcoretot=mcoretot+mfenode mscore=mscore+msnode dTs=dTs+dTsnode msilnode=mnode-mfenode mcontot=mcontot+msilnode dT_m=msilnode*T[i,j] T_m_mixed=T_m_mixed+dT_m if i>=i_con: break m_l_fe[j]=mcoretot frac_S_core[j]=mscore/m_l_fe[j] Xs_at=frac_S_core[j]/(frac_S_core[j]+((m_mol_S/m_mol_Fe)*(1-frac_S_core[j]))) Xs_at=Xs_at/100 mcontotold=mcontot Tmixed=T_m_mixed/mcontot rcore[j]=((3*m_l_fe[j])/(4*np.pi*rho_l_fe_0))**(1/3) Tpeak[j]=np.amax(T[:,j]) i_s0=i_s1 #if j>=500: if i_s1>=(i_max-1): j_fullsize=j break if i_con>=0: j_con=j j_congrow=j i_core_old=0 if i_con==0: dr_ocean=radius[i_con,j] else: dr_ocean=radius[i_con,j]/(i_con+1) dr_con=dr_ocean db=dr_ocean i_core=int(rcore[j]/dr_con) Tdiff=T[i_con,j] Sdiff=frac_S_core[j] i_con_old=0 TM[j]=Tmixed TC[j]=Tmixed T[:i_con,j]=Tmixed if i_core>=1: TCMB[j_con]=(TM[j_con]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j_con]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) else: TCMB[j_con]=TM[j_con] eta_m=eta[i_con,j] eta_m_old=eta_m if i_core>=1: db=radius[i_core-1,j] else: db=0 exp_data[6,0]=time[j] exp_data[7,0]=Tdiff exp_data[8,0]=Sdiff exp_data[9,0]=1810-(18*Sdiff) exp_data[31,0]=radius[i_s1,j] break if time[j]>6: if np.all(T[:,j]<1600): exp_data[17,0]=time[j] exp_data[19,0]=np.amax(T[:,:j]) exp_data[18,0]=np.amax(chi_fe[:,:j]) j_undiff=j print('does not differentiate') break if Rp_u>=Rf: j_fullsize=j print(i_s1,i_max-1) break drc_sum=0 i_core_old=0 jcount=0 E_f1=0 TM[j_end]=T[i_con,j_end] if j_fullsize<0: #convection starts before fully accreted while j0: break j=j+1 modj=j%1000 t=t0+(dt*j) time[j]=t/(1E6*365*24*60*60) #increase in uncompacted radius this time step dRp_u=(R0/dt_acc)*np.log(Rf/R0)*dt*np.exp(np.log(Rf/R0)*(1/dt_acc)*(t-t0)) dRp_u_sum=dRp_u+dRp_u_sum Rp_u=Rp_u+dRp_u #how many nodes do we add this time step? discriminant=dRp_u_sum-dr_max if discriminant<0: #no nodes added dRp_u_sum=dRp_u_sum i_s1=i_s0 mod_dis=discriminant%dr_max if mod_dis==0: #exact integer number of nodes added - no remainder to carry over i_s1=i_s0+int(dRp_u_sum//dr_max) dRp_u_sum=0 if mod_dis>0: i_s1=i_s0+int(dRp_u_sum//dr_max) dRp_u_sum=mod_dis if i_s1>=i_max-1: i_s1=i_max-1 #for i in range(0,i_max-1,1): for i in range(i_s0,i_s1+1,1): #adding new material radius[i,j-1]=radius[i_s0-1,j-1]+(dr_max*(i-(i_s0-1))) rho[i,j-1]=rho_0 rho_sil[i,j-1]=rho_0 phi[i,j-1]=phi_0 T[i,j-1]=T0 core_dis=rcore[j-1]-dr_con if core_dis<0: drc_sum=drc_sum i_core=i_core_old mod_core=core_dis%dr_con if mod_core==0: i_core=int(rcore[j-1]//dr_con) drc_sum=0 if mod_core>0: i_core=int(rcore[j-1]//dr_con) drc_sum=mod_core for i in range(0,i_core,1): if i_core==i_core_old: rho[i,j]=rho[i,j-1] #TC[j]=TC[j-1] core_T=TC[j] else: #TC[j]=TC[j-1] core_T=TC[j] if ii_core_old: for i in range(i_core_old,i_core,1): T[i,j-1]=TC[j-1] if i_core<1: TC[j]=TM[j-1] chi_fe[:i_core,j-1]=1 chi_fe_tot[:i_core,j-1]=1 k[:i_core,j]=k_fe phi[:i_core,j]=0 rho_sil[:i_core,j]=0 eta[:i_core,j]=eta_fe for i in range(0,i_con,1): radius[i,j]=dr_con*float(i+1) #changing magma ocean properties chi_fe[i_core:i_con,j-1]=0 chi_fe_tot[i_core:i_con,j]=0 rho[i_core:i_con,j]=rhom rho_sil[i_core:i_con,j]=rhom k[i_core:i_con,j]=kb phi[i_core:i_con,j]=0 T_c_old=TC[j-1] T_m_old=TM[j-1] T_cmb_old=TCMB[j-1] #gravity profile for i in range(0,i_s1+1,1): r[i]=radius[i,j-1] rho_old2[i]=rho[i,j-1] if i==0: g[i,j]=(bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max) else: g[i,j]=(((r[i-1]/r[i])**2)*g[i-1,j])+((bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max)) for i in range((-i_max+i_s1),-(i_max+1),-1): if i==(-i_max+i_s1): P[i,j]=0 #surface pressure =0 else: dP[i]=(rho[i,j-1]*g[i,j]*dr_max) P[i,j]=P[i+1,j]+dP[i] #compaction - only occurring in lid portion for i in range(i_con,i_s1+1,1): if phi[i,j-1]>0: phi_old=phi[i,j-1] #old value of phi g_phi=((1-phi_0)/(1-phi_old))**(2/3) f_phi=0.5*((3/(np.pi*((2*(g_phi**(1/2))*(3-g_phi))-3)))**(1/3)) P_grain=P[i,j]*((np.pi)/(2*(3**0.5)*((4*(3**0.5)*((1-phi_old)**(2/3))*f_phi*f_phi)-1))) dln_phi=A*(P_grain**(2/3))*(b**(-3))*np.exp(-Ea_comp/(Rgas*T[i,j-1]))*dt phi_old=phi[i,j-1] ln_phi_old=np.log(1-phi_old) ln_phi=ln_phi_old+dln_phi phi[i,j]=1-np.exp(ln_phi) if phi[i,j]<=0: phi[i,j]=0 else: phi[i,j]=phi[i,j-1] rho_sil[i,j]=rhom*(1-phi[i,j]) k[i,j]=kb*((np.exp((-4*phi[i,j])/phi_1)+np.exp(-4.4-((4*phi[i,j])/phi_2)))**(1/4)) #Krause et al, 2011, LPSC expression k[i,j]=kb*np.exp(100*phi[i,j]*(-0.1246)) rho[i,j]=((chi_fe_tot[i,j-1]/rho_fe)+((1-chi_fe_tot[i,j-1])/rho_sil[i,j]))**(-1) for i in range(i_con,i_s1+1,1): #now just over compacting area if i==0: radius[i,j]=((rho[i,j-1]/rho[i,j])**(1/3))*radius[i,j-1] dr[i,j]=(radius[i,j]) else: Vold=(4/3)*np.pi*((radius[i,j-1]**3)-(radius[i-1,j-1]**3)) Vbelow=(4/3)*np.pi*(radius[i-1,j]**3) radius[i,j]=((3/(4*np.pi))*(((rho[i,j-1]/rho[i,j])*Vold)+Vbelow))**(1/3) dr[i,j]=(radius[i,j]-radius[i-1,j]) #compacted radius of asteroid Rp_c[j]=radius[i_s1-1,j] T[i_s1,j]=T0 #set surface temp eta[i_s1,j]=4E30 i_s0=i_s1 #for next time step #thermal evolution #first of lid for i in range(i_con,i_s1,1): h=h0*c0*np.exp((-1)*np.log(2)*t/thalf) if i==0: dTdr=0 d2Tdr2=(T[i+1,j-1]-T[i,j-1])/((radius[i+1,j-1])**2) dkdr=(k[i+1,j]-k[i,j])/(radius[i+1,j]-radius[i,j]) else: dTdr=(T[i+1,j-1]-T[i-1,j-1])/(radius[i+1,j-1]-radius[i-1,j-1]) d2Tdr2=(1/dr[i,j-1])*((T[i+1,j-1]-T[i,j-1])/(radius[i+1,j-1]-radius[i,j-1])-(T[i,j-1]-T[i-1,j-1])/(radius[i,j-1]-radius[i-1,j-1])) dkdr=(k[i+1,j]-k[i-1,j])/(radius[i+1,j]-radius[i-1,j]) E_cond=((dkdr*dTdr)+(k[i,j]*((((1/radius[i,j-1])*dTdr)+d2Tdr2)))) #contribution from conduction if T[i,j-1]=0.9999: chi[i,j]=1 if T[i,j]>Tl_fe: chi[i,j]=1 chi_fe[i,j]=x0_fe*chi[i,j] if T[i,j]>=Tl_fe0: T[i,j]=Tl_fe0 #viscosity stuff #silicate melt fraction - assuming linear melting - will need to modify the specific heat capacity - should probably add into cfl check if T[i,j]ceu: frac_S[i,j]=ceu #magma ocean evolution h=h0*c0*((mcontotold+mcoretot)/(mcontotold))*np.exp((-1)*np.log(2)*t/thalf) racrit=1000 d0=(racrit**(1/3))*(((gamma*(TM[j-1]-T0))/8)**(4/3))*(((kb*eta_m_old)/(rhom*rhom*2*cpm*a_sil_exp*g[i_con,j]*(TM[j-1]-T0)))**(1/3)) fs=kb*((TM[j-1]-T0)/d0) if T_m_old>=T_c_old: #conductive fluxes, passing heat into the core f2=kb*(T_cmb_old-T_m_old)/dr_ocean f1=k_fe*(T_c_old-T_cmb_old)/dr_ocean elif b_int>=0: f1=k_fe*((T_c_old-T_cmb_old)**(4/3))*(((g[i_core,j]*rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/3)) f2=kb*((T_cmb_old-T_m_old)**(4/3))*(((g[i_core,j]*rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/3)) else: f2=kb*(T_cmb_old-T_m_old)/dr_ocean f1=k_fe*(T_c_old-T_cmb_old)/dr_ocean rm=radius[i_con,j] surfflux[j]=fs Alid=4*np.pi*(rm**2) Asurf=4*np.pi*(radius[i_s1-1,j]**2) if i_core==0: Vm=4/3*np.pi*(rm**3) Acmb=0 else: Vm=4/3*np.pi*((rm**3)-(radius[i_core-1,j]**3)) RADC=radius[i_core,j] Acmb=4*np.pi*(RADC**2) E_change=((-fs*Alid)+(h*rhom*Vm)+(f2*Acmb)) F1[j]=f1 F2[j]=f2 if TM[j-1]>=T_sil_s: cpm=2*850 else: cpm=850 dT=(E_change*dt)/(Vm*rhom*cpm) TM[j]=TM[j-1]+dT if TM[j]=0: Vc=((4/3)*np.pi)*((radius[i_core,j]**3)-(Rb**3)) delTC=(dt/(rho_l_fe_0*cpc*Vc))*((-f1)*Acmb) TC[j]=delTC+TC[j-1] for i in range(b_int,i_core,1): T[i,j]=TC[j] chi_sil[i,j]=0 #core - diffusive for i in range(0,b_int,1): if i==0: dTdr=0 d2Tdr2=(T[i+1,j-1]-T[i,j-1])/((radius[i+1,j-1])**2) else: dTdr=(T[i+1,j-1]-T[i-1,j-1])/(radius[i+1,j-1]-radius[i-1,j-1]) d2Tdr2=(1/dr_ocean)*((T[i+1,j-1]-T[i,j-1])/(radius[i+1,j-1]-radius[i,j-1])-(T[i,j-1]-T[i-1,j-1])/(radius[i,j-1]-radius[i-1,j-1])) T[i,j]=((k[i,j]/(rho_l_fe_0*cpc))*dt*((((1/radius[i,j-1])*dTdr)+d2Tdr2)))+T[i,j-1] chi_sil[i,j]=0 else: #core - diffusive for i in range(0,i_core,1): if i==0: dTdr=0 d2Tdr2=(T[i+1,j-1]-T[i,j-1])/((radius[i+1,j-1])**2) else: dTdr=(T[i+1,j-1]-T[i-1,j-1])/(radius[i+1,j-1]-radius[i-1,j-1]) d2Tdr2=(1/dr_ocean)*((T[i+1,j-1]-T[i,j-1])/(radius[i+1,j-1]-radius[i,j-1])-(T[i,j-1]-T[i-1,j-1])/(radius[i,j-1]-radius[i-1,j-1])) T[i,j]=((k[i,j]/(rho_l_fe_0*cpc))*dt*((((1/radius[i,j-1])*dTdr)+d2Tdr2)))+T[i,j-1] chi_sil[i,j]=0 radius_convecting[j]=radius[i_con,j] TC[j]=T[i_core-1,j] #for rayleigh number check i_old_2=i_con_old #looking for convection con_check=0 i_con_old=i_con i_con=-1 i=i_s1 convecting[i_s1,j]=-5 eta[i_s1,j]=1E30 while i_con<=0: i=i-1 critra[i,j]=20.9*((gamma*(T[i,j]-1800))**(4)) rayleigh[i,j]=(g[i,j]*a_sil_exp*rho[i,j]*rho[i,j]*(TM[j]-T[i,j])*(radius[i,j]**3)*cpm)/(eta[i,j]*k[i,j]) if rayleigh[i,j]>critra[i,j]: i_con=i if i_con>=i_s1: i_con=i_s1-1 break elif i<=i_con_old: i_con=i_con_old break dRa_check=rayleigh[i_con,j]-rayleigh[i_old_2,j-1] if T[i_con,j]>T[i_con,j-1]: dRa_check=0 #mass of core added and new core radius if i_con-i_con_old>0: i=i_con_old-1 m_l_fe[j]=0 mcoretot_old=0 mcoretot=0 mscore=0 dTs=0 mcontot=0 T_m_mixed=0 while i<=i_con: i=i+1 if i==0: vnode=(4/3)*np.pi*(radius[i,j]**3) else: vnode=(4/3)*np.pi*((radius[i,j]**3)-(radius[i-1,j]**3)) mnode=rho[i,j]*vnode mfenode=chi_fe[i,j]*mnode msilnode=mnode-mfenode dTsnode=mfenode*TM[j] msnode=chi_fe[i,j]*frac_S[i,j]*mnode mcoretot=mcoretot+mfenode mscore=mscore+msnode dTs=dTs+dTsnode mcontot=mcontot+msilnode dTm=msilnode*T[i,j] T_m_mixed=T_m_mixed+dTm if i>=i_con: break Tmixed=((mcontotold*TM[j])+T_m_mixed)/(mcontot+mcontotold) TM[j]=Tmixed m_l_fe[j]=mcoretot+m_l_fe[j-1] frac_S_core[j]=((frac_S_core[j-1]*m_l_fe[j-1])+mscore)/m_l_fe[j] frac_S_new=mscore/mcoretot Xs_at=frac_S_new/(frac_S_new+((m_mol_S/m_mol_Fe)*(1-frac_S_new))) Xs_at=Xs_at/100 mcontotold=mcontot+mcontotold else: m_l_fe[j]=m_l_fe[j-1] frac_S_core[j]=frac_S_core[j-1] mcontotold=mcontotold if m_l_fe[j]==0: rcore[j]=0 else: rcore[j]=((3*m_l_fe[j])/(4*np.pi*rho_l_fe_0))**(1/3) frac_S[:i_core,j]=frac_S_core[j] i_core_old=i_core Tpeak[j]=np.amax(T[:,j]) #new bit added - 29/05-works out TCMB before core material added if i_core==0: TC[j]=TM[j] if TM[j]>=TC[j]: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) elif b_int>=0: const_TC=(k_fe**(3/4))*(((rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/4)) const_TM=(kb**(3/4))*(((rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/4)) TCMB[j]=((const_TC*TC[j])+(const_TM*TM[j]))/(const_TC+const_TM) else: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) eta_cmb_old=10**(a0+(a1*TCMB[j])+a2*(np.tanh((TCMB[j]+a3)/a4))) #what does core look like after addition of new material? if i_core>=1: if i_con>i_con_old: #only matters if adding new material to the core dT_core[:i_core,j]=TM[j]-T[:i_core,j] if np.any(dT_core[:,j]<0): b_stable=np.transpose(np.asarray(np.where(dT_core[:,j]<0))) b_int=np.asscalar(b_stable[0]) Rb=radius[b_int,j] db=radius[i_core,j]-Rb #starting volume average to get temp of convecting layer TSUM=0 i=b_int-1 mcoremix=0 while i=i_core-1: break Tc_mixed=(TSUM+(TCMB[j-1]*mcoretot))/(mcoremix+mcoretot) for i in range(b_int,i_core,1): T[i,j]=Tc_mixed TC[j]=Tc_mixed else: b_int=-1 TC[j]=TM[j] #general instability in the core if i_core>=1: if i_con==i_con_old: unstable[:i_core-1,j]=T[i_core-1,j]-T[:i_core-1,j] if np.any(unstable[:,j]<0): b_stable=np.transpose(np.asarray(np.where(unstable[:,j]<0))) b_int=np.asscalar(b_stable[0]) Rb=radius[b_int,j] db=radius[i_core,j]-Rb #starting volume average to get temp of convecting layer TSUM=0 i=b_int-1 mcoremix=0 while i=i_core-1: break Tc_mixed=TSUM/mcoremix for i in range(b_int,i_core,1): T[i,j]=Tc_mixed TC[j]=Tc_mixed junstable=j else: b_int=-1 TC[j]=TC[j] dcon[j]=db #new bit added - 29/05 if i_core==0: TC[j]=TM[j] if TM[j]>=TC[j]: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) elif b_int>=0: const_TC=(k_fe**(3/4))*(((rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/4)) const_TM=(kb**(3/4))*(((rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/4)) TCMB[j]=((const_TC*TC[j])+(const_TM*TM[j]))/(const_TC+const_TM) else: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) eta_cmb_old=10**(a0+(a1*TCMB[j])+a2*(np.tanh((TCMB[j]+a3)/a4))) f_ad=(k_fe*a_exp*T[i_core-1,j]*g[i_core-1,j])/cpc if i_con==i_con_old: E1_ave=(f1*dt)+E1_ave else: E1_ave=(f1*dt)+E1_ave x=x+1 j_ave=int(round((j+j_congrow)/2)) time_ave[x]=time[j_ave] dt_ave=(time[j]-time[j_congrow])*(1E6*365*24*3600) if dt_avef_ad: fdrive=f1_ave[x]-f_ad u_mac=((t_spin*bigg*a_exp*RADC*fdrive)/(cpc))**(1/2) mag_rey_ave[x]=(u_mac*db)/mag_diff else: mag_rey_ave[x]=0 if f1>f_ad: fdrive=f1-f_ad u_mac=((t_spin*bigg*a_exp*RADC*fdrive)/(cpc))**(1/2) mag_rey[j]=(u_mac*db)/mag_diff else: mag_rey[j]=0 i_s0=i_s1 if i_s1>=(i_max-1): #only breaking when fully accreted. j_end=j j_fullsize=j break dm=rm-RADC rayleigh_mantle=(g[i_con,j]*a_sil_exp*(TM[j]-T0)*rhom*(dm**3))/(diffm*eta_m) if rayleigh_mantle<1000: j_cond=j exp_data[15,0]=time[j] break if Rp_u>=Rf: j_end=j j_fullsize=j break #else is fully accreted before onset of convection else: while j0: break j=j+1 modj=j%1000 t=t0+(dt*j) time[j]=t/(1E6*365*24*60*60) #i_s1=i_max-1 #gravity profile for i in range(0,i_s1+1,1): r[i]=radius[i,j-1] rho_old2[i]=rho[i,j-1] if i==0: g[i,j]=(bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max) else: g[i,j]=(((r[i-1]/r[i])**2)*g[i-1,j])+((bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max)) for i in range((-i_max+i_s1),-(i_max+1),-1): if i==(-i_max+i_s1): P[i,j]=0 #surface pressure =0 else: dP[i]=(rho[i,j-1]*g[i,j]*dr_max) P[i,j]=P[i+1,j]+dP[i] #compaction for i in range(0,i_s1+1,1): phi_old=phi[i,j-1] #old value of phi g_phi=((1-phi_0)/(1-phi_old))**(2/3) f_phi=0.5*((3/(np.pi*((2*(g_phi**(1/2))*(3-g_phi))-3)))**(1/3)) P_grain=P[i,j]*((np.pi)/(2*(3**0.5)*((4*(3**0.5)*((1-phi_old)**(2/3))*f_phi*f_phi)-1))) dln_phi=A*(P_grain**(2/3))*(b**(-3))*np.exp(-Ea_comp/(Rgas*T[i,j-1]))*dt phi_old=phi[i,j-1] ln_phi_old=np.log(1-phi_old) ln_phi=ln_phi_old+dln_phi phi[i,j]=1-np.exp(ln_phi) if phi[i,j]<=0: phi[i,j]=0 rho_sil[i,j]=rhom*(1-phi[i,j]) k[i,j]=kb*((np.exp((-4*phi[i,j])/phi_1)+np.exp(-4.4-((4*phi[i,j])/phi_2)))**(1/4)) #Krause et al, 2011, LPSC expression k[i,j]=kb*np.exp(100*phi[i,j]*(-0.1246)) rho[i,j]=((chi_fe_tot[i,j-1]/rho_fe)+((1-chi_fe_tot[i,j-1])/rho_sil[i,j]))**(-1) for i in range(0,i_s1+1,1): #over entire body including surface node - need to get planetary radius if i==0: radius[i,j]=((rho[i,j-1]/rho[i,j])**(1/3))*radius[i,j-1] dr[i,j]=(radius[i,j]) else: Vold=(4/3)*np.pi*((radius[i,j-1]**3)-(radius[i-1,j-1]**3)) Vbelow=(4/3)*np.pi*(radius[i-1,j]**3) radius[i,j]=((3/(4*np.pi))*(((rho[i,j-1]/rho[i,j])*Vold)+Vbelow))**(1/3) dr[i,j]=(radius[i,j]-radius[i-1,j]) #compacted radius of asteroid Rp_c[j]=radius[i_s1-1,j] T[i_s1,j]=T0 #set surface temp eta[i_s1,j]=4E30 i_s0=i_s1 #for next time step for i in range(0,i_s1,1): h=h0*c0*np.exp((-1)*np.log(2)*t/thalf) if i==0: dTdr=0 d2Tdr2=(T[i+1,j-1]-T[i,j-1])/((radius[i+1,j-1])**2) dkdr=(k[i+1,j]-k[i,j])/(radius[i+1,j]-radius[i,j]) else: dTdr=(T[i+1,j-1]-T[i-1,j-1])/(radius[i+1,j-1]-radius[i-1,j-1]) d2Tdr2=(1/dr[i,j-1])*((T[i+1,j-1]-T[i,j-1])/(radius[i+1,j-1]-radius[i,j-1])-(T[i,j-1]-T[i-1,j-1])/(radius[i,j-1]-radius[i-1,j-1])) dkdr=(k[i+1,j]-k[i-1,j])/(radius[i+1,j]-radius[i-1,j]) E_cond=((dkdr*dTdr)+(k[i,j]*((((1/radius[i,j-1])*dTdr)+d2Tdr2)))) #contribution from conduction if T[i,j-1]=0.9999: chi[i,j]=1 if T[i,j]>Tl_fe: chi[i,j]=1 chi_fe[i,j]=x0_fe*chi[i,j] if T[i,j]>=Tl_fe0: T[i,j]=Tl_fe0 #viscosity stuff #silicate melt fraction - assuming linear melting - will need to modify the specific heat capacity - should probably add into cfl check if T[i,j]ceu: frac_S[i,j]=ceu #looking for convection con_check=0 #i_con_old=i_con i_con=-1 i=i_s1+1 convecting[i_s1,j]=-5 eta[i_s1,j]=1E30 while i_con<=0: i=i-1 critra[i,j]=20.9*((gamma*(T[i,j]-1800))**(4)) rayleigh[i,j]=(g[i,j]*a_sil_exp*rho[i,j]*rho[i,j]*(T[0,j]-T[i,j])*(radius[i,j]**3)*cpm)/(eta[i,j]*k[i,j]) if rayleigh[i,j]>critra[i,j]: i_con=i if i_con>=i_s1: i_con=i_s1-1 break elif i<=1: break if i_con>=0: #mass of liquid iron present in body i=-1 m_l_fe[j]=0 mcoretot=0 mscore=0 dTs=0 T_m_mixed=0 mcontot=0 while i<=i_con: i=i+1 if i==0: vnode=(4/3)*np.pi*(radius[i,j]**3) else: vnode=(4/3)*np.pi*((radius[i,j]**3)-(radius[i-1,j]**3)) mnode=rho[i,j]*vnode mfenode=chi_fe[i,j]*mnode dTsnode=mfenode*T[i,j] msnode=chi_fe[i,j]*frac_S[i,j]*mnode mcoretot=mcoretot+mfenode mscore=mscore+msnode dTs=dTs+dTsnode #m_l_fe[j]=m_l_fe[j]+mfenode msilnode=mnode-mfenode mcontot=msilnode+mcontot dT_m=msilnode*T[i,j] T_m_mixed=T_m_mixed+dT_m if i>=i_con: break m_l_fe[j]=mcoretot Tmixed=T_m_mixed/mcontot mcontotold=mcontot if m_l_fe[j]==0: rcore[j]=0 frac_S_core[j]=0 else: rcore[j]=((3*m_l_fe[j])/(4*np.pi*rho_l_fe_0))**(1/3) frac_S_core[j]=mscore/m_l_fe[j] frac_S_new=mscore/mcoretot Xs_at=frac_S_new/(frac_S_new+((m_mol_S/m_mol_Fe)*(1-frac_S_new))) Xs_at=Xs_at/100 Tpeak[j]=np.amax(T[:,j]) if i_con>=0: j_con=j j_congrow=j i_core_old=0 if i_con==0: dr_ocean=radius[i_con,j] Tmixed=TM[j] else: dr_ocean=radius[i_con,j]/(i_con+1) if i_con<=i_core: i_con=i_core+1 dr_con=dr_ocean i_core=int(rcore[j]/dr_con) dr_con=dr_ocean db=dr_ocean i_con_old=0 TC[j]=Tmixed TM[j]=Tmixed T[:i_con,j]=Tmixed if i_core>=1: TCMB[j_con]=(TM[j_con]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j_con]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) else: TCMB[j_con]=TM[j_con] Tdiff=T[i_con,j] Sdiff=frac_S_core[j] exp_data[6,0]=time[j] exp_data[7,0]=Tdiff exp_data[8,0]=Sdiff exp_data[9,0]=1810-(18*Sdiff) exp_data[31,0]=radius[i_s1,j] break elif time[j]>6: if np.all(T[:,j]<1600): exp_data[17,0]=time[j] exp_data[19,0]=np.amax(T[:,:j]) exp_data[18,0]=np.amax(chi_fe[:,:j]) j_undiff=j print('does not differentiate') break else: continue if j_undiff>0: pass elif j_cond>0: pass else: eta_m=eta[i_con,j] eta_m_old=eta_m while j<=j_max: j=j+1 modj=j%1000 t=t0+(dt*j) time[j]=t/(1E6*365*24*60*60) core_dis=rcore[j-1]-dr_con if core_dis<0: drc_sum=drc_sum i_core=i_core_old mod_core=core_dis%dr_con if mod_core==0: i_core=int(rcore[j-1]//dr_con) drc_sum=0 if mod_core>0: i_core=int(rcore[j-1]//dr_con) drc_sum=mod_core #setting up 'new' thermal structure in core for i in range(0,i_core,1): if i_core==i_core_old: rho[i,j]=rho[i,j-1] else: if ii_core_old: for i in range(i_core_old,i_core,1): T[i,j-1]=TC[j-1] if i_core<1: TC[j]=TM[j-1] #new stuff T_c_old=TC[j-1] T_m_old=TM[j-1] T_cmb_old=TCMB[j-1] chi_fe[:i_core,j-1]=1 chi_fe_tot[:i_core,j-1]=1 k[:i_core,j]=k_fe phi[:i_core,j]=0 rho_sil[:i_core,j]=0 eta[:i_core,j]=eta_fe for i in range(0,i_con,1): radius[i,j]=dr_con*float(i+1) #changing magma ocean properties chi_fe[i_core:i_con,j-1]=0 chi_fe_tot[i_core:i_con,j]=0 rho[i_core:i_con,j]=rhom rho_sil[i_core:i_con,j]=rhom k[i_core:i_con,j]=kb phi[i_core:i_con,j]=0 #gravity profile for i in range(0,i_s1+1,1): r[i]=radius[i,j-1] rho_old2[i]=rho[i,j-1] if i==0: g[i,j]=(bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max) else: g[i,j]=(((r[i-1]/r[i])**2)*g[i-1,j])+((bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max)) #pressure profile for i in range((-i_max+i_s1),-(i_max+1),-1): if i==(-i_max+i_s1): P[i,j]=0 #surface pressure =0 else: dP[i]=(rho[i,j-1]*g[i,j]*dr_max) P[i,j]=P[i+1,j]+dP[i] #compaction - only occurring in lid portion for i in range(i_con,i_s1+1,1): if phi[i,j-1]>0: phi_old=phi[i,j-1] #old value of phi g_phi=((1-phi_0)/(1-phi_old))**(2/3) f_phi=0.5*((3/(np.pi*((2*(g_phi**(1/2))*(3-g_phi))-3)))**(1/3)) P_grain=P[i,j]*((np.pi)/(2*(3**0.5)*((4*(3**0.5)*((1-phi_old)**(2/3))*f_phi*f_phi)-1))) dln_phi=A*(P_grain**(2/3))*(b**(-3))*np.exp(-Ea_comp/(Rgas*T[i,j-1]))*dt phi_old=phi[i,j-1] ln_phi_old=np.log(1-phi_old) ln_phi=ln_phi_old+dln_phi phi[i,j]=1-np.exp(ln_phi) if phi[i,j]<=0: phi[i,j]=0 else: phi[i,j]=phi[i,j-1] rho_sil[i,j]=rhom*(1-phi[i,j]) k[i,j]=kb*((np.exp((-4*phi[i,j])/phi_1)+np.exp(-4.4-((4*phi[i,j])/phi_2)))**(1/4)) #Krause et al, 2011, LPSC expression k[i,j]=kb*np.exp(100*phi[i,j]*(-0.1246)) rho[i,j]=((chi_fe_tot[i,j-1]/rho_fe)+((1-chi_fe_tot[i,j-1])/rho_sil[i,j]))**(-1) for i in range(i_con,i_s1+1,1): #now just over the top compacting bit. if i==0: radius[i,j]=((rho[i,j-1]/rho[i,j])**(1/3))*radius[i,j-1] dr[i,j]=(radius[i,j]) else: Vold=(4/3)*np.pi*((radius[i,j-1]**3)-(radius[i-1,j-1]**3)) Vbelow=(4/3)*np.pi*(radius[i-1,j]**3) radius[i,j]=((3/(4*np.pi))*(((rho[i,j-1]/rho[i,j])*Vold)+Vbelow))**(1/3) dr[i,j]=(radius[i,j]-radius[i-1,j]) Rp_c[j]=radius[i_s1,j] T[i_s1,j]=T0 #thermal evolution #first in the conductive lid for i in range(i_con,i_s1,1): h=h0*c0*np.exp((-1)*np.log(2)*t/thalf) if i==0: dTdr=0 d2Tdr2=(T[i+1,j-1]-T[i,j-1])/((radius[i+1,j-1])**2) dkdr=(k[i+1,j]-k[i,j])/(radius[i+1,j]-radius[i,j]) else: dTdr=(T[i+1,j-1]-T[i-1,j-1])/(radius[i+1,j-1]-radius[i-1,j-1]) d2Tdr2=(1/dr[i,j-1])*((T[i+1,j-1]-T[i,j-1])/(radius[i+1,j-1]-radius[i,j-1])-(T[i,j-1]-T[i-1,j-1])/(radius[i,j-1]-radius[i-1,j-1])) dkdr=(k[i+1,j]-k[i-1,j])/(radius[i+1,j]-radius[i-1,j]) E_cond=((dkdr*dTdr)+(k[i,j]*((((1/radius[i,j-1])*dTdr)+d2Tdr2)))) #contribution from conduction if T[i,j-1]=0.9999: chi[i,j]=1 if T[i,j]>Tl_fe: chi[i,j]=1 chi_fe[i,j]=x0_fe*chi[i,j] if T[i,j]>=Tl_fe0: T[i,j]=Tl_fe0 #viscosity stuff #silicate melt fraction - assuming linear melting - will need to modify the specific heat capacity - should probably add into cfl check if T[i,j]ceu: frac_S[i,j]=ceu #magma ocean evolution h=h0*c0*((mcontotold+mcoretot)/(mcontotold))*np.exp((-1)*np.log(2)*t/thalf) racrit=1000 d0=(racrit**(1/3))*(((gamma*(TM[j-1]-T0))/8)**(4/3))*(((kb*eta_m_old)/(rhom*rhom*2*cpm*a_sil_exp*g[i_con,j]*(TM[j-1]-T0)))**(1/3)) fs=kb*((TM[j-1]-T0)/d0) if T_m_old>=T_c_old: #conductive fluxes, passing heat into the core f2=kb*(T_cmb_old-T_m_old)/dr_ocean f1=k_fe*(T_c_old-T_cmb_old)/dr_ocean elif b_int>=0: f1=k_fe*((T_c_old-T_cmb_old)**(4/3))*(((g[i_core,j]*rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/3)) f2=kb*((T_cmb_old-T_m_old)**(4/3))*(((g[i_core,j]*rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/3)) else: f2=kb*(T_cmb_old-T_m_old)/dr_ocean f1=k_fe*(T_c_old-T_cmb_old)/dr_ocean rm=radius[i_con,j] surfflux[j]=fs Alid=4*np.pi*(rm**2) Asurf=4*np.pi*(radius[i_s1,j]**2) if i_core==0: Vm=4/3*np.pi*(rm**3) Acmb=0 else: Vm=4/3*np.pi*((rm**3)-(radius[i_core-1,j]**3)) RADC=radius[i_core,j] Acmb=4*np.pi*(RADC**2) E_change=((-fs*Alid)+(h*rhom*Vm)+(f2*Acmb)) F1[j]=f1 F2[j]=f2 if TM[j-1]>=T_sil_s: cpm=2*850 else: cpm=850 dT=(E_change*dt)/(Vm*rhom*cpm) TM[j]=TM[j-1]+dT if TM[j]=0: Vc=((4/3)*np.pi)*((radius[i_core,j]**3)-(Rb**3)) delTC=(dt/(rho_l_fe_0*cpc*Vc))*((-f1)*Acmb) TC[j]=delTC+TC[j-1] for i in range(b_int,i_core,1): T[i,j]=TC[j] chi_sil[i,j]=0 #core - diffusive for i in range(0,b_int,1): if i==0: dTdr=0 d2Tdr2=(T[i+1,j-1]-T[i,j-1])/((radius[i+1,j-1])**2) else: dTdr=(T[i+1,j-1]-T[i-1,j-1])/(radius[i+1,j-1]-radius[i-1,j-1]) d2Tdr2=(1/dr_ocean)*((T[i+1,j-1]-T[i,j-1])/(radius[i+1,j-1]-radius[i,j-1])-(T[i,j-1]-T[i-1,j-1])/(radius[i,j-1]-radius[i-1,j-1])) T[i,j]=((k[i,j]/(rho_l_fe_0*cpc))*dt*((((1/radius[i,j-1])*dTdr)+d2Tdr2)))+T[i,j-1] chi_sil[i,j]=0 else: #core - diffusive for i in range(0,i_core,1): if i==0: dTdr=0 d2Tdr2=(T[i+1,j-1]-T[i,j-1])/((radius[i+1,j-1])**2) else: dTdr=(T[i+1,j-1]-T[i-1,j-1])/(radius[i+1,j-1]-radius[i-1,j-1]) d2Tdr2=(1/dr_ocean)*((T[i+1,j-1]-T[i,j-1])/(radius[i+1,j-1]-radius[i,j-1])-(T[i,j-1]-T[i-1,j-1])/(radius[i,j-1]-radius[i-1,j-1])) T[i,j]=((k[i,j]/(rho_l_fe_0*cpc))*dt*((((1/radius[i,j-1])*dTdr)+d2Tdr2)))+T[i,j-1] chi_sil[i,j]=0 radius_convecting[j]=radius[i_con,j] TC[j]=T[i_core-1,j] #for rayleigh number check i_old_2=i_con_old #looking for convection con_check=0 i_con_old=i_con i_con=-1 i=i_s1+1 convecting[i_s1,j]=-5 eta[i_s1,j]=1E30 while i_con<=0: i=i-1 critra[i,j]=20.9*((gamma*(T[i,j]-1800))**(4)) rayleigh[i,j]=(g[i,j]*a_sil_exp*rho[i,j]*rho[i,j]*(TM[j]-T[i,j])*(radius[i,j]**3)*cpm)/(eta[i,j]*k[i,j]) if rayleigh[i,j]>critra[i,j]: i_con=i if i_con>=i_s1: i_con=i_s1-1 break elif i<=i_con_old: i_con=i_con_old break dRa_check=rayleigh[i_con,j]-rayleigh[i_old_2,j-1] if T[i_con,j]>T[i_con,j-1]: dRa_check=0 #mass of core added and new core radius if i_con-i_con_old>0: i=i_con_old-1 m_l_fe[j]=0 mcoretot_old=0 mcoretot=0 mscore=0 dTs=0 mcontot=0 T_m_mixed=0 while i<=i_con: i=i+1 if i==0: vnode=(4/3)*np.pi*(radius[i,j]**3) else: vnode=(4/3)*np.pi*((radius[i,j]**3)-(radius[i-1,j]**3)) mnode=rho[i,j]*vnode mfenode=chi_fe[i,j]*mnode msilnode=mnode-mfenode dTsnode=mfenode*TM[j] msnode=chi_fe[i,j]*frac_S[i,j]*mnode mcoretot=mcoretot+mfenode mscore=mscore+msnode dTs=dTs+dTsnode mcontot=mcontot+msilnode dTm=msilnode*T[i,j] T_m_mixed=T_m_mixed+dTm if i>=i_con: break Tmixed=((mcontotold*TM[j])+T_m_mixed)/(mcontot+mcontotold) TM[j]=Tmixed m_l_fe[j]=mcoretot+m_l_fe[j-1] frac_S_core[j]=((frac_S_core[j-1]*m_l_fe[j-1])+mscore)/m_l_fe[j] frac_S_new=mscore/mcoretot Xs_at=frac_S_new/(frac_S_new+((m_mol_S/m_mol_Fe)*(1-frac_S_new))) Xs_at=Xs_at/100 mcontotold=mcontot+mcontotold else: m_l_fe[j]=m_l_fe[j-1] frac_S_core[j]=frac_S_core[j-1] mcontotold=mcontotold if m_l_fe[j]==0: rcore[j]=0 else: rcore[j]=((3*m_l_fe[j])/(4*np.pi*rho_l_fe_0))**(1/3) frac_S[:i_core,j]=frac_S_core[j] i_core_old=i_core Tpeak[j]=np.amax(T[:,j]) #new bit added - 29/05 -TCMB pre added material if i_core==0: TC[j]=TM[j] if TM[j]>=TC[j]: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) elif b_int>=0: const_TC=(k_fe**(3/4))*(((rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/4)) const_TM=(kb**(3/4))*(((rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/4)) TCMB[j]=((const_TC*TC[j])+(const_TM*TM[j]))/(const_TC+const_TM) else: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) eta_cmb_old=10**(a0+(a1*TCMB[j])+a2*(np.tanh((TCMB[j]+a3)/a4))) #what does core look like after addition of new material? if i_core>=1: if i_con>i_con_old: #only matters if adding new material to the core dT_core[:i_core,j]=TM[j]-T[:i_core,j] if np.any(dT_core[:,j]<0): b_stable=np.transpose(np.asarray(np.where(dT_core[:,j]<0))) b_int=np.asscalar(b_stable[0]) Rb=radius[b_int,j] db=radius[i_core,j]-Rb #starting volume average to get temp of convecting layer TSUM=0 i=b_int-1 mcoremix=0 while i=i_core-1: break Tc_mixed=(TSUM+(mcoretot*TCMB[j]))/(mcoremix+mcoretot) for i in range(b_int,i_core,1): T[i,j]=Tc_mixed TC[j]=Tc_mixed else: b_int=-1 TC[j]=TM[j] #general instability in the core if i_core>=1: if i_con==i_con_old: unstable[:i_core-1,j]=T[i_core-1,j]-T[:i_core-1,j] if np.any(unstable[:,j]<0): b_stable=np.transpose(np.asarray(np.where(unstable[:,j]<0))) b_int=np.asscalar(b_stable[0]) Rb=radius[b_int,j] db=radius[i_core,j]-Rb #starting volume average to get temp of convecting layer TSUM=0 i=b_int-1 mcoremix=0 while i=i_core-1: break Tc_mixed=(TSUM)/(mcoremix) for i in range(b_int,i_core,1): T[i,j]=Tc_mixed TC[j]=Tc_mixed junstable=j else: b_int=-1 TC[j]=TC[j] dcon[j]=db f_ad=(k_fe*a_exp*T[i_core-1,j]*g[i_core-1,j])/cpc if i_con==i_con_old: E1_ave=(f1*dt)+E1_ave else: E1_ave=(f1*dt)+E1_ave x=x+1 j_ave=int(round((j+j_congrow)/2)) time_ave[x]=time[j_ave] dt_ave=(time[j]-time[j_congrow])*(1E6*365*24*3600) if dt_avef_ad: fdrive=f1_ave[x]-f_ad u_mac=((t_spin*bigg*a_exp*RADC*fdrive)/(cpc))**(1/2) mag_rey_ave[x]=(u_mac*db)/mag_diff else: mag_rey_ave[x]=0 if f1>f_ad: fdrive=f1-f_ad u_mac=((t_spin*bigg*a_exp*RADC*fdrive)/(cpc))**(1/2) mag_rey[j]=(u_mac*db)/mag_diff else: mag_rey[j]=0 #new bit added - 29/05 if i_core==0: TC[j]=TM[j] if TM[j]>=TC[j]: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) elif b_int>=0: const_TC=(k_fe**(3/4))*(((rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/4)) const_TM=(kb**(3/4))*(((rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/4)) TCMB[j]=((const_TC*TC[j])+(const_TM*TM[j]))/(const_TC+const_TM) else: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) eta_cmb_old=10**(a0+(a1*TCMB[j])+a2*(np.tanh((TCMB[j]+a3)/a4))) if i_con==i_old_2: if dRa_check<0: j_corefull=j Tl_fe_core=Tl_fe0-(18*frac_S_core[j]) S_core_final=frac_S_core[j] E1_ave=(f1*dt)+E1_ave x=x+1 j_ave=int(round((j+j_congrow)/2)) time_ave[x]=time[j_ave]#time at which grow magma ocean dt_ave=(time[j]-time[j_congrow])*(1E6*365*24*3600) if dt_avef_ad: fdrive=f1_ave[x]-f_ad u_mac=((t_spin*bigg*a_exp*RADC*fdrive)/(cpc))**(1/2) mag_rey_ave[x]=(u_mac*db)/mag_diff else: mag_rey_ave[x]=0 x=x+1 time_ave[x]=time[j] f1_ave[x]=F1[j] mag_rey_ave[x]=mag_rey[j] dTcdt[x]=((TC[j]-TC[j-1])/dt)*(1E6*365*24*3600)#dTC/dt in K/Myr dcon_x[x]=db j_congrow=j xcount=0 exp_data[10,0]=radius[i_core-1,j]/radius[i_s1,j] exp_data[11,0]=radius[i_con-1,j]/radius[i_s1,j] exp_data[12,0]=time[j] exp_data[13,0]=S_core_final exp_data[14,0]=Tl_fe_core break if j>=j_max-1: break #now fully formed core - resuming previous code. #break if becomes unstable i.e. onset of partial convection in code #or if any node in core reaches liquidus temperature i_con_final=i_con junstable=-10 while j<=j_max: if b_int>0: j_unstable=j break j=j+1 modj=j%1000 t=t0+(dt*j) time[j]=t/(1E6*365*24*60*60) if i_core<1: TC[j]=TM[j-1] #new stuff T_c_old=TC[j-1] T_cmb_old=TCMB[j-1] T_m_old=TM[j-1] chi_fe[:i_core,j-1]=1 chi_fe_tot[:i_core,j-1]=1 k[:i_core,j]=k_fe phi[:i_core,j]=0 rho_sil[:i_core,j]=0 eta[:i_core,j]=eta_fe #changing magma ocean properties chi_fe[i_core:i_con,j-1]=0 chi_fe_tot[i_core:i_con,j]=0 rho[i_core:i_con,j]=rhom rho_sil[i_core:i_con,j]=rhom k[i_core:i_con,j]=kb phi[i_core:i_con,j]=0 rho[:i_core,j]=rho[:i_core,j-1] #gravity profile for i in range(0,i_s1+1,1): r[i]=radius[i,j-1] rho_old2[i]=rho[i,j-1] if i==0: g[i,j]=(bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max) else: g[i,j]=(((r[i-1]/r[i])**2)*g[i-1,j])+((bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max)) for i in range((-i_max+i_s1),-(i_max+1),-1): if i==(-i_max+i_s1): P[i,j]=0 #surface pressure =0 else: dP[i]=(rho[i,j-1]*g[i,j]*dr_max) P[i,j]=P[i+1,j]+dP[i] #compaction - only occurring in lid portion for i in range(i_con,i_s1+1,1): if phi[i,j-1]>0: phi_old=phi[i,j-1] #old value of phi g_phi=((1-phi_0)/(1-phi_old))**(2/3) f_phi=0.5*((3/(np.pi*((2*(g_phi**(1/2))*(3-g_phi))-3)))**(1/3)) P_grain=P[i,j]*((np.pi)/(2*(3**0.5)*((4*(3**0.5)*((1-phi_old)**(2/3))*f_phi*f_phi)-1))) dln_phi=A*(P_grain**(2/3))*(b**(-3))*np.exp(-Ea_comp/(Rgas*T[i,j-1]))*dt phi_old=phi[i,j-1] ln_phi_old=np.log(1-phi_old) ln_phi=ln_phi_old+dln_phi phi[i,j]=1-np.exp(ln_phi) if phi[i,j]<=0: phi[i,j]=0 else: phi[i,j]=phi[i,j-1] rho_sil[i,j]=rhom*(1-phi[i,j]) k[i,j]=kb*((np.exp((-4*phi[i,j])/phi_1)+np.exp(-4.4-((4*phi[i,j])/phi_2)))**(1/4)) #Krause et al, 2011, LPSC expression k[i,j]=kb*np.exp(100*phi[i,j]*(-0.1246)) rho[i,j]=((chi_fe_tot[i,j-1]/rho_fe)+((1-chi_fe_tot[i,j-1])/rho_sil[i,j]))**(-1) for i in range(0,i_con,1): radius[i,j]=dr_con*float(i+1) for i in range(i_con,i_s1+1,1): if i==0: radius[i,j]=((rho[i,j-1]/rho[i,j])**(1/3))*radius[i,j-1] dr[i,j]=(radius[i,j]) else: Vold=(4/3)*np.pi*((radius[i,j-1]**3)-(radius[i-1,j-1]**3)) Vbelow=(4/3)*np.pi*(radius[i-1,j]**3) radius[i,j]=((3/(4*np.pi))*(((rho[i,j-1]/rho[i,j])*Vold)+Vbelow))**(1/3) dr[i,j]=(radius[i,j]-radius[i-1,j]) Rp_c[j]=radius[i_s1,j] T[i_s1,j]=T0 #thermal evolution #magma ocean evolution h=h0*c0*((mcontotold+mcoretot)/(mcontotold))*np.exp((-1)*np.log(2)*t/thalf) racrit=1000 d0=(racrit**(1/3))*(((gamma*(TM[j-1]-T0))/8)**(4/3))*(((kb*eta_m_old)/(rhom*rhom*2*cpm*a_sil_exp*g[i_con,j]*(TM[j-1]-T0)))**(1/3)) #how does d0 affect lid thickness... is it bigger or smaller than node? d_i=int(d0//dr_ocean) i_con=i_con_final-d_i if i_con<=i_core: j_cond=j exp_data[15,0]=round(time[j],2) break fs=kb*((TM[j-1]-T0)/d0) #very different expression to corestratv3_3. if T_m_old>=T_c_old: #conductive fluxes, passing heat into the core f2=kb*(T_cmb_old-T_m_old)/dr_ocean f1=k_fe*(T_c_old-T_cmb_old)/dr_ocean elif b_int>=0: f1=k_fe*((T_c_old-T_cmb_old)**(4/3))*(((g[i_core,j]*rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/3)) f2=kb*((T_cmb_old-T_m_old)**(4/3))*(((g[i_core,j]*rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/3)) else: f2=kb*(T_cmb_old-T_m_old)/dr_ocean f1=k_fe*(T_c_old-T_cmb_old)/dr_ocean rm=radius[i_con,j] surfflux[j]=fs Alid=4*np.pi*(rm**2) Asurf=4*np.pi*(radius[i_s1,j]**2) if i_core==0: Vm=4/3*np.pi*(rm**3) Acmb=1 f2=0 else: Vm=4/3*np.pi*((rm**3)-(radius[i_core-1,j]**3)) RADC=radius[i_core,j] Acmb=4*np.pi*(RADC**2) E_change=((-fs*Alid)+(h*rhom*Vm)+(f2*Acmb)) F1[j]=f1 F2[j]=f2 if TM[j-1]>=T_sil_s: cpm=2*850 else: cpm=850 dT=(E_change*dt)/(Vm*rhom*cpm) TM[j]=TM[j-1]+dT if TM[j]=0.9999: chi[i,j]=1 if T[i,j]>Tl_fe: chi[i,j]=1 chi_fe[i,j]=x0_fe*chi[i,j] if T[i,j]>=Tl_fe0: T[i,j]=Tl_fe0 #viscosity stuff #silicate melt fraction - assuming linear melting - will need to modify the specific heat capacity - should probably add into cfl check if T[i,j]ceu: frac_S[i,j]=ceu #core - diffusive for i in range(0,i_core,1): if i==0: dTdr=0 d2Tdr2=(T[i+1,j-1]-T[i,j-1])/((radius[i+1,j-1])**2) else: dTdr=(T[i+1,j-1]-T[i-1,j-1])/(radius[i+1,j-1]-radius[i-1,j-1]) d2Tdr2=(1/dr_ocean)*((T[i+1,j-1]-T[i,j-1])/(radius[i+1,j-1]-radius[i,j-1])-(T[i,j-1]-T[i-1,j-1])/(radius[i,j-1]-radius[i-1,j-1])) T[i,j]=((k[i,j]/(rho_l_fe_0*cpc))*dt*((((1/radius[i,j-1])*dTdr)+d2Tdr2)))+T[i,j-1] chi_sil[i,j]=0 radius_convecting[j]=radius[i_con,j] TC[j]=T[i_core-1,j] #for rayleigh number check i_old_2=i_con_old #looking for convection con_check=0 i_con_old=i_con i_con=i_con i_core_old=i_core i_core=i_core Tpeak[j]=np.amax(T[:,j]) dT_core[:i_core-1,j]=T[i_core-1,j]-T[:i_core-1,j] if np.any(dT_core[:,j]<0): b_stable=np.transpose(np.asarray(np.where(dT_core[:,j]<0))) b_int=np.asscalar(b_stable[0]) Rb=radius[b_int,j] db=radius[i_core,j]-Rb Vc=((4/3)*np.pi)*((radius[i_core,j]**3)-(Rb**3)) #starting volume average to get temp of convecting layer TSUM=0 i=b_int-1 mcoremix=0 while i=i_core-1: break Tc_mixed=TSUM/mcoremix for i in range(b_int,i_core,1): T[i,j]=Tc_mixed TC[j]=Tc_mixed junstable=j eta_m_old=eta_m dcon[j]=db if TM[j]>=TC[j]: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) elif b_int>=0: const_TC=(k_fe**(3/4))*(((rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/4)) const_TM=(kb**(3/4))*(((rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/4)) TCMB[j]=((const_TC*TC[j])+(const_TM*TM[j]))/(const_TC+const_TM) else: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) eta_cmb_old=10**(a0+(a1*TCMB[j])+a2*(np.tanh((TCMB[j]+a3)/a4))) f_ad=(k_fe*a_exp*T[i_core-1,j]*g[i_core-1,j])/cpc xcount=xcount+1 if xcount==xcount_max: E1_ave=(f1*dt)+E1_ave x=x+1 j_ave=int(round((j+j_congrow)/2)) time_ave[x]=time[j_ave] dt_ave=(time[j]-time[j_congrow])*(1E6*365*24*3600) if dt_avef_ad: fdrive=f1_ave[x]-f_ad u_mac=((t_spin*bigg*a_exp*RADC*fdrive)/(cpc))**(1/2) mag_rey_ave[x]=(u_mac*db)/mag_diff else: mag_rey_ave[x]=0 else: E1_ave=(f1*dt)+E1_ave #magnetic Reynolds number if f1>f_ad: fdrive=f1-f_ad u_mac=((t_spin*bigg*a_exp*RADC*fdrive)/(cpc))**(1/2) mag_rey[j]=(u_mac*db)/mag_diff #x=x+1 #mag_rey_ave[x]=mag_rey[j] #time_ave[x]=time[j] else: mag_rey[j]=0 #if j>284: # print(fs,d0,TM[j-1]) if junstable>0: j_unstable=j #print(fs,d0,TM[j-1]) break if np.any(T[:i_core,j]<=Tl_fe_core): jfreeze=j exp_data[16,0]=round(time[j],2) break if j>=j_max-1: break if b_int>0: #partial core convection - changing heat fluxes in and out of cmb to convective fluxes... #const_TC=(k_fe**(3/4))*(((rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/4)) #const_TM=(kb**(3/4))*(((rhom*a_sil_exp)/(diffm*eta_m_old))**(1/4)) #TCMB[j]=((const_TC*TC[j])+(const_TM*TM[j]))/(const_TC+const_TM) while j<=j_max: j=j+1 modj=j%1000 t=t0+(dt*j) time[j]=t/(1E6*365*24*60*60) if i_core<1: TC[j]=TM[j-1] #new stuff T_c_old=TC[j-1] T_cmb_old=TCMB[j-1] T_m_old=TM[j-1] chi_fe[:i_core,j-1]=1 chi_fe_tot[:i_core,j-1]=1 k[:i_core,j]=k_fe phi[:i_core,j]=0 rho_sil[:i_core,j]=0 eta[:i_core,j]=eta_fe #changing magma ocean properties chi_fe[i_core:i_con,j-1]=0 chi_fe_tot[i_core:i_con,j]=0 rho[i_core:i_con,j]=rhom rho_sil[i_core:i_con,j]=rhom k[i_core:i_con,j]=kb phi[i_core:i_con,j]=0 rho[:i_core,j]=rho[:i_core,j-1] for i in range(0,i_s1+1,1): r[i]=radius[i,j-1] rho_old2[i]=rho[i,j-1] if i==0: g[i,j]=(bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max) else: g[i,j]=(((r[i-1]/r[i])**2)*g[i-1,j])+((bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max)) for i in range((-i_max+i_s1),-(i_max+1),-1): if i==(-i_max+i_s1): P[i,j]=0 #surface pressure =0 else: dP[i]=(rho[i,j-1]*g[i,j]*dr_max) P[i,j]=P[i+1,j]+dP[i] #compaction - only occurring in lid portion for i in range(i_con,i_s1+1,1): if phi[i,j-1]>0: phi_old=phi[i,j-1] #old value of phi g_phi=((1-phi_0)/(1-phi_old))**(2/3) f_phi=0.5*((3/(np.pi*((2*(g_phi**(1/2))*(3-g_phi))-3)))**(1/3)) P_grain=P[i,j]*((np.pi)/(2*(3**0.5)*((4*(3**0.5)*((1-phi_old)**(2/3))*f_phi*f_phi)-1))) dln_phi=A*(P_grain**(2/3))*(b**(-3))*np.exp(-Ea_comp/(Rgas*T[i,j-1]))*dt phi_old=phi[i,j-1] ln_phi_old=np.log(1-phi_old) ln_phi=ln_phi_old+dln_phi phi[i,j]=1-np.exp(ln_phi) if phi[i,j]<=0: phi[i,j]=0 else: phi[i,j]=phi[i,j-1] rho_sil[i,j]=rhom*(1-phi[i,j]) k[i,j]=kb*((np.exp((-4*phi[i,j])/phi_1)+np.exp(-4.4-((4*phi[i,j])/phi_2)))**(1/4)) #Krause et al, 2011, LPSC expression k[i,j]=kb*np.exp(100*phi[i,j]*(-0.1246)) rho[i,j]=((chi_fe_tot[i,j-1]/rho_fe)+((1-chi_fe_tot[i,j-1])/rho_sil[i,j]))**(-1) for i in range(0,i_con,1): radius[i,j]=dr_con*float(i+1) for i in range(i_con,i_s1+1,1): #over entire body including surface node - need to get planetary radius if i==0: radius[i,j]=((rho[i,j-1]/rho[i,j])**(1/3))*radius[i,j-1] dr[i,j]=(radius[i,j]) else: Vold=(4/3)*np.pi*((radius[i,j-1]**3)-(radius[i-1,j-1]**3)) Vbelow=(4/3)*np.pi*(radius[i-1,j]**3) radius[i,j]=((3/(4*np.pi))*(((rho[i,j-1]/rho[i,j])*Vold)+Vbelow))**(1/3) dr[i,j]=(radius[i,j]-radius[i-1,j]) Rp_c[j]=radius[i_s1,j] T[i_s1,j]=T0 #thermal evolution #magma ocean evolution h=h0*c0*((mcontotold+mcoretot)/(mcontotold))*np.exp((-1)*np.log(2)*t/thalf) racrit=1000 d0=(racrit**(1/3))*(((gamma*(TM[j-1]-T0))/8)**(4/3))*(((kb*eta_m_old)/(rhom*rhom*2*cpm*a_sil_exp*g[i_con,j]*(TM[j-1]-T0)))**(1/3)) d_i=int(d0//dr_ocean) i_con=i_con_final-d_i fs=kb*((TM[j-1]-T0)/d0) #very different expression to corestratv3_3. if T_m_old>=T_c_old: #conductive fluxes, passing heat into the core f2=kb*(T_cmb_old-T_m_old)/dr_ocean f1=k_fe*(T_c_old-T_cmb_old)/dr_ocean elif b_int>=0: f1=k_fe*((T_c_old-T_cmb_old)**(4/3))*(((g[i_core,j]*rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/3)) f2=kb*((T_cmb_old-T_m_old)**(4/3))*(((g[i_core,j]*rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/3)) else: f2=kb*(T_cmb_old-T_m_old)/dr_ocean f1=k_fe*(T_c_old-T_cmb_old)/dr_ocean rm=radius[i_con,j] surfflux[j]=fs Alid=4*np.pi*(rm**2) Asurf=4*np.pi*(radius[i_s1,j]**2) if i_core==0: Vm=4/3*np.pi*(rm**3) Acmb=1 f2=0 else: Vm=4/3*np.pi*((rm**3)-(radius[i_core-1,j]**3)) RADC=radius[i_core,j] Acmb=4*np.pi*(RADC**2) E_change=((-fs*Alid)+(h*rhom*Vm)+(f2*Acmb)) F1[j]=f1 F2[j]=f2 if TM[j-1]>=T_sil_s: cpm=2*850 else: cpm=850 dT=(E_change*dt)/(Vm*rhom*cpm) TM[j]=TM[j-1]+dT if TM[j]=0.9999: chi[i,j]=1 if T[i,j]>Tl_fe: chi[i,j]=1 chi_fe[i,j]=x0_fe*chi[i,j] if T[i,j]>=Tl_fe0: T[i,j]=Tl_fe0 #viscosity stuff #silicate melt fraction - assuming linear melting - will need to modify the specific heat capacity - should probably add into cfl check if T[i,j]ceu: frac_S[i,j]=ceu #core - diffusive for i in range(0,b_int,1): if i==0: dTdr=0 d2Tdr2=(T[i+1,j-1]-T[i,j-1])/((radius[i+1,j-1])**2) else: dTdr=(T[i+1,j-1]-T[i-1,j-1])/(radius[i+1,j-1]-radius[i-1,j-1]) d2Tdr2=(1/dr[i,j-1])*((T[i+1,j-1]-T[i,j-1])/(radius[i+1,j-1]-radius[i,j-1])-(T[i,j-1]-T[i-1,j-1])/(radius[i,j-1]-radius[i-1,j-1])) T[i,j]=((k[i,j]/(rho_l_fe_0*cpc))*dt*((((1/radius[i,j-1])*dTdr)+d2Tdr2)))+T[i,j-1] chi_sil[i,j]=0 radius_convecting[j]=radius[i_con,j] Vc=((4/3)*np.pi)*((radius[i_core,j]**3)-(Rb**3)) delTC=(dt/(rho_l_fe_0*cpc*Vc))*((-f1)*Acmb) TC[j]=delTC+TC[j-1] for i in range(b_int,i_core,1): T[i,j]=TC[j] chi_sil[i,j]=0 i_old_2=i_con_old con_check=0 i_con_old=i_con i_con=i_con i_core_old=i_core i_core=i_core Rb_old=Rb Tpeak[j]=np.amax(T[:,j]) dT_core[:i_core-1,j]=T[i_core-1,j]-T[:i_core-1,j] if np.any(dT_core[:,j]<0): b_stable=np.transpose(np.asarray(np.where(dT_core[:,j]<0))) b_int=np.asscalar(b_stable[0]) Rb=radius[b_int,j] db=radius[i_core,j]-Rb Vc=((4/3)*np.pi)*((radius[i_core,j]**3)-(Rb**3)) #starting volume average to get temp of convecting layer TSUM=0 i=b_int-1 mcoremix=0 while i=i_core-1: break Tc_mixed=TSUM/mcoremix for i in range(b_int,i_core,1): T[i,j]=Tc_mixed TC[j]=Tc_mixed junstable=j eta_m_old=eta_m dcon[j]=db #new bit added - 29/05 if TM[j]>=TC[j]: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) elif b_int>=0: const_TC=(k_fe**(3/4))*(((rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/4)) const_TM=(kb**(3/4))*(((rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/4)) TCMB[j]=((const_TC*TC[j])+(const_TM*TM[j]))/(const_TC+const_TM) else: TCMB[j]=(TM[j]+((k_fe/kb)*(dr_ocean/dr_ocean)*TC[j]))/(1+((k_fe/kb)*(dr_ocean/dr_ocean))) eta_cmb_old=10**(a0+(a1*TCMB[j])+a2*(np.tanh((TCMB[j]+a3)/a4))) f_ad=(k_fe*a_exp*T[i_core-1,j]*g[i_core-1,j])/cpc if Rb==Rb_old: E1_ave=(f1*dt)+E1_ave else: E1_ave=(f1*dt)+E1_ave x=x+1 j_ave=int(round((j+j_congrow)/2)) time_ave[x]=time[j_ave] dt_ave=(time[j]-time[j_congrow])*(1E6*365*24*3600) if dt_avef_ad: fdrive=f1_ave[x]-f_ad u_mac=((t_spin*bigg*a_exp*RADC*fdrive)/(cpc))**(1/2) mag_rey_ave[x]=(u_mac*db)/mag_diff else: mag_rey_ave[x]=0 dm=rm-RADC rayleigh_mantle=(g[i_con,j]*a_sil_exp*(TM[j]-T0)*rhom*(dm**3))/(diffm*eta_m) #print(rayleigh_mantle) #if j<300: # print(fs,d0) #magnetic Reynolds number fdrive=f1-f_ad if f1>f_ad: fdrive=f1-f_ad u_mac=((t_spin*bigg*a_exp*RADC*fdrive)/(cpc))**(1/2) mag_rey[j]=(u_mac*db)/mag_diff # x=x+1 # mag_rey_ave[x]=mag_rey[j] # time_ave[x]=time[j] else: mag_rey[j]=0 if b_int<=0: j_fullcon=j break if rayleigh_mantle<1000: j_cond=j exp_data[15,0]=time[j] break if np.any(T[:i_core,j]0: break j=j+1 modj=j%1000 t=t0+(dt*j) time[j]=t/(1E6*365*24*60*60) if i_core<1: TC[j]=TM[j-1] #new stuff T_c_old=TC[j-1] T_cmb_old=TCMB[j-1] T_m_old=TM[j-1] chi_fe[:i_core,j-1]=1 chi_fe_tot[:i_core,j-1]=1 k[:i_core,j]=k_fe phi[:i_core,j]=0 rho_sil[:i_core,j]=0 eta[:i_core,j]=eta_fe #changing magma ocean properties chi_fe[i_core:i_con,j-1]=0 chi_fe_tot[i_core:i_con,j]=0 rho[i_core:i_con,j]=rhom rho_sil[i_core:i_con,j]=rhom k[i_core:i_con,j]=kb phi[i_core:i_con,j]=0 rho[:i_core,j]=rho[:i_core,j-1] for i in range(0,i_s1+1,1): r[i]=radius[i,j-1] rho_old2[i]=rho[i,j-1] if i==0: g[i,j]=(bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max) else: g[i,j]=(((r[i-1]/r[i])**2)*g[i-1,j])+((bigg/(r[i]*r[i]))*(r[i]*r[i]*rho_old2[i]*dr_max)) for i in range((-i_max+i_s1),-(i_max+1),-1): if i==(-i_max+i_s1): P[i,j]=0 #surface pressure =0 else: dP[i]=(rho[i,j-1]*g[i,j]*dr_max) P[i,j]=P[i+1,j]+dP[i] #compaction - only occurring in lid portion for i in range(i_con,i_s1+1,1): if phi[i,j-1]>0: phi_old=phi[i,j-1] #old value of phi g_phi=((1-phi_0)/(1-phi_old))**(2/3) f_phi=0.5*((3/(np.pi*((2*(g_phi**(1/2))*(3-g_phi))-3)))**(1/3)) P_grain=P[i,j]*((np.pi)/(2*(3**0.5)*((4*(3**0.5)*((1-phi_old)**(2/3))*f_phi*f_phi)-1))) dln_phi=A*(P_grain**(2/3))*(b**(-3))*np.exp(-Ea_comp/(Rgas*T[i,j-1]))*dt phi_old=phi[i,j-1] ln_phi_old=np.log(1-phi_old) ln_phi=ln_phi_old+dln_phi phi[i,j]=1-np.exp(ln_phi) if phi[i,j]<=0: phi[i,j]=0 else: phi[i,j]=phi[i,j-1] rho_sil[i,j]=rhom*(1-phi[i,j]) k[i,j]=kb*((np.exp((-4*phi[i,j])/phi_1)+np.exp(-4.4-((4*phi[i,j])/phi_2)))**(1/4)) #Krause et al, 2011, LPSC expression k[i,j]=kb*np.exp(100*phi[i,j]*(-0.1246)) rho[i,j]=((chi_fe_tot[i,j-1]/rho_fe)+((1-chi_fe_tot[i,j-1])/rho_sil[i,j]))**(-1) for i in range(0,i_con,1): radius[i,j]=dr_con*float(i+1) for i in range(i_con,i_s1+1,1): #over entire body including surface node - need to get planetary radius if i==0: radius[i,j]=((rho[i,j-1]/rho[i,j])**(1/3))*radius[i,j-1] dr[i,j]=(radius[i,j]) else: Vold=(4/3)*np.pi*((radius[i,j-1]**3)-(radius[i-1,j-1]**3)) Vbelow=(4/3)*np.pi*(radius[i-1,j]**3) radius[i,j]=((3/(4*np.pi))*(((rho[i,j-1]/rho[i,j])*Vold)+Vbelow))**(1/3) dr[i,j]=(radius[i,j]-radius[i-1,j]) Rp_c[j]=radius[i_s1,j] T[i_s1,j]=T0 #thermal evolution #magma ocean evolution h=h0*c0*((mcontotold+mcoretot)/(mcontotold))*np.exp((-1)*np.log(2)*t/thalf) racrit=1000 d0=(racrit**(1/3))*(((gamma*(TM[j-1]-T0))/8)**(4/3))*(((kb*eta_m_old)/(rhom*rhom*2*cpm*a_sil_exp*g[i_con,j]*(TM[j-1]-T0)))**(1/3)) d_i=int(d0//dr_ocean) i_con=i_con_final-d_i if i_con<=i_core: j_cond=j exp_data[15,0]=time[j] break fs=kb*((TM[j-1]-T0)/d0) #very different expression to corestratv3_3. rm=radius[i_con,j] surfflux[j]=fs Alid=4*np.pi*(rm**2) Asurf=4*np.pi*(radius[i_s1,j]**2) if i_core==0: Vm=4/3*np.pi*(rm**3) Acmb=1 f2=0 Vc=0 else: Vm=4/3*np.pi*((rm**3)-(radius[i_core-1,j]**3)) RADC=radius[i_core,j] Acmb=4*np.pi*(RADC**2) Vc=(4/3)*np.pi*(RADC**3) f1=k_fe*((T_c_old-T_cmb_old)**(4/3))*(((g[i_core,j-1]*rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/3)) f2=kb*((T_cmb_old-T_m_old)**(4/3))*(((g[i_core,j-1]*rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/3)) E_change=((-fs*Alid)+(h*rhom*Vm)+(f2*Acmb)) F1[j]=f1 F2[j]=f2 if TM[j-1]>=T_sil_s: cpm=2*850 else: cpm=850 dT=(E_change*dt)/(Vm*rhom*cpm) TM[j]=TM[j-1]+dT if TM[j]=0.9999: chi[i,j]=1 if T[i,j]>Tl_fe: chi[i,j]=1 chi_fe[i,j]=x0_fe*chi[i,j] if T[i,j]>=Tl_fe0: T[i,j]=Tl_fe0 #viscosity stuff #silicate melt fraction - assuming linear melting - will need to modify the specific heat capacity - should probably add into cfl check if T[i,j]ceu: frac_S[i,j]=ceu radius_convecting[j]=radius[i_con,j] #change in T of convecting core portion delTC=(dt/(rho_l_fe_0*cpc*Vc))*((-f1)*Acmb) TC[j]=delTC+TC[j-1] for i in range(0,i_core,1): T[i,j]=TC[j] chi_sil[i,j]=0 i_old_2=i_con_old #looking for convection con_check=0 i_con_old=i_con i_con=i_con i_core_old=i_core i_core=i_core Tpeak[j]=np.amax(T[:,j]) #new bit added - 29/05 eta_m_old=eta_m const_TC=(k_fe**(3/4))*(((rho_l_fe_0*a_exp)/(diffc*eta_fe_l))**(1/4)) const_TM=(kb**(3/4))*(((rhom*a_sil_exp)/(diffm*eta_cmb_old))**(1/4)) TCMB[j]=((const_TC*TC[j])+(const_TM*TM[j]))/(const_TC+const_TM) eta_cmb_old=10**(a0+(a1*TCMB[j])+a2*(np.tanh((TCMB[j]+a3)/a4))) dcon[j]=RADC dm=rm-RADC rayleigh_mantle=(g[i_con,j]*a_sil_exp*(TM[j]-T0)*rhom*(dm**3))/(diffm*eta_m) # print(rayleigh_mantle) #print(rayleigh_mantle) xcount=xcount+1 if xcount==xcount_max: E1_ave=(f1*dt)+E1_ave x=x+1 j_ave=int(round((j+j_congrow)/2)) time_ave[x]=time[j_ave] dt_ave=(time[j]-time[j_congrow])*(1E6*365*24*3600) if dt_avef_ad: fdrive=f1_ave[x]-f_ad u_mac=((t_spin*bigg*a_exp*RADC*fdrive)/(cpc))**(1/2) mag_rey_ave[x]=(u_mac*db)/mag_diff else: mag_rey_ave[x]=0 else: E1_ave=(f1*dt)+E1_ave #magnetic Reynolds number f_ad=(k_fe*a_exp*T[i_core-1,j]*g[i_core-1,j])/cpc fdrive=f1-f_ad if f1>f_ad: fdrive=f1-f_ad u_mac=((t_spin*bigg*a_exp*RADC*fdrive)/(cpc))**(1/2) mag_rey[j]=(u_mac*db)/mag_diff #mag_rey_ave[x]=mag_rey[j] #time_ave[x]=time[j] else: mag_rey[j]=0 if np.any(T[:i_core,j]0: cutoff_time=np.transpose(np.asarray(np.where(time_ave>4))) if np.any(time_ave>4): x0=int(cutoff_time[0]) dynamo1=np.transpose(np.asarray(np.where(mag_rey_ave[x0:x]>=10))) if np.all(mag_rey_ave[x0:x]<10): td1start=0 td1end=0 else: d1start=int(np.min(dynamo1)+x0) d1end=int(np.max(dynamo1)+x0) td1start=round(time_ave[d1start],2) td1end=round(time_ave[d1end],2) else: td1start=0 td1end=0 x0=0 exp_data[20,0]=td1start exp_data[21,0]=td1end exp_data[22,0]=np.amax(mag_rey_ave[x0:x]) xmax=int(np.argmax(mag_rey_ave[x0:x])) dTcdt_max=-1*dTcdt[xmax+x0] dcon_max=dcon_x[xmax+x0] print(dTcdt_max,dcon_max) exp_data[27,0]=dTcdt_max exp_data[28,0]=dcon_max print(exp_data[22,0]) if td1start<4: if exp_data[22,0]>=10: spikes[:x,0]=time_ave[:x] spikes[:x,1]=mag_rey_ave[:x] spikes[:x,2]=dcon_x[:x] spikes[:x,3]=dTcdt[:x] spikes_file=('12E-7_TM_'+str(R0_int)+'_'+str(round(R/1E3))+'_'+str(round(exp_data[0,0],3))+'_'+str(round(exp_data[1,0],3))) else: exp_data[20,0]=0 exp_data[21,0]=0 exp_data[22,0]=0 exp_data[30,0]=radius[i_s1,j] end=TIME.time() sim_time=round((end-start)) exp_data[23,0]=sim_time print('time taken', sim_time)