-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathMain_code_callable.py
1846 lines (1544 loc) · 105 KB
/
Main_code_callable.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# This script contains the forward model
#####################
import numpy as np
import pylab
from scipy.integrate import *
from scipy.interpolate import interp1d
from scipy import optimize
import pdb
import scipy.optimize
from radiative_functions import *
from other_functions import *
from stellar_funs import main_sun_fun
from carbon_cycle_model import *
from escape_functions import *
from all_classes import *
#from outgassing_module import *
from outgassing_module_fast import *
from Albedo_module import *
from thermodynamic_variables import *
import time
from numba import jit
#####################
def forward_model(Switch_Inputs,Planet_inputs,Init_conditions,Numerics,Stellar_inputs,MC_inputs,max_time_attempt):
plot_switch = "n" # change to "y" to plot individual model runs for diagnostic purposes
print_switch = Switch_Inputs.print_switch # This controls whether outputs print during calculations (slows things down, but useful for diagnostics)
speedup_flag = Switch_Inputs.speedup_flag # Redundant - does not do anything in current version
start_speed = Switch_Inputs.start_speed # Redundant - does not do anything in current version
fin_speed = Switch_Inputs.fin_speed # Redundant - does not do anything in current version
heating_switch = Switch_Inputs.heating_switch # Controls locus of internal heating, keep default values
C_cycle_switch = Switch_Inputs.C_cycle_switch # Turns carbon cycle on or off, keep default values
RE = Planet_inputs.RE #Planet radius relative Earth
ME = Planet_inputs.ME #Planet mass relative Earth
pm = Planet_inputs.pm #Average mantle density
rc = Planet_inputs.rc #Metallic core radius (m)
Total_Fe_mol_fraction = Planet_inputs.Total_Fe_mol_fraction # iron mol fraction in mantle
Planet_sep = Planet_inputs.Planet_sep #planet-star separation (AU)
albedoC = Planet_inputs.albedoC #cold state albedo
albedoH = Planet_inputs.albedoH #hot state albedo
#Stellar parameters
tsat_XUV = Stellar_inputs.tsat_XUV #XUV saturation time
Stellar_Mass = Stellar_inputs.Stellar_Mass #stellar mass (relative sun)
fsat = Stellar_inputs.fsat
beta0 = Stellar_inputs.beta0
epsilon = Stellar_inputs.epsilon
#generate random seed for this forward model call
np.random.seed(int(time.time()))
seed_save = np.random.randint(1,1e9)
## Initial volatlie and redox conditions:
Init_solid_H2O = Init_conditions.Init_solid_H2O
Init_fluid_H2O = Init_conditions.Init_fluid_H2O
Init_solid_O= Init_conditions.Init_solid_O
Init_fluid_O = Init_conditions.Init_fluid_O
Init_solid_FeO1_5 = Init_conditions.Init_solid_FeO1_5
Init_solid_FeO = Init_conditions.Init_solid_FeO
Init_fluid_CO2 = Init_conditions.Init_fluid_CO2
Init_solid_CO2= Init_conditions.Init_solid_CO2
#Oxidation parameters
wet_oxid_eff = MC_inputs.interiord
MFrac_hydrated = MC_inputs.interiorb
dry_oxid_frac = MC_inputs.interiorc
surface_magma_fr = MC_inputs.surface_magma_frac
#ocean chemistry and weathering parameters
ocean_Ca = MC_inputs.ocean_a
omega_ocean = MC_inputs.ocean_b
efold_weath = MC_inputs.ccycle_a
alpha_exp = MC_inputs.ccycle_b
supp_lim = MC_inputs.supp_lim
#Escape parameters
mult = MC_inputs.esc_c
mix_epsilon = MC_inputs.esc_d
Te_input_escape = MC_inputs.Tstrat
#Interior parameters
transition_time_stag = MC_inputs.interiorg
visc_offset = MC_inputs.interiora
heatscale = MC_inputs.interiore
#impact parameters
imp_coef = MC_inputs.esc_a
tdc = MC_inputs.esc_b
MEarth = 5.972e24 #Mass of Earth (kg)
kCO2 = 2e-3 #Crystal-melt partition coefficent for CO2
#kCO2 = 0.99 #sensitivity test reduced mantle (CO2 retained in interior)
G = 6.67e-11 #gravitational constant
cp = 1.2e3 # silicate heat capacity
rp = RE * 6.371e6 #Planet radius (m)
Mp = ME * MEarth #Planet mass (kg)
delHf = 4e5 #Latent heat of silicates
g = G*Mp/(rp**2) # gravity (m/s2)
Tsolidus = sol_liq(rp,g,pm,rp,0.0,0.0) #Solidus for magma ocean evolution
Tliquid = Tsolidus + 600 #Liquidus for magma ocean evolution
alpha = 2e-5 #Thermal expansion coefficient (per K)
k = 4.2 #Thermal conductivity, W/m/K
kappa = 1e-6 #Thermal diffusivity of silicates, m2/s
Racr = 1.1e3 #Critical Rayeligh number
kH2O = 0.01 #Crystal-melt partition coefficent for water
a1 = 104.42e-9 #Solidus coefficient
b1 = 1420+0.000-80.0 #Solidus coefficient
a2 = 26.53e-9 #Solidus coefficient
b2 = 1825+0.000 #Solidus coefficient
min_Te = 150.0 ## Minimum Te for purposes of OLR/ASR calculations and escape calculations
min_ASR = 5.67e-8 * (min_Te/(0.5**0.25))**4.0 ## Minimum Absorbed Shortwave Radiation (ASR)
TMoffset = 0.0 ## Redundant
min_Te = 207.14285714 # Threshold to prevent skin temperature from getting too low where OLR grid contains errors. Note this lower limit does not apply to stratosphere temperatures used for escape calculations.
# Define radiogenic inventory
Uinit = heatscale*22e-9 #Uranium abundance
K_over_U = MC_inputs.K_over_U #K/U ratio
core_flow_coefficient = 12e12 #Assumed core heatflow
K40_over_K = 1.165e-4
lam_Ar = 0.0581
lam_Ca = 0.4962
init40K = (pm * 4.*np.pi * (rp**3 - rc**3)/3.0) * Uinit * K_over_U *K40_over_K * np.exp (4.5*(lam_Ar+lam_Ca))
initAr40 = 0.0
init_U238 = (pm * 4.*np.pi * (rp**3 - rc**3)/3.0) * 0.9928 * Uinit * np.exp((np.log(2)/4.46e9)*(4.5e9))
init_U235 = (pm * 4.*np.pi * (rp**3 - rc**3)/3.0) * 0.0072 * Uinit * np.exp((np.log(2)/7.04e8)*(4.5e9))
init_Th = (pm * 4.*np.pi * (rp**3 - rc**3)/3.0) *(0.069e-6) * heatscale * np.exp((np.log(2)/1.4e10)*(4.5e9))
init_He_mantle = 0.0
init_He_atmo = 0.0
He_escape_loss = 0.0
## D/H initialization - D/H evolution has not been fully implemented, keep switched off.
Init_D_to_H = 1.4e-4
DH_switch = 0.0 # 0 = off, 1 = on
Max_mantle_H2O = 1.4e21 * MC_inputs.interiorf * (rp**3 - rc**3) / ((6.371e6)**3 - (3.4e6)**3) ## Max mantle water content (kg)
Start_time = Switch_Inputs.Start_time #Model start time (relative to stellar evolution track)
Max_time=np.max([Numerics.tfin0,Numerics.tfin1,Numerics.tfin2,Numerics.tfin3,Numerics.tfin4]) #Model end time
test_time = np.linspace(Start_time*365*24*60*60,Max_time*365*24*60*60,10000)
new_t = np.linspace(Start_time/1e9,Max_time/1e9,100000)
[Relative_total_Lum,Relative_XUV_lum,Absolute_total_Lum,Absolute_XUV_Lum] = main_sun_fun(new_t,Stellar_Mass,tsat_XUV,beta0,fsat) #Calculate stellar evolution
ASR_new = (Absolute_total_Lum/(16*3.14159*(Planet_sep*1.496e11)**2) ) #ASR flux through time (not accounting for bond albedo)
for ij in range(0,len(ASR_new)): # do not permit ASR outside of interpolation grid
if ASR_new[ij] < min_ASR:
ASR_new[ij] = min_ASR
Te_ar = (ASR_new/5.67e-8)**0.25
Tskin_ar = Te_ar*(0.5**0.25) ## Skin temperature through time
for ij in range(0,len(Tskin_ar)): #Don't permit skin temperature to exceed range min_Te - 350 due to errors in grid (does not apply to stratospheric temperature used to calculate escape fluxes)
if Tskin_ar[ij] > 350:
Tskin_ar[ij] = 350.0
if Tskin_ar[ij] < min_Te:
Tskin_ar[ij] = min_Te
Te_fun = interp1d(new_t*1e9*365*24*60*60,Tskin_ar) #Skin temperature function, used in OLR calculations
ASR_new_fun = interp1d(new_t*1e9*365*24*60*60, ASR_new) #ASR function, used to calculate shortwave radiation fluxes through time
AbsXUV = interp1d(new_t*1e9*365*24*60*60 , Absolute_XUV_Lum/(4*np.pi*(Planet_sep*1.496e11)**2)) #XUV function, used to calculate XUV-driven escape
#@jit(nopython=True) # function for finding surface temperature that balances ASR and interior heatflow
def funTs_general(Ts,Tp,ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,MMW,PO2_surf):
if max(Ts)<Tp:
Ra =alpha * g * (Tp -Ts) * ll**3 / (kappa * visc)
qm = (k/ll) * (Tp - Ts) * (Ra/Racr)**beta
else:
qm = - 2.0 * (Ts - Tp) / (rp-rc)
Ts_in= max(Ts)
[OLR_new,newpCO2,ocean_pH,ALK,Mass_oceans_crude,DIC_check] = correction(Ts_in,Te_input,H2O_Pressure_surface,CO2_Pressure_surface,rp,g,ocean_CO3,1e5,MMW,PO2_surf)
heat_atm = OLR_new - ASR_input
return (qm - heat_atm)**2
@jit(nopython=True) # as above but precompiled
def funTs_general2(Ts,Tp,ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,MMW,PO2_surf):
if Ts[0]<Tp:
Ra =alpha * g * (Tp -Ts[0]) * ll**3 / (kappa * visc)
qm = (k/ll) * (Tp - Ts[0]) * (Ra/Racr)**beta
else:
qm = - 2.0 * (Ts[0] - Tp) / (rp-rc)
Ts_in= Ts[0]
[OLR_new,newpCO2,ocean_pH,ALK,Mass_oceans_crude,DIC_check] = correction(Ts_in,Te_input,H2O_Pressure_surface,CO2_Pressure_surface,rp,g,ocean_CO3,1e5,MMW,PO2_surf)
heat_atm = OLR_new - ASR_input
return -(qm - heat_atm)**2
#Radiation balance accommodating stagnant lid and plate tectonics:
def funTs_general_stag(Ts,Tp,ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,MMW,PO2_surf,MV,PT):
if ((max(Ts)<Tp)and(PT < 1)):
Rai = alpha * g * (Tp -Ts) * ll**3 / (kappa * visc)
theta = 300000*(Tp - Ts)/(8.314*Tp**2)
qa = 0.5 * k/ll * (Tp - Ts) * theta**(-4.0/3.0) * Rai**(1./3.)
deltaTm = Tp -Ts
volc_term = MV*3000*(cp*deltaTm + delHf)
qm = qa + volc_term/(4*np.pi*rp**2)
elif ((max(Ts)<Tp)and(PT >= 1)):
Ra =alpha * g * (Tp -Ts) * ll**3 / (kappa * visc)
qm = (k/ll) * (Tp - Ts) * (Ra/Racr)**beta
else:
qm = - 2.0 * (Ts - Tp) / (rp-rc)
Ts_in= max(Ts)
[OLR_new,newpCO2,ocean_pH,ALK,Mass_oceans_crude,DIC_check] = correction(Ts_in,Te_input,H2O_Pressure_surface,CO2_Pressure_surface,rp,g,ocean_CO3,1e5,MMW,PO2_surf)
heat_atm = OLR_new - ASR_input
return (qm - heat_atm)**2
@jit(nopython=True) # as above but precompiled
def funTs_general_stag2(Ts,Tp,ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,MMW,PO2_surf,MV,PT):
if Ts[0]<Tp:#and(PT<1.0)):
if PT < 1.0:
Rai = alpha * g * (Tp -Ts[0]) * ll**3 / (kappa * visc)
theta = 300000*(Tp - Ts[0])/(8.314*Tp**2)
qa = 0.5 * k/ll * (Tp - Ts[0]) * theta**(-4.0/3.0) * Rai**(1./3.)
deltaTm = Tp -Ts[0]
volc_term = MV*3000*(cp*deltaTm + delHf)
qm = qa + volc_term/(4*np.pi*rp**2)
if PT >= 1.0:
Ra =alpha * g * (Tp -Ts[0]) * ll**3 / (kappa * visc)
qm = (k/ll) * (Tp - Ts[0]) * (Ra/Racr)**beta
else:
qm = - 2.0 * (Ts[0] - Tp) / (rp-rc)
Ts_in= Ts[0]
[OLR_new,newpCO2,ocean_pH,ALK,Mass_oceans_crude,DIC_check] = correction(Ts_in,Te_input,H2O_Pressure_surface,CO2_Pressure_surface,rp,g,ocean_CO3,1e5,MMW,PO2_surf)
heat_atm = OLR_new - ASR_input
return -(qm - heat_atm)**2
#@jit(nopython=True) # alternative function for finding surface temperature that balances ASR and interior heatflow (does exact same thing, hardly used)
def funTs_scalar(Ts,Tp,ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,MMW,PO2_surf):
if Ts<Tp:
Ra =alpha * g * (Tp -Ts) * ll**3 / (kappa * visc)
qm = (k/ll) * (Tp - Ts) * (Ra/Racr)**beta
else:
qm = - 2.0 * (Ts - Tp) / (rp-rc)
Ts_in= Ts
[OLR_new,newpCO2,ocean_pH,ALK,Mass_oceans_crude,DIC_check] = correction(Ts_in,Te_input,H2O_Pressure_surface,CO2_Pressure_surface,rp,g,ocean_CO3,1e5,MMW,PO2_surf)
heat_atm = OLR_new - ASR_input
return (qm - heat_atm)**2
def FeO_mass_frac(Total_Fe): # Convert total iron mole fraction to mass fraction
XAl2O3 = 0.022423
XCaO = 0.0335
XNa2O = 0.0024
XK2O = 0.0001077
XMgO = 0.478144
XSiO2 = 0.4034
m_sil = XMgO * (16.+25.) + XSiO2 * (28.+32.) + XAl2O3 * (27.*2.+3.*16.) + XCaO * (40.+16.) + (Total_Fe) * (56.0+16.0)
return (Total_Fe) * (56.0+16.0)/m_sil
Total_Fe_mass_fraction = FeO_mass_frac(Total_Fe_mol_fraction)
global solid_switch,phi_final, ynewIC,switch_counter,speedup,liquid_switch,solid_counter,last_magma
# Define switches used to keep track of volatile transfers between magma ocean and solid mantle convection phases (and vice versa)
solid_switch = 0
solid_switch = 1 #for starting as solid only
liquid_switch = 0
liquid_switch_worked = 0
speedup = 0.0
switch_counter = 0
solid_counter = 0
last_magma = 0.0
model_run_time = time.time()
def system_of_equations(t0,y):
PltTech = 1 #Plate tectonics marker
global last_magma
crt = t0/(365*24*60*60)
if transition_time_stag > 0: # if negative, always plate tectoncis
if (crt > transition_time_stag)and((crt - last_magma)>50e6):
PltTech = 0 # Denotes stagnant lid
## For experimenting with plate tectonics regime change
#if ((crt >3.9e9) or ((crt < 3.0e9)and(crt>2.0e9)) )and((crt - last_magma)>50e6):# switching
#if (crt > 13e6)and((crt - last_magma)>50e6):# and (t0/(365*24*60*60) < 20e6): #staggers
#if (crt > 3.9e9)and((crt - last_magma)>50e6):# and (crt < 20e6): # PT until 4.0 Ga
# PltTech = 0
## For experimenting with artifical heatflow injection
Qr_cofactor = 1.0
#if (crt >3.5e9)and(crt<4.0e9):
# Qr_cofactor = 3.0
#if (crt >3.9e9)and(crt<3.95e9): #more realistic LIP, injects 50 TW for 50 Myrs USING THIS ONE
# Qr_cofactor = 4.0
tic = time.time()
if (tic - model_run_time)>60*60*max_time_attempt: #limit maximum time spent before abandon attempt
print ("TIMED OUT")
return np.nan*[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
Va = 0.0
actual_visc = 0.0
global solid_switch,speedup,liquid_switch,liquid_switch_worked
global phi_final,ynewIC,switch_counter,solid_counter
depletion_fraction = y[52]
speedup=0.0
ocean_CO3 = float(omega_ocean * Sol_prod(y[8]) / ocean_Ca)
Mantle_mass = (4./3. * np.pi * pm * (y[2]**3 - rc**3)) #Mass solid mantle
Mantle_mass0 = (4./3. * np.pi * pm * (rp**3 - rc**3)) #Mass solid mantle
Vmantle0 = (4.0*np.pi/3.0) * (rp**3 - rc**3) #Volume total mantle
if print_switch == "y":
print (t0/(365*24*60*60))
#################################################################################
#### If in magma ocean phase
if (y[8] > Tsolidus):
last_magma = t0/(365*24*60*60)
beta = 1./3. #Convective heatflow exponent
if print_switch== "y":
print('still molten',t0/(365*24*60*60),np.transpose(y))
#For switching from solid to magma ocean
if liquid_switch == 1:
if y[2]+1 < rp:
liquid_switch_worked = 1.0
liquid_switch = 0.0
else:
T_partition = np.max([y[7],Tsolidus+0.01])
rad_check = optimize.minimize(find_r,x0=float(y[2]),args = (T_partition,alpha,g,cp,pm,rp,0.0,0.0,0.0))
y[2] = np.max([rad_check.x[0],0.0])
#Calculate surface melt fraction
if y[8] > Tliquid:
actual_phi_surf = 1.0
elif y[8] < Tsolidus:
actual_phi_surf = 0.0
else:
actual_phi_surf =( y[8] -Tsolidus)/(Tliquid - Tsolidus)
ll = np.max([rp - y[2],1.0]) ## length scale is depth of magma ocean pre-solidification (even if melt <0.4)
Qr = Qr_cofactor*qr(t0,Start_time,heatscale,K_over_U)
Qcore = Qr_cofactor*np.exp(-(t0/(1e9*365*24*60*60)-4.5)/7.0)*core_flow_coefficient
Mliq = Mliq_fun(y[1],rp,y[2],pm)
Mcrystal = (1-actual_phi_surf)*Mliq
phi_final = actual_phi_surf
[FH2O,H2O_Pressure_surface] = H2O_partition_function( y[1],Mliq,Mcrystal,rp,g,kH2O,y[24])
[FCO2,CO2_Pressure_surface] = CO2_partition_function( y[12],Mliq,Mcrystal,rp,g,kCO2,y[24]) #molten so can ignore aqueous CO2
AB = AB_fun(float(y[8]),H2O_Pressure_surface,float(y[1]+y[4]+y[12]),albedoC,albedoH)
Tsolidus_Pmod = sol_liq(rp,g,pm,rp,float(H2O_Pressure_surface+CO2_Pressure_surface+1e5),float(0*y[0]/Mantle_mass0))
Tsolidus_visc = sol_liq(rp,g,pm,rp,float(H2O_Pressure_surface+CO2_Pressure_surface+1e5),0.0)
visc = viscosity_fun(y[7],pm,visc_offset,y[8],float(Tsolidus_visc))
if print_switch=="y":
print ('visc',visc)
if np.isnan(visc):
return [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
ASR_input = float((1-AB)*ASR_new_fun(t0))
if ASR_input < min_ASR:
ASR_input = min_ASR
Te_ar = (ASR_input/5.67e-8)**0.25
Te_input = Te_ar*(0.5**0.25)
if Te_input > 350:
Te_input = 350.0
if Te_input < min_Te:
Te_input = min_Te
if (2>1):
initialize_fast = np.min([y[8],y[7]+TMoffset,y[7]+TMoffset-1])
if y[8] > y[7]+TMoffset:
initialize_fast = y[8]
initialize_fast = np.array(initialize_fast)
ace1 = nelder_mead(funTs_general2, x0=initialize_fast, bounds=np.array([[100.0], [4500.0]]).T, args = (float(y[7]),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22])), tol_f=0.0001,tol_x=0.0001, max_iter=1000)
SurfT = ace1.x[0]
new_abs = abs(ace1.fun)
if new_abs > 0.1:
lower = 180.0
upper = float(y[7]+TMoffset)+10
ace1b=scipy.optimize.minimize_scalar(funTs_scalar, args=(float(y[7]+TMoffset),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22])),bounds=[lower,upper],tol=1e-10,method='bounded',options={'maxiter':1000,'xatol':1e-10})
if abs(ace1b.fun) < new_abs:
SurfT = float(ace1b.x)
new_abs = abs(ace1b.fun)
y[8] = SurfT
if new_abs > 1.0:
ace2= optimize.minimize(funTs_general,x0=y[8],args = (float(y[7]),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22])),method='L-BFGS-B',bounds = ((150,1500),))
if abs(ace2.fun) < new_abs:
SurfT = ace2.x[0]
new_abs = abs(ace2.fun)
if new_abs > 1.0:
rand_start = 150 + 1850*np.random.uniform()
ace3= optimize.minimize(funTs_general,x0=rand_start,args = (float(y[7]),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22])),method='L-BFGS-B',bounds = ((150,2000),))
if abs(ace3.fun) < new_abs:
SurfT =ace3.x[0]
new_abs = abs(ace3.fun)
if new_abs > 1.0:
differ_ace = 10.0
counter = 0
initialize_fast = np.array(initialize_fast)
budget_differ = 10.0
aceCOMP = nelder_mead(funTs_general2, x0=initialize_fast, bounds=np.array([[150.0], [4500.0]]).T, args = (float(y[7]),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22])), tol_f=1e-10,tol_x=1e-10, max_iter=1000)
if abs(aceCOMP.fun)<new_abs:
SurfT = aceCOMP.x[0]
budget_differ = aceCOMP.fun
budget_counter = 0
while abs(budget_differ) > 0.1:
budget_counter = budget_counter + 1
ran_num = np.min([150 + 2000*np.random.uniform(),3999])
ran_numar=np.array(ran_num)
aceCOMP2 = nelder_mead(funTs_general2, x0=ran_numar, bounds=np.array([[150.0], [4500.0]]).T, args = (float(y[7]),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22])), tol_f=1e-10,tol_x=1e-10, max_iter=2000)
if abs(aceCOMP2.fun)<abs(budget_differ):
SurfT = aceCOMP2.x[0]
budget_differ = aceCOMP2.fun
if budget_counter == 10:
budget_differ = -2e-3
Ra =np.max([-1e15,alpha * g * (y[7]+TMoffset -SurfT) * ll**3 / (kappa * visc) ])
qm = (k/ll) * (y[7]+TMoffset - SurfT) * (Ra/Racr)**beta
thermal_cond = 2 # W/m2/K
qc = thermal_cond * (SurfT - y[7]- TMoffset) / (rp-rc)
delta_u = k * (y[7]+TMoffset - SurfT) / qm
## check to see if anywhere near solidus
if adiabat(rc,y[7],alpha,g,cp,rp) > sol_liq(rc,g,pm,rp,0.0,0.0): #ignore pressure overburden when magma ocean solidifying (volatiles mostly in mantle)
rs_term = 0.0
else: # if at solidus, start increasing it according to mantle temperature cooling
rs_term = rs_term_fun(float(y[2]),a1,b1,a2,b2,g,alpha,cp,pm,rp,y[7],0.0) #doesnt seem to affect things
if (2>1):
if y[7]+TMoffset>SurfT:
if heating_switch == 1:
numerator = (4./3.)* np.pi * pm * Qr * (rp**3.0-y[2]**3.0) + Qcore - 4.0*np.pi*((rp)**2)*qm
elif heating_switch == 0 :
numerator = (4./3.)* np.pi * pm * Qr * (rp**3.0-rc**3.0) + Qcore - 4*np.pi*((rp)**2)*qm
else:
print ('ERROR')
denominator = (4./3.) * np.pi * pm * cp *(rp**3 - y[2]**3) - 4*np.pi*(y[2]**2)* delHf*pm * rs_term
dy7_dt = numerator/denominator #this is Tp
else:
if heating_switch == 1:
numerator = (4./3.)* np.pi * pm * Qr * (rp**3.0-y[2]**3.0) + Qcore + 4.0*np.pi*((rp)**2)*qc
elif heating_switch == 0:
numerator = (4./3.)* np.pi * pm * Qr * (rp**3.0-rc**3.0) + Qcore + 4*np.pi*((rp)**2)*qc
else:
print ('ERROR')
denominator = (4./3.) * np.pi * pm * cp *(rp**3 - y[2]**3) - 4*np.pi*(y[2]**2)* delHf*pm * rs_term
dy7_dt = numerator/denominator
dy16_dt = 0.0
[OLR,newpCO2,ocean_pH,ALK,Mass_oceans_crude,DIC_check] = correction(float(SurfT),float(Te_input),H2O_Pressure_surface,CO2_Pressure_surface,rp,g,ocean_CO3,1e5,float(y[24]),float(y[22]))
CO2_Pressure_surface = newpCO2
ddelta_dt = 0.0
dVcrust_dt = 0.0
dQmantle_dt = 0.0 #Qr # per unit volume
dQcrust_dt = 0.0 #
dK40mantle_dt = - (lam_Ar+lam_Ca) * y[35] / (1e9*365*24*60*60)
dAr40mantle_dt = (lam_Ar/(lam_Ar+lam_Ca)) * (y[2]**3 - rc**3)/(rp**3 - rc**3) * (lam_Ar+lam_Ca) * y[35] / (1e9*365*24*60*60)
dK40lid_dt = 0.0
dAr40atmo_dt = (lam_Ar/(lam_Ar+lam_Ca)) * (rp**3 - y[2]**3)/(rp**3 - rc**3) * (lam_Ar+lam_Ca) * y[35] / (1e9*365*24*60*60)
dU238_mantle_dt = -(np.log(2)/4.46e9) * y[43] / (365*24*60*60)
dU235_mantle_dt = -(np.log(2)/7.04e8) * y[44] / (365*24*60*60)
dTh_mantle_dt = - (np.log(2)/1.4e10) * y[45] / (365*24*60*60)
dU238_lid_dt = 0.0
dU235_lid_dt = 0.0
dTh_lid_dt = 0.0
dHe_mantle_dt = (y[2]**3 - rc**3)/(rp**3 - rc**3) * ( (4./238) * 8.0 * (np.log(2)/4.46e9) * y[43] + (4./235) * 7.0 * (np.log(2)/7.04e8) * y[44]+ (4./232) * 6.0 * (np.log(2)/1.4e10) * y[45]) / (365*24*60*60)
dHe_atmo_dt = (rp**3 - y[2]**3)/(rp**3 - rc**3) * ( (4./238) * 8.0 * (np.log(2)/4.46e9) * y[43] + (4./235) * 7.0 * (np.log(2)/7.04e8) * y[44]+ (4./232) * 6.0 * (np.log(2)/1.4e10) * y[45]) / (365*24*60*60) #- He_escape_loss
ASR = ASR_input
y[8] = SurfT
heat_atm = OLR - ASR
if (2>1):
if SurfT > y[7]+TMoffset:
true_balance = - heat_atm - qc
else:
true_balance = - heat_atm + qm
if print_switch == "y":
print ("phi",actual_phi_surf,"visc",visc)
print ("time",t0/(365*24*60*60*1e6),"Ra",Ra)
print (OLR,ASR)
print ("Heat balance",true_balance)
print (" ")
solid_switch = 0.0
if liquid_switch ==1 : ## If transitioning from liquid to solid, adjust inventories
T_partition = np.max([y[7],Tsolidus+0.01])
[XFeO,XFe2O3,fO2,F_FeO1_5,F_FeO] = solve_fO2_F_redo(y[4],H2O_Pressure_surface,T_partition,Total_Fe_mol_fraction,Mliq,y[2],rp,y[24])
mu_O = 16.0
mu_FeO_1_5 = 56.0 + 1.5*16.0
if liquid_switch_worked==0: #first time or hasn't worked yet
y[3] = y[3] - Mliq * (y[3] / (4./3. * np.pi * pm * (y[2]**3 - rc**3)))
y[4] = y[4] + Mliq * (y[3] / (4./3. * np.pi * pm * (y[2]**3 - rc**3)))
y[5] = y[5] - Mliq * (y[5] / (4./3. * np.pi * pm * (y[2]**3 - rc**3)))
y[6] = y[6] - Mliq * (y[6] / (4./3. * np.pi * pm * (y[2]**3 - rc**3)))
y[0] = y[0] - Mliq * (y[0] / (4./3. * np.pi * pm * (y[2]**3 - rc**3)))
y[1] = y[1] + Mliq * (y[0] / (4./3. * np.pi * pm * (y[2]**3 - rc**3)))
y[42] = y[42] - DH_switch*Mliq * (y[42] / (4./3. * np.pi * pm * (y[2]**3 - rc**3)))
y[41] = y[41] + DH_switch*Mliq * (y[42] / (4./3. * np.pi * pm * (y[2]**3 - rc**3)))
y[13] = y[13] - Mliq * (y[13] / (4./3. * np.pi * pm * (y[2]**3 - rc**3)))
y[12] = y[12] + Mliq * (y[13] / (4./3. * np.pi * pm * (y[2]**3 - rc**3)))
y[26] = 0.0
y[27] = 0.0
if y[29] > 0: # put everything back in the mantle. Everythin gets mixed up again during initial plate tectonics?
y[28] = y[28] + y[29]
y[29] = 0.0
deltacrust = rp - (rp**3 - 3*y[27]/(4.0*np.pi))**(1./3.)
if rp - y[2] > deltacrust:
y[38] = y[38] + y[36] * ((rp-deltacrust)**3 - y[2]**3)/((rp-deltacrust)**3 - rc **3 )
y[36] = y[36] - y[36] * ((rp-deltacrust)**3 - y[2]**3)/((rp-deltacrust)**3 - rc **3 )
y[35] = y[35] + y[37]
y[37] = 0.0
#dy35_dt = dK40mantle_dt
#dy36_dt = dAr40mantle_dt
#dy37_dt = dK40lid_dt
#dy38_dt = dAr40atmo_dt
y[50] = y[50] + y[49] * ((rp-deltacrust)**3 - y[2]**3)/((rp-deltacrust)**3 - rc **3 )
y[49] = y[49] - y[49] * ((rp-deltacrust)**3 - y[2]**3)/((rp-deltacrust)**3 - rc **3 )
y[43] = y[43] + y[46]
y[46] = 0.0
y[44] = y[44] + y[47]
y[47] = 0.0
y[45] = y[45] + y[48]
y[48] = 0.0
else:
y[38] = y[38] ## no unmelted rock gets melted, no new Ar added atmosphere
y[35] = y[35]+ y[37] ## crust gets mixed back into mantle
y[37] = 0.0
y[39] = 0.0
y[50] = y[50]
y[43] = y[43] + y[46]
y[44] = y[44] + y[47]
y[45] = y[45] + y[48]
y[46] = 0.0
y[47] = 0.0
y[48] = 0.0
switch_name = "switch_garbage/switch_IC_%d" %seed_save
load_name2 = switch_name+".npy"
if switch_counter == 0 :
np.save(switch_name,y)
else:
y = np.load(load_name2)
switch_counter = switch_counter + 1
return [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
else: #switch has worked, back to normal, reset for next magma ocean transition
liquid_switch_worked = 0.0
liquid_switch = 0.0
switch_counter = 0.0
return [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
#################################################################################
#### If in solid mantle phase
else:
if print_switch =="y":
print('Mantle fully solid',t0/(365*24*60*60),y)
beta = 1./3.
ll = rp - rc #this is just mantle thickness
if solid_switch == 0:
Mliq = Mliq_fun(y[1],rp,y[2],pm)
Mcrystal = (1-phi_final)* Mliq
[FH2O,H2O_Pressure_surface] = H2O_partition_function( y[1],Mliq,Mcrystal,rp,g,kH2O,y[24])
[FCO2,CO2_Pressure_surface] = CO2_partition_function( y[12],Mliq,Mcrystal,rp,g,kCO2,y[24])
T_partition = np.max([y[7],Tsolidus+0.01])
[XFeO,XFe2O3,fO2,F_FeO1_5,F_FeO] = solve_fO2_F_redo(y[4],H2O_Pressure_surface,T_partition,Total_Fe_mol_fraction,Mliq,y[2],rp,y[24])
mu_O = 16.0
mu_FeO_1_5 = 56.0 + 1.5*16.0
if y[2] < rp:
y[3] = y[3] + (4./3.) *np.pi * pm *F_FeO1_5 * (rp**3 - y[2]**3) * 0.5*mu_O / (mu_FeO_1_5)
y[4] = y[4] - (4./3.) *np.pi * pm *F_FeO1_5 * (rp**3 - y[2]**3) * 0.5*mu_O / (mu_FeO_1_5)
y[5] = y[5] + (4./3.) *np.pi * pm * F_FeO1_5 *(rp**3 - y[2]**3) ## mass solid FeO1_5
y[6] = y[6] + (4./3.) *np.pi * pm * F_FeO * (rp**3 - y[2]**3) #mass solid F_FeO
Water_transfer = np.min([Max_mantle_H2O-y[0],FH2O * kH2O * Mcrystal + FH2O * (Mliq-Mcrystal)])
CO2_transfer = np.min([Max_mantle_H2O-y[13],kCO2 * FCO2 * Mcrystal + FCO2 * (Mliq-Mcrystal) ])
DHO_transfer = DH_switch*Water_transfer*y[41]/y[1]
y[42] = y[42] + DHO_transfer
y[41] = y[41] - DHO_transfer
y[0] = y[0] + Water_transfer
y[1] = y[1] - Water_transfer
y[13] = y[13] + CO2_transfer
y[12] = y[12] - CO2_transfer
y[26] = 0.0
y[27] = 0.0
y[39] = 0.0
y[37] = pm*y[27]*(y[35]/Mantle_mass)
y[35] = y[35] - pm*y[27]*(y[35]/Mantle_mass)
# dy35_dt = dK40mantle_dt
# dy36_dt = dAr40mantle_dt
# dy37_dt = dK40lid_dt
# dy38_dt = dAr40atmo_dt
y[46] = pm*y[27]*(y[43]/Mantle_mass)
y[47] = pm*y[27]*(y[44]/Mantle_mass)
y[48] = pm*y[27]*(y[45]/Mantle_mass)
y[43] = y[43] - pm*y[27]*(y[43]/Mantle_mass)
y[44] = y[44] - pm*y[27]*(y[44]/Mantle_mass)
y[45] = y[45] - pm*y[27]*(y[45]/Mantle_mass)
y[2] = rp
liquid_name = "switch_garbage/liquid_IC_%d" %seed_save
load_name = liquid_name + ".npy"
if solid_counter == 0:
np.save(liquid_name,y)
else:
y = np.load(load_name)
solid_counter = solid_counter + 1
return [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
else:
solid_switch = 1.0
solid_counter = 0.0
return [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
y[2] = rp
Mliq = 0.0
Mcrystal = 0.0
[FH2O,H2O_Pressure_surface] = H2O_partition_function( y[1],Mliq,Mcrystal,rp,g,kH2O,y[24])
[FCO2,CO2_Pressure_surface] = CO2_partition_function( y[12],Mliq,Mcrystal,rp,g,kCO2,y[24])
Tsolidus_unaltered = sol_liq(rp,g,pm,rp,float(H2O_Pressure_surface+CO2_Pressure_surface+1e5),float(0*y[0]/Mantle_mass0))
Tsolidus_Pmod = Tsolidus_unaltered + depletion_fraction * 600
if PltTech == 1:
deltardck = 0
else:
deltardck = y[26]
Tsolidus_visc = sol_liq(rp-deltardck,g,pm,rp,float(H2O_Pressure_surface+CO2_Pressure_surface+1e5),0.0)
visc = viscosity_fun(y[7],pm,visc_offset,y[8],float(Tsolidus_visc))
initialize_fast = np.min([y[8],y[7]+TMoffset,y[7]+TMoffset-1])
AB = AB_fun(float(y[8]),H2O_Pressure_surface,float(y[1]+y[4]+y[12]),albedoC,albedoH)
ASR_input = float((1-AB)*ASR_new_fun(t0))
if ASR_input < min_ASR:
ASR_input = min_ASR
Te_ar = (ASR_input/5.67e-8)**0.25
Te_input = Te_ar*(0.5**0.25)
if Te_input > 350:
Te_input = 350.0
if Te_input < min_Te:
Te_input = min_Te
initialize_fast = np.array(initialize_fast)
ace1 = nelder_mead(funTs_general_stag2, x0=initialize_fast, bounds=np.array([[10.0], [4000.0]]).T, args = (float(y[7]),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22]),float(y[25]),PltTech), tol_f=0.0001,tol_x=0.0001, max_iter=1000)
SurfT = ace1.x[0]
if abs(ace1.fun) > 1.0:#0.1:
ace2 = optimize.minimize(funTs_general_stag,x0=initialize_fast,args = (float(y[7]),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22]),float(y[25]),PltTech),method='Nelder-Mead',bounds = ((10,4000),))
if (abs(ace2.fun) < abs(ace1.fun))and(abs(ace2.fun)<1.0):
SurfT = ace2.x[0]
else:
ace3= optimize.minimize(funTs_general_stag,x0=y[8],args = (float(y[7]),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22]),float(y[25]),PltTech),method='L-BFGS-B',bounds = ((600,2000),))
if (abs(ace3.fun) < abs(ace1.fun))and(abs(ace3.fun)<1.0):
SurfT = ace3.x[0]
else:
ace4= optimize.minimize(funTs_general_stag,x0=y[8],args = (float(y[7]),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22]),float(y[25]),PltTech),method='L-BFGS-B',bounds = ((10,600),))
if (abs(ace4.fun) < abs(ace1.fun))and(abs(ace4.fun)<1.0):
SurfT = ace4.x[0]
if abs(ace1.fun) > 1.0:
differ_ace = 10.0
counter = 0
initialize_fast = np.array(initialize_fast)
budget_differ = 10.0
aceCOMP = nelder_mead(funTs_general_stag2, x0=initialize_fast, bounds=np.array([[10.0], [4000.0]]).T, args = (float(y[7]),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22]),float(y[25]),PltTech), tol_f=1e-10,tol_x=1e-10, max_iter=1000)
if abs(aceCOMP.fun)<ace1.fun:
SurfT = aceCOMP.x[0]
budget_differ = aceCOMP.fun
budget_counter = 0
while abs(budget_differ) > 0.1:
budget_counter = budget_counter + 1
ran_num = np.min([150 + 2000*np.random.uniform(),3999])
ran_numar=np.array(ran_num)
aceCOMP2 = nelder_mead(funTs_general_stag2, x0=ran_numar, bounds=np.array([[10.0], [4000.0]]).T, args = (float(y[7]),ll,visc,beta,Te_input,ASR_input,H2O_Pressure_surface,CO2_Pressure_surface,ocean_CO3,float(y[24]),float(y[22]),float(y[25]),PltTech), tol_f=1e-10,tol_x=1e-10, max_iter=2000)
if abs(aceCOMP2.fun)<abs(budget_differ):
SurfT = aceCOMP2.x[0]
budget_differ = aceCOMP2.fun
if budget_counter >= 10:
budget_differ = -2e-3
else:
SurfT = ace1.x[0]
T_diff = y[7]+TMoffset-SurfT
Qr = Qr_cofactor*qr(t0,Start_time,heatscale,K_over_U)
Qcore = Qr_cofactor*np.exp(-(t0/(1e9*365*24*60*60)-4.5)/7.0)*core_flow_coefficient
if PltTech ==1:
Vmantle = Vmantle0
Ra = np.max([-1e15,alpha * g * T_diff * ll**3 / (kappa * visc)])
qm = (k/ll) * T_diff * (Ra/Racr)**beta
volc_term = 0.0
dmelt_eq = y[16]
dcrust = y[26]
deltac = dmelt_eq
Vcrust_eq = (rp**3 - (rp-dmelt_eq)**3)*4*np.pi/3.0
dVcrust_dt = 0.0
y[27] = Vcrust_eq
ddelta_dt = 0.0
y[26] = y[16]
dQmantle_dt = 0.0 #Qr # per unit volume
dQcrust_dt = 0.0 #
else:
Vmantle = Vmantle0 - y[27]
delta = y[26]
Rai =alpha * g * T_diff* ll**3 / (kappa * visc)
theta = 300000*T_diff/(8.314*float(y[7])**2)
qm = 0.5 * k/ll * T_diff * theta**(-4.0/3.0) * Rai**(1./3.)
deltaTm = y[7] - SurfT ## adiabatic loss if stagnant lid
Vcrust = y[27]
deltac = rp - (rp**3 - 3*Vcrust/(4.0*np.pi))**(1./3.)
T1 = float(y[7]) - 2.5 * 8.314*float(y[7])**2/300000.0
fm = float(y[25])
term1 = (SurfT * (delta - deltac) + T1 * deltac)/delta
xm = Qr * pm * y[28]*Vmantle0/Vmantle # W/kg rock * kg rock/m3 * Vmantle0/Vmantle
xc = Qr * pm * y[29]*Vmantle0/Vcrust
term2 = (xc*deltac**2 * (delta-deltac) + xm * deltac * (delta - deltac)**2 )/ (2*k*delta)
Tc = term1 + term2
if (abs(Tc-T1)>1):
kdtdz_at_Rp_minus_delta = - k*(T1 - Tc)/(delta - deltac) + xm*0.5*(delta - deltac)
else:
kdtdz_at_Rp_minus_delta = - k*(T1 - SurfT)/deltac + xc*0.5 * deltac
ddelta_dt = (-qm - kdtdz_at_Rp_minus_delta) / (pm*cp*(y[7] - T1))
dVcrust_dt = fm - (fm - 4*np.pi*(rp - delta)**2 * np.min([0.0,ddelta_dt])) * (np.tanh((deltac -delta)*20) + 1)
Ra = Rai
volc_term = float(y[25])*3000*(cp*deltaTm + delHf)
y[32] = qm * 4.0 * np.pi*rp**2
y[33] = volc_term
[OLR,newpCO2,ocean_pH,ALK,Mass_oceans_crude,DIC_check] = correction(float(SurfT),float(Te_input),H2O_Pressure_surface,CO2_Pressure_surface,rp,g,ocean_CO3,1e5,float(y[24]),float(y[22]))
CO2_Pressure_surface = newpCO2
ASR = ASR_input
heat_atm = OLR - ASR
delta_u = k * (y[7]+TMoffset - SurfT) / qm
numerator = - 4*np.pi*(rp**2)*qm + (Qr * pm * y[28]*Vmantle0 ) + Qcore - volc_term
denominator = pm * cp * Vmantle
thermal_cond = 2 # W/m2/K
qc = thermal_cond * (SurfT - y[7]-TMoffset) / (rp-rc)
if SurfT > y[7]+TMoffset:
numerator = 4*np.pi*(rp**2)*qc + (Qr * pm * y[28]*Vmantle0 ) + Qcore
denominator = pm * cp * Vmantle
y[32] = -qc * 4.0 * np.pi*rp**2
y[33] = 0
ddelta_dt = 0.0
dVcrust_dt = 0.0
rs_term = 0.0
dy7_dt = numerator/denominator
[dK40mantle_dt,dAr40mantle_dt,dK40lid_dt,dAr40atmo_dt] =[ 0.0,0.0,0.0,0.0]
[dU238_mantle_dt,dU235_mantle_dt,dTh_mantle_dt,dU238_lid_dt,dU235_lid_dt,dTh_lid_dt,dHe_mantle_dt,dHe_atmo_dt] = [0,0,0,0,0,0,0,0]
y[8] = SurfT
liquid_switch = 1.0
liquid_switch_worked = 0.0
####################################################################################
## end magma ocean/solid mantle portion. The rest of the code applies to both phases
####################################################################################
toc1 = time.time()
y[30] = (Qr * pm * y[28]*Vmantle0)
y[31] = (Qr * pm * y[29]*Vmantle0 )
y[9] = OLR
y[10] = ASR
if y[8]<=y[7]+TMoffset:
y[11] = qm
else:
y[11] = -qc
if (2>1): #same thing seems to work
dy2_dt = rs_term * dy7_dt
drs_dt = dy2_dt
rs = np.min([y[2],rp])
[FH2O,H2O_Pressure_surface] = H2O_partition_function( y[1],Mliq,Mcrystal,rp,g,kH2O,y[24])
if y[1]< 0:
H2O_Pressure_surface = 0.0
FH2O = 0.0
water_frac = my_water_frac(float(y[8]),Te_input,H2O_Pressure_surface,CO2_Pressure_surface)
atmo_H2O = np.max([H2O_Pressure_surface*water_frac,0.0])
if solid_switch == 0:
T_partition = np.max([y[7],Tsolidus+0.01])
[XFeO,XFe2O3,fO2,F_FeO1_5,F_FeO] = solve_fO2_F_redo(y[4],H2O_Pressure_surface,T_partition,Total_Fe_mol_fraction,Mliq,rs,rp,y[24]) #
else:
[XFeO,XFe2O3,fO2,F_FeO1_5,F_FeO] = solve_fO2_F_redo(y[4],H2O_Pressure_surface,y[8],Total_Fe_mol_fraction,Mliq,rs,rp,y[24]) #
fO2_pos = np.max([0,fO2])
Pressure_surface = fO2_pos +atmo_H2O + CO2_Pressure_surface + 1e5
#####################################
### OPTIONAL CO2 dep Tstrat (AND [51] below)
#CO2_upper_atmosphere = y[51]
#Te_input_escape = 214 - 44*CO2_upper_atmosphere
##############################
frac_h2o_upper = my_fH2O(float(SurfT),Te_input_escape,H2O_Pressure_surface,CO2_Pressure_surface)
frac_h2o_upper = np.min([atmo_H2O / Pressure_surface, frac_h2o_upper]) #fudged this for 5x less escape
if (H2O_Pressure_surface<1e-5)and(frac_h2o_upper<1e-9): # check to see if it matters
frac_h2o_upper = 0.0
atmo_H2O = 0.0
H2O_Pressure_surface = 0.0
#######################
## Atmosphsic escape calculations
## diffusion limited escape:
fCO2_p = (1- frac_h2o_upper)*CO2_Pressure_surface / (CO2_Pressure_surface+fO2_pos+1e5)
y[51] = fCO2_p
fO2_p = (1- frac_h2o_upper)*fO2_pos / (CO2_Pressure_surface+fO2_pos+1e5)
fN2_p = (1- frac_h2o_upper)*1e5 / (CO2_Pressure_surface+fO2_pos+1e5)
mol_diff_H2O_flux = better_atomic_diffusion(frac_h2o_upper,Te_input_escape,g,fCO2_p,fO2_p,fN2_p) #mol H2O/m2/s I think
#XUV-driven escape
XH_upper = 2*frac_h2o_upper / (3*frac_h2o_upper + 2*fO2_p + fCO2_p + fN2_p) ## assumes CO2 and N2 don't dissociate, hmm
## now calculate atomic ratios and do XUV limited
XH = 2*frac_h2o_upper / (3*frac_h2o_upper + 2*fO2_p + fCO2_p + fN2_p)
XO = (2*fO2_p+frac_h2o_upper) / (3*frac_h2o_upper + 2*fO2_p + fCO2_p + fN2_p)
XC = fCO2_p / (3*frac_h2o_upper + 2*fO2_p + fCO2_p + fN2_p)
true_epsilon = find_epsilon(Te_input,RE,ME,float(AbsXUV(t0)), XO, XH, XC,epsilon,mix_epsilon)
#true_epsilon = find_epsilon(1000.0,RE,ME,float(AbsXUV(t0)), XO, XH, XC,epsilon,mix_epsilon)
if (XC < 0)or(y[12]<1e10):
XC = 0.0
[mH_Odert3,mO_Odert3,mOdert_build3,mC_Odert3] = Odert_three(Te_input_escape,RE,ME,true_epsilon,float(AbsXUV(t0)), XO, XH, XC) #kg/m2 /s
#[mH_Odert3,mO_Odert3,mOdert_build3,mC_Odert3] = Odert_three(1000.0,RE,ME,true_epsilon,float(AbsXUV(t0)), XO, XH, XC) #kg/m2 /s
numH = ( mH_Odert3 / 0.001 ) # mol H / m2/ s
numO = ( mO_Odert3 / 0.016 ) # mol O / m2/ s
#numC = ( mC_Odert3 / 0.012 ) # mol C / m2/ s
numC = ( mC_Odert3 / 0.044 ) # mol CO2 / m2/ s
if 2*mol_diff_H2O_flux> numH: ## if diffusion limit exceeds XUV-driven, shift diffusion-limit downward
mol_diff_H2O_flux = 0.5*np.copy(numH)
## Combined escape flux, weighted by H abundance:
w1 = mult*(2./3 - XH_upper)**4
w2 = XH_upper**4
Mixer_H = (w1*2*mol_diff_H2O_flux + w2 * numH ) / (w1+w2) # mol H / m2 /s
Mixer_O = (w1*0.0 + w2 * numO ) / (w1+w2)
Mixer_C = (w1*0.0 + w2 * numC ) / (w1+w2)
#Mixer_Build = 0.5* Mixer_H - Mixer_O + 2*Mixer_C # mol O / m2 /s C atmo drag affects redox
Mixer_Build = 0.5* Mixer_H - Mixer_O ## CO2 drag doesn't affect redox
## nominal
escape = 4*np.pi*rp**2 * Mixer_H*0.018/2 ## kg H2O /s
net_escape = 4*np.pi*rp**2 * Mixer_Build*0.016 ## kg O2/s
CO2_loss = 4*np.pi*rp**2 * Mixer_C*0.044 ## kg CO2/s
## ion escape test
### 2.0243e-14 * np.exp(8.05905*tgyr) # bar/ gyr
### 2.0243e-14 * np.exp(8.05905*tgyr) * (1e5*4*np.pi*rp**2)/(g*365*24*60*60*1e9) # kg/s
#tgyr = 4.6 - t0/(365*24*60*60*1e9)
#if tgyr >3.3:
# dOdt = 2.0243e-14 * np.exp(8.05905*tgyr) * (1e5*4*np.pi*rp**2)/(g*365*24*60*60*1e9) # kg/s
#else:
# dOdt = 2.0243e-14 * np.exp(8.05905*3.3) * (1e5*4*np.pi*rp**2)/(g*365*24*60*60*1e9)
#escape = 4*np.pi*rp**2 * Mixer_H*0.018/2
#net_escape = 4*np.pi*rp**2 * Mixer_Build*0.016 - dOdt
#CO2_loss = 4*np.pi*rp**2 * Mixer_C*0.044 ## kg CO2/s
## end ion escape test
if y[50]>0:
He_mixing_ratio = (y[50]/0.004) / ((Pressure_surface*(4*np.pi*rp**2)/g)/y[24])
else:
He_mixing_ratio = 0.0
He_escape_loss = 9.6e5/(365*24*60*60) * (AbsXUV(t0)/AbsXUV(4.5*1e9*365*24*60*60)) * (He_mixing_ratio/10e-6)
# done with escape calculations
#######################
## Find ocean depth and land fraction:
Ocean_depth = (0.018/y[24]) * (1-water_frac) * H2O_Pressure_surface / (g*1000) ## max ocean depth continents 11.4 * gEarth/gplanet (COwan and Abbot 2014)
Max_depth = 11400 * (9.8 / g)
if Ocean_depth > Max_depth:
Linear_Ocean_fraction = 1.0
Ocean_fraction = 1.0
else:
Linear_Ocean_fraction = (Ocean_depth/Max_depth) ## Approximation to Earth hypsometric curve
Ocean_fraction = (Ocean_depth/Max_depth)**0.25
## Melt and crustal oxidation variables:
actual_phi_surf_melt = 0.0
Va = 0.0
F_CO2 = 0.0
F_H2O = 0.0
O2_consumption = 0.0
OG_O2_consumption =0.0
Abs_crustal_prod = 0.0
crustal_depth = rp - (rp**3 - y[27]*3./(4.*np.pi))**(1./3.)
Poverburd = fO2_pos +H2O_Pressure_surface + CO2_Pressure_surface + 1e5
if PltTech == 1:
deltardck = 0
else:
deltardck = y[26]
Tsolidus_unaltered = sol_liq(rp-deltardck,g,pm,rp,float(Poverburd),float(0*y[0]/Mantle_mass0))
Tsolidus_Pmod = Tsolidus_unaltered + depletion_fraction * 600
y[34] = 0.0
dQmantle_dt = 0.0
dQcrust_dt = 0.0 #
ddt_depletion_fraction = 0.0
## If solid mantle, calculate crutal production and melt volume:
if (y[2]>=rp)and(Ra>0)and(y[8]<=Tsolidus):
T_for_melting = float(y[7])
if T_for_melting <= Tsolidus_Pmod: # no melting if potential temperature is below solidus
[actual_phi_surf_melt, Va] = [0.0, 0.0]
Melt_volume = 0.0
Abs_crustal_prod = 0.0
else: #mantle temperature above solidus, calculate melt production
mantleH2Ofrac = float(0*y[0]/Mantle_mass)
rdck = optimize.minimize(find_r2,x0=float(y[2]),args = (T_for_melting,alpha,g,cp,pm,rp,float(Poverburd),mantleH2Ofrac,deltardck,depletion_fraction))
rad_check = float(rdck.x[0]) # find radius where partial melting begins
y[34] = rp - rad_check
if rad_check>rp:
rad_check = rp
if PltTech == 1:
rlid = rp
#calculate melt production
[actual_phi_surf_melt,actual_visc,Va] = temp_meltfrac2(0.99998*rad_check,rp,alpha,pm,T_for_melting,cp,g,Poverburd,mantleH2Ofrac,rlid,depletion_fraction)
crustal_depth = rp - (rp**3 - actual_phi_surf_melt * Va*3./(4.*np.pi))**(1./3.) #calulate crustal depth (depth of melt)
Ra_melt =alpha * g * (y[7]+TMoffset - SurfT) * ll**3 / (kappa * visc)
Q = (k/ll) * (y[7]+TMoffset - SurfT) * (Ra_melt/Racr)**beta
Aoc = 4*np.pi*rp**2
Melt_volume = (Q*4*np.pi*rp**2)**2/(2*k*(y[7]+TMoffset - SurfT))**2 * (np.pi*kappa)/(Aoc) * crustal_depth
Abs_crustal_prod = Melt_volume
ddt_depletion_fraction = 0
elif (PltTech <1)and(0.99998*rad_check < rp-y[26]):
rlid = rp - y[26]
[actual_phi_surf_melt,actual_visc,Va] = temp_meltfrac2(0.99998*rad_check,rp,alpha,pm,T_for_melting,cp,g,Poverburd,mantleH2Ofrac,rlid,depletion_fraction )
crustal_depth = rp - (rp**3 - y[27]*3./(4.*np.pi))**(1./3.)
dm = rp - 0.99998*rad_check #where melting begins
velocity = 0.05 * kappa /(rp - rc) * (Rai/theta)**(2.0/3.0)
Melt_volume = 17.8 * np.pi * rp * velocity * (dm - y[26]) * actual_phi_surf_melt
ddt_depletion_fraction = 0.0 #Melt_volume*(1 - actual_phi_surf_melt) * ( 1 - depletion_fraction) / Vmantle # Uncomment for mantle depletion test
Abs_crustal_prod = 0.1*Melt_volume
elif (PltTech <1) and (0.99998*rad_check >= rp-y[26]):
[actual_phi_surf_melt,actual_visc,Va] = [0,0,0]
Melt_volume = 0.0
Abs_crustal_prod = 0.0
crustal_depth = rp - (rp**3 - y[27]*3./(4.*np.pi))**(1./3.)
ddt_depletion_fraction = 0.0
Dradio = 0.002
D_K_potassium = 1.0