-
Notifications
You must be signed in to change notification settings - Fork 0
/
solliq.py
1163 lines (1061 loc) · 99.5 KB
/
solliq.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
# coding: utf-8
# ## solliq - Solidus and liquidus parameterizations for geomaterials
# #### Thomas Ruedas (Thomas dot Ruedas at dlr dot de) - solliq v.1.3, 14/10/2018
# The solidi and liquidi of geomaterials are phase boundaries and depend on pressure $p$, temperature $T$, and composition. Solidi and liquidi of peridotite and basalt/eclogite (both anhydrous) as well as iron and iron-sulfur alloys are given here as parameterizations in the form $T(p)$, possibly with additional terms that account for compositional variations. $T$ is always given in K, $p$ in GPa.
#
# #### Note on citation
# Some of the parameterizations implemented here are published in the scientific literature:
# - [Terrestrial peridotite upper-mantle solidus](#Tsol_E): Hirschmann (2000), Hirschmann et al. (2009), with the modifications described in the corresponding section below
# - [Martian peridotite solidus](#Tsol_M): Ruedas & Breuer (2017)
# - [CMAS solidus](#Tsol_CMAS): Ruedas & Breuer (2018)
# - [Fe and alkali corrections](#Fealk): Ruedas & Breuer (2018)
# - [Chondritic lower-mantle solidus](#Tsol_chon): Andrault et al. (2018)
#
# If you use this software for your work, please acknowledge it and also cite the corresponding literature where applicable.
#
# This software is released under the GNU General Public License v.3.
#
# #### Contents
# - [Fitting approach](#fit)
# - [Peridotite solidus](#Tsol_per)
# - [Terrestrial peridotite](#Tsol_E)
# - [CMAS solidus](#Tsol_CMAS)
# - [Martian solidus](#Tsol_M)
# - [Fe and alkali corrections](#Fealk)
# - [Chondritic solidus](#Tsol_chon)
# - [Peridotite liquidus](#Tliq_per)
# - [Basalt/eclogite solidus](#Tsol_bas)
# - [Basalt/eclogite liquidus](#Tliq_bas)
# - [Iron and binary iron alloys](#alloy)
# - [Iron](#Fe)
# - [FeS](#FeS)
# - [Eutectic](#eutect)
# - [Interpolated melting curves](#intpol)
# - [Main program (user interface)](#main)
# - [Version history](#hist)
# - [References](#ref)
# ### <a id='fit'></a>Fitting approach
# The parameterizations of the melting curves are constructed by performing least-squares fits of the available experimental data to empirical model functions, in particular polynomials of first to third degree or to the Simon-Glatzel equation (Simon & Glatzel, 1929), which is usually given as
# \begin{equation}
# T_\mathrm{m}(p)=T_\mathrm{ref}\left(\frac{p-p_\mathrm{ref}}{a'}+1\right)^\frac{1}{c'},
# \end{equation}
# where $p_\mathrm{ref}$ and $T_\mathrm{ref}$ are a reference pressure and temperature often chosen to correspond to the melting point at surface (zero) pressure. The equation is used here in the modified form
# \begin{equation}
# T_\mathrm{m}(p)=\frac{T_\mathrm{ref}}{a'^\frac{1}{c'}}(p-p_\mathrm{ref}+a')^\frac{1}{c'}=a(p+b)^c,
# \end{equation}
# with $a=T_\mathrm{ref}/a'^\frac{1}{c'}$, $b=a'-p_\mathrm{ref}$, and $c=\frac{1}{c'}$. For each dataset, several model functions were tried, and the one with the least rms residual was usually chosen if it met certain requirements. The most important requirements reflect experimental experience and are that the curve should be monotonically increasing and have no inflexion points in the $p$ range of interest and that its slope should decrease with increasing $p$. The functions are generally not meant to be used for extrapolation much beyond the range of data coverage.
#
# If the data warranted it, the possible effects of phase transitions of the solid phase were taken into account by carrying out piecewise fits; such solid-solid phase transitions may manifest themselves for instance in the form of kinks in the melting curve, although sharp kinks are expected in chemically pure systems rather than solid-solution systems such as several of those of interest here. In such piecewise fits, continuity conditions were imposed on solidus fits to data from natural compositions that included solid solutions where necessary; the practice to allow jumps in the solidus followed in older versions is abandoned now, because smooth transitions are expected in all systems that feature solid solutions. The most important example is the transition from a terrestrial upper-mantle mineral assemblage to a lower-mantle assemblage.
# ### <a id='Tsol_per'></a>Peridotite solidus
# #### <a id='Tsol_E'></a>Terrestrial peridotite
# Data for terrestrial peridotite are available for the entire depth range of the Earth's mantle. The following parameterization is based on Hirschmann (2000, tab.2 and pers.comm.) and Hirschmann et al. (2009) for the upper mantle, and the downward shift by 35 K suggested by Katz et al. (2003) is applied. In the lower mantle, the linear Hirschmann et al. fit is replaced with a curved model function fitted to some of the data also used by Hirschmann et al. (2009), supplemented by additional data from Hirose & Fei (2002), Fiquet et al. (2010), and Zerr et al. (1997,1998), covering 22.5 GPa $< p <$ 122 GPa.
# \begin{equation}
# T_\mathrm{s,E}(p)=
# \begin{cases}
# 1358.81061+132.899012p-5.1404654p^2&p<10\,\mathrm{GPa}\\
# -1.092(p-10)^2+32.39(p-10)+2173.15&10\,\mathrm{GPa}\leq p < 23.5\,\mathrm{GPa}\\
# -1647.15-5.85p+1329.126\ln(p)&\text{else.}
# \end{cases}
# \end{equation}
# A Simon-Glatzel-type fit could not be carried out for the lower mantle data with the constraint to match the Hirschmann et al. (2009) curve at 23.5 GPa. The high-$p$ study on a peridotitic composition by Nomura et al. (2014) was excluded, because they used a synthetic composition (pyrolite) and arrived at a substantially lower solidus; futhermore, their data have strong scatter and provide only rather wide brackets on the solidus temperature.
#
# The lower-mantle solidus permits plausible extrapolation to higher $p$, e.g., for application in exoplanets, but it must be kept in mind that a transition of bridgmanite to postperovskite is expected somewhere between 130 and 170 GPa, and there are no experimental constraints on the solidus in a postperovskite lithology.
# In[ ]:
import numpy as np
def Tsol_Earth(p):
if p < 10:
Tsol=1358.81061+(132.899012-5.1404654*p)*p # w/ Katz shift
elif p >= 10 and p < 23.5:
pr=p-10.
Tsol=(-1.092*pr+32.39)*pr+2173.15 # w/ Katz shift
else:
Tsol=-1647.15020384-5.85124635*p+1329.1263647*np.log(p)
return Tsol
# #### <a id='Tsol_CMAS'></a>CMAS solidus
# The CaO-MgO-Al<sub>2</sub>O<sub>3</sub>-SiO<sub>2</sub> (CMAS) system is the simplest oxide system capable of reproducing all four phases of a fertile lherzolite, i.e., olivine (forsterite Mg<sub>2</sub>SiO<sub>4</sub>), orthopyroxene (enstatite MgSiO<sub>3</sub>+wollastonite CaSiO<sub>3</sub>), clinopyroxene (diopside CaMgSi<sub>2</sub>O<sub>6</sub>), and either plagioclase (anorthite CaAl<sub>2</sub>Si<sub>2</sub>O<sub>8</sub>), spinel (MgAl<sub>2</sub>O<sub>4</sub>), or garnet (pyrope Mg<sub>3</sub>Al<sub>2</sub>Si<sub>3</sub>O<sub>12</sub>+grossular Ca<sub>3</sub>Al<sub>2</sub>Si<sub>3</sub>O<sub>12</sub>). It can also serve as a proxy for the presumably very Fe-poor peridotite of Mercury. Data are available for $p$ up to 21.5 GPa (Asahara et al., 1998; Asahara & Ohtani, 2001; Gudfinnsson & Presnall, 1996, 2000; Herzberg et al., 1990; Klemme & O'Neill, 2000; Litasov & Ohtani, 2002; Liu & O'Neill, 2004; Milholland & Presnall, 1998; Presnall et al., 1979) and have been fitted with a uniform polynomial, i.e., without considering phase transitions that may introduce kinks (which would probably be rather sharp in this case due to the absence of solid solutions in some of the phases):
# \begin{equation}
# T_\mathrm{s,CMAS}(p)=1477.54+139.391p-7.7545p^2+0.160258p^3
# \end{equation}
# (Ruedas & Breuer, 2018, App.B). For $p>23.5$ GPa, which would correspond to the lower mantle, we use the terrestrial peridotite solidus, tentatively shifted up by about 139 K to ensure continuity, but there are no experimental constraints.
# In[ ]:
def Tsol_CMAS(p):
if p < 23.5:
Tsol=1477.54+(139.391-(7.7545-0.160258*p)*p)*p
else:
Tsol=Tsol_Earth(p)+139.2162
return Tsol
# #### <a id='Tsol_M'></a>Martian solidus
# The peridotite of the martian mantle has a different composition than its terrestrial counterpart, with the substantially lower Mg# being the most conspicuous difference. Ruedas & Breuer (2017, App. A) fitted experiments on martian model compositions (Bertka & Holloway, 1994; Collinet et al., 2015; Matsukage et al., 2013; Schmerr et al., 2001) to polynomials to produce the parameterization
# \begin{equation}
# T_\mathrm{s,M}(p)=
# \begin{cases}
# 0.118912p^3-6.37695p^2+130.33p+1340.38&\text{for $p<23$ GPa}\\
# 62.5p+975&\mathrm{else.}
# \end{cases}
# \end{equation}
# With the exception of the Collinet et al. data, all data used for the fit have been shifted downward by 30 K to correct for the absence of K in the sample, using an estimate on the effect of K on the solidus derived from Wang & Takahashi (2000). As Collinet et al. did not provide subsolidus data, the melting degree-temperature relations from their experiments were extrapolated linearly at each pressure, and the resulting solidus estimate was weighted doubly. Note that the function for $p>23$ GPa is constrained by only four points and only to 25 GPa, and so extrapolation to much higher $p$ is strongly discouraged. For Mars itself, this is not a problem, because Mars is not expected to have a perovskitic basal layer, but it would be a concern for exoplanets that are larger than Mars and have an iron-rich mantle.
# In[ ]:
def Tsol_Mars(p):
if p < 23:
Tsol=((0.118912*p-6.37695)*p+130.33)*p+1340.38
else:
Tsol=62.5*p+975
return Tsol
# #### <a id='Fealk'></a>Fe and alkali corrections
# For peridotite with compositions not covered by experimental data, we can try to interpolate between better constrained solidi. The reference is the terrestrial solidus, and interpolation is carried out by linear interpolation based on Mg# between it and either the martian or the CMAS solidus corrected for the difference with Earth in their alkali (Na+K) contents. For these corrections, the oxide composition (specifically, the MgO, FeO, Na<sub>2</sub>O, and K<sub>2</sub>O content) of the target material needs to be known; the terrestrial and martian reference oxide compositions are taken from Palme & O'Neill (2014, Tabs. 3,4) and Taylor & McLennan (2009, Table 5.4), respectively. An alkali effect $\mathrm{d}T_\mathrm{s}/\mathrm{d}X_\mathrm{Na+K}=-14951.52$ K has been estimated from values from Hirschmann (2000) and Wang & Takahashi (2000). The interpolated solidus is then
# \begin{equation}
# T_\mathrm{s}(p)=(1-w_\mathrm{Mg\#})T_\mathrm{s,E}+w_\mathrm{Mg\#}T_\mathrm{sr}+
# \frac{\mathrm{d}T_\mathrm{s}}{\mathrm{d}X_\mathrm{Na+K}}(X_\mathrm{Na+K}-X_\mathrm{Na+K,E}),
# \end{equation}
# where
# \begin{align}
# w_\mathrm{Mg\#}&=\frac{\mathrm{Mg\#}_\mathrm{E}-\mathrm{Mg\#}}{\mathrm{Mg\#}_\mathrm{E}-\mathrm{Mg\#}_2}\\
# T_\mathrm{sr}&=T_\mathrm{s,2}-\frac{\mathrm{d}T_\mathrm{s}}{\mathrm{d}X_\mathrm{Na+K}}(X_\mathrm{Na+K,2}-X_\mathrm{Na+K,E}).
# \end{align}
# The subscripts E and 2 indicate the Earth and the other reference system (i.e., Mars for $\mathrm{Mg\#}<\mathrm{Mg\#}_\mathrm{E}$ and CMAS otherwise), respectively. $T_\mathrm{sr}$ is hence the solidus of the other reference system with the alkali effect reduced to that of the Earth.
# In[ ]:
# oxide contents of reference bulk silicate planet compositions
ox_E={'MgO': 36.77e-2,'FeO': 8.1e-2,'Na2O': 3.5e-3,'K2O': 3e-4} # Earth
ox_M={'MgO': 30.2e-2,'FeO': 17.9e-2,'Na2O': 5e-3,'K2O': 4e-4} # Mars
def Tsol_intp(p,ox):
global Mg0,Mg0_E
uMgO=u['Mg']+u['O']
uFeO=u['Fe']+u['O']
el1=ox_E['MgO']/uMgO
el2=ox_E['FeO']/uFeO
Mg0_E=el1/(el1+el2)
xwNaK_E=ox_E['Na2O']+ox_E['K2O']
dTsdNaK=-14951.52 # K per mass frac.
# compositional variables
el1=ox['MgO']/uMgO
el2=ox['FeO']/uFeO
Mg0=el1/(el1+el2)
# interpolation
TsE=Tsol_Earth(p)
if Mg0 < Mg0_E:
# iron-rich, interpolate between Earth and Mars solidi
el1=ox_M['MgO']/uMgO
el2=ox_M['FeO']/uFeO
w_dMg0=(Mg0_E-Mg0)/(Mg0_E-el1/(el1+el2))
TsM=Tsol_Mars(p)-dTsdNaK*(ox_M['Na2O']+ox_M['K2O']-xwNaK_E)
else:
# iron-poor, interpolate between Earth and CMAS solidi
w_dMg0=(Mg0-Mg0_E)/(1-Mg0_E)
TsM=Tsol_CMAS(p)+dTsdNaK*xwNaK_E
# weighted linear interpolation plus alkali correction relative to terrestrial
return (1.-w_dMg0)*TsE+w_dMg0*TsM+dTsdNaK*(ox['Na2O']+ox['K2O']-xwNaK_E)
# #### <a id='Tsol_chon'></a>Chondritic solidus
# Some melting experiments have also been conducted on various natural and synthetic "chondritic" compositions. The fit implemented here follows Andrault et al. (2018) in treating the data for the (terrestrial) upper and lower mantle pressure range separately. In the latter, only the data from Andrault et al. (2011, synthetic CMASF) are available, as a single additional point from Ohtani & Sawamoto (1987) used a simplified (CMASF) composition and is probably not reliable. For lower pressures, we supplement the extensive dataset from Andrault et al. (2018, synthetic CMASFNCrT) with some points from Berthet et al. (2009, Indarch EH4), Ohtani (1987, synthetic CMASF), and Takahashi (1983, Yamato Y-74191 L3 chondrite) that seem to fit in well. Additional points from Agee & Draper (2004, natural and synthetic Homestead) and Cartier et al. (2014, Hvittis enstatitic chondrite EL6) were too far above the general trend to be considered. Data for carbonaceous chondrites were excluded. In summary, the implemented solidus is
# \begin{equation}
# T_\mathrm{s,ch}(p)=
# \begin{cases}
# 1337.159(p+0.9717)^{0.1663}&\text{for $p<27.5$ GPa}\\
# 520.735(p+9.6535)^{0.4155}&\mathrm{else,}
# \end{cases}
# \end{equation}
# whereby the latter is simply a rewritten form of the lower-mantle fit given by Andrault et al. (2018). This solidus has data coverage up to 140 GPa.
# In[ ]:
def Tsol_chon(p):
if p < 27.5:
Tsol=1337.159*(p+0.9717)**0.1663
else:
Tsol=520.735*(p+9.6535)**0.4155
return Tsol
# ### <a id='Tliq_per'></a>Peridotite liquidus
# For the liquidus, a distinction should be made between a system that retains all of the melt produced and thus has a constant bulk composition (batch melting) and a system in which the produced melt is removed immediately, leaving behind a solid system with a different bulk composition that changes continuously as melting progresses (fractional melting). For many applications in mantle dynamics, the latter case is considered more appropriate and implies, at least in a simplified treatment, that the liquidus to be used is that of the last phase to remain solid as the melting degree increases (cf. Iwamori et al., 1995, p.257). This "liquidus phase" varies with pressure. Based on the phase diagram for peridotite by Litasov & Ohtani (2007), the function switches through the succession forsterite-pyrope-periclase as $p$ increases; this phase diagram is derived from batch melting experiments, but due to lack of experimental work representing fractional melting, it is assumed that the same liquidus phases occur in fractionally melting mantle in approximately the same $p$ ranges. The switch to periclase is triggered by `wpc`, which represents the mass fraction of periclase, becoming non-zero; this is currently hard-wired to occur linearly between 23 and 24 GPa. The individual liquidus curves are a polynomial parameterization derived from experimental data for forsterite (Davis & England, 1964; Navrotsky et al., 1989; Ohtani & Kumazawa, 1981; Presnall & Walter, 1993; Richet et al., 1993) and Simon-Glatzel-type fits to experimental data for pyrope (Boyd & England, 1962; Irifune & Ohtani, 1986; Kudo & Ito, 1996; Ohtani et al., 1981; Shen & Lazor, 1995; Zhang & Herzberg, 1994) and to first-principles simulations of periclase by Alfè (2005), as no usable experimental data were available for the latter; experimentally determined melting curves of periclase still have large uncertainty and may vary by hundreds or even more than 1000 K at high $p$.
# \begin{equation}
# T_\mathrm{l,per}=
# \begin{cases}
# 2160.6+64.7109p-3.97463p^2+0.0957894p^3&p<13.244\,\text{GPa (forsterite)}\\
# 589.691(p+11.9076)^{0.453188}&p\geq 13.244\,\text{GPa and no periclase present (pyrope)}\\
# 1705.67(p+5.20372)^{0.315027}&\text{(periclase)}
# \end{cases}
# \end{equation}
# The melting curve of pyrope could in theory be extrapolated to normal pressure, but nominal melting temperature values at normal pressure have not been used for fitting, because pyrope melts incongruently under these conditions.
# In[ ]:
def Tliq_per(p):
pxfopy=13.244 # GPa, transition p from fo to py
if p > pxfopy:
wpc=min(max(0.,p-23),1.)
if wpc == 1:
Tliq_per=Tliqpc(p)
elif wpc > 0:
Tliq_per=(1.-wpc)*Tliqpy(p)+wpc*Tliqpc(p)
else:
Tliq_per=Tliqpy(p)
else:
Tliq_per=Tliqfo(p)
return Tliq_per
# forsterite melting curve
def Tliqfo(p):
pmaxfit=14.64
if p > pmaxfit:
print("WARNING: p>%.2f GPa, Tliq (forsterite) assumed const." % pmaxfit)
Tliqfo=2557.6807
else:
Tliqfo=2160.6+(64.7109-(3.97463-0.0957894*p)*p)*p
return Tliqfo
# pyrope melting curve
def Tliqpy(p):
return 589.691*(p+11.9076)**0.453188
# periclase melting curve
def Tliqpc(p):
return 1705.67*(p+5.20372)**0.315027
# For possible future use with peridotite or basalt/eclogite, a melting curve for bridgmanite was also constructed from experimental data (Zerr & Boehler, 1993):
# \begin{equation}
# T_\mathrm{l,br}(p)=-0.338678p^2+79.9414p+1169.02.
# \end{equation}
# It is not being used in the current version.
# In[ ]:
# bridgmanite melting curve
def Tliqbr(p):
return (-0.338678*p+79.9414)*p+1169.02
# For completeness and because it is frequently used, we also construct a "batch melting liquidus" for peridotite from experimental data on natural terrestrial compositions (Fiquet et al., 2010, KLB-1; Herzberg et al., 1990, KLB-1; Ito & Kennedy, 1967, KA 64-16; Ito & Takahashi, 1987, PHN1611; Scarfe & Takahashi, 1986, PHN1611; Takahashi, 1986, KLB-1; Takahashi & Scarfe, 1985, KLB-1; Takahashi et al., 1993, KLB-1; Trønnes & Frost, 2002, pyrolite/KLB-1):
# \begin{equation}
# T_\mathrm{lb,per}=
# \begin{cases}
# 2067.59+25.1327p-0.2455p^2&p<22.76\,\text{GPa}\\
# 2207.71(p-20.7372)^{0.182}&\text{else.}
# \end{cases}
# \end{equation}
# Some other datasets from the literature were not used, because they use synthetic, simplified or otherwise deviating compositions.
# In[ ]:
def Tliqb_per(p):
if p < 22.76:
Tm=2067.59+(25.1327-0.2455*p)*p
else:
Tm=2207.71*(p-20.7372)**0.182
return Tm
# For a chondritic composition, we use the datasets by Agee & Draper (2004, natural and synthetic Homestead), Andrault et al. (2011, synthetic CMASF), Berthet et al. (2009, Indarch EH4), Ohtani (1987, synthetic CMASF), and Takahashi (1983, Yamato Y-74191 L3 chondrite); again, carbonaceous chondrites were avoided. As the data do not warrant making a distinction between the upper and lower mantle pressure range, we fit the whole dataset with a Simon-Glatzel-type equation:
# \begin{equation}
# T_\mathrm{lb,ch}=310.317(p+29.102)^{0.5337}.
# \end{equation}
# In[ ]:
def Tliqb_chon(p):
return 310.317*(p+29.102)**0.5337
# ### <a id='Tsol_bas'></a>Basalt/eclogite solidus
# Basalt consists mostly of clinopyroxene and plagioclase. When the plagioclase transforms into garnet, the basalt becomes eclogite. Data for the solidus of basalt are more scarce, and so we interpolate between a few points measured from the (terrestrial) MORB basalt/eclogite phase diagram provided by Litasov & Ohtani (2007) for the stability field of basalt, which we extend to 3 GPa in order to achieve continuity with the eclogite parameterization, for which experimental data were fitted. When the eclogite transforms into a bridgmanite-bearing high-$p$ assemblage, the function shifts linearly to the high-$p$ melting curve over a $p$ interval from 26 to 30 GPa. The bound of 26 GPa separates the subsets of experimental data used for fitting the eclogite and the high-$p$ solidi and should therefore not be changed; the 30 GPa bound is somewhat arbitrary and may be changed if desired. The experimental data are from Andrault et al. (2014), Hirose et al. (1999), Hirose & Fei (2002), Pradhan et al. (2015), and Yasuda et al. (1994); from the Andrault et al. data, only those above ~56 GPa are used, because the lower-$p$ data lie significantly beneath those of the other studies.
# \begin{equation}
# T_\mathrm{s,bas}(p)=
# \begin{cases}
# 1.83261p^4+13.0994p^2+1340&p<3\,\text{GPa}\\
# -1.273p^2+78.676p+1380.92&3\,\text{GPa}\leq p<26\,\text{GPa}\\
# -0.0503p^2+19.277p+2165.09&p\geq 30\,\text{GPa.}
# \end{cases}
# \end{equation}
# In[ ]:
def Tsol_bas(p):
pgbu=30. # upper boundary of transition
Tsgt=lambda p: (-1.273*p+78.676)*p+1380.92
Tsbr=lambda p: (-0.0503*p+19.277)*p+2165.09
if p < 3:
# plagioclase stability field
p2=p*p
Ts=(1.83261*p2+13.0994)*p2+1340.
elif p >= 3 and p < 26:
# garnet/majorite stability field
Ts=Tsgt(p)
elif p >= 26 and p < pgbu:
# transition: linear shift from Tsgt to Tsbr
Ts=(Tsgt(p)*(pgbu-p)+Tsbr(p)*(p-26))/(pgbu-26)
else:
# perovskite stability field
Ts=Tsbr(p)
return Ts
# ### <a id='Tliq_bas'></a>Basalt/eclogite liquidus
# According to the phase diagram from Litasov & Ohtani (2007), the liquidus phase in basalt melting also changes with $p$. We switch between liquidus phases as $p$ increases and use again parameterized melting curves from experimental data, this time for forsterite as above, for diopside (Boyd & England, 1963; Gasparik, 1996; Richet & Bottinga, 1984; Shen & Lazor, 1995; Yoder, 1952), for pyrope as above, and for Ca-perovskite (Gasparik et al., 1994; Zerr et al., 1997).
# \begin{equation}
# T_\mathrm{l,bas}(p)=
# \begin{cases}
# 2160.6+64.7109p-3.97463p^2+0.0957894p^3&p<1\,\text{GPa (forsterite)}\\
# 1665.25+133.484p-9.02117p^2+0.243279p^3&1\,\mathrm{GPa}\leq p < 3.728\,\text{GPa (diopside)}\\
# 589.691(p+11.9076)^{0.453188}&3.728\,\mathrm{GPa}\leq p < 24\,\text{GPa (pyrope)}\\
# 394.881(p+11.1079)^{0.593651}&\text{else (Ca-perovskite).}
# \end{cases}
# \end{equation}
# This parameterization has some particularly large steps, e.g. at 1 GPa, which imply strong changes in the width of the partial melting interval and correspondingly strong variations with isobaric melt productivity in different stability fields.
# In[ ]:
def Tliq_bas(p):
if p < 1:
Tliq=Tliqfo(p)
elif p < 3.728:
Tliq=Tliqdi(p)
elif p < 24:
Tliq=Tliqpy(p)
else:
Tliq=TliqCapv(p)
return Tliq
# diopside melting curve
def Tliqdi(p):
return 1665.25+(133.484-(9.02117-0.243279*p)*p)*p
# Ca-perovskite melting curve
def TliqCapv(p):
return 394.881*(p+11.1079)**0.593651
# The phase diagram from Litasov & Ohtani (2007) and some experimental data from Hirose & Fei (2002) and Yasuda et al. (1994) are also used to construct a "batch melting liquidus" for basalt in a way similar to the solidus; for $p\gtrsim 22$ GPa, however, the melting curve of Ca-perovskite is used as eclogite transforms into a perovskitic assemblage and Ca-perovskite has been found to be the likely liquidus phase up to pressures corresponding to the terrestrial CMB:
# \begin{equation}
# T_\mathrm{lb,bas}(p)=
# \begin{cases}
# 21.5909p+1488&p<0.88\,\mathrm{GPa}\\
# -2.24727p^2+97.6753p+1422.79&0.88\,\mathrm{GPa}\leq p<3.12\,\mathrm{GPa}\\
# 1581.57(p-1.468)^{0.151}&3.12\,\mathrm{GPa}\leq p<22\,\mathrm{GPa}\\
# T_\mathrm{l,Capv}(p)&\text{else.}
# \end{cases}
# \end{equation}
# The interval between 0.88 and 3.12 GPa is derived from Litasov & Ohtani (2007) and serves as a link between the lowest-$p$ segment and the range constrained by experimental data. Note, however, that around 22 GPa this liquidus gets very close to the solidus and that it is not at all constrained at higher $p$; its use at $p\gtrsim 20$ GPa should probably be avoided.
# In[ ]:
def Tliqb_bas(p):
if p < 0.88:
# plagioclase stability field
Tliq=21.5909*p+1488.
elif p < 3.12:
# garnet stability field up to intersection with data fit
Tliq=(-2.24727*p+97.6753)*p+1422.79
elif p < 22:
# approx. garnet stability field, from exp. data
Tliq=1581.57*(p-1.468)**0.151
else:
# Ca-perovskite stability field
Tliq=TliqCapv(p)
return Tliq
# ### <a id='alloy'></a>Iron and binary iron alloys
# #### <a id='Fe'></a>Iron
# With increasing pressure, different Fe allotropes are on the solid side of the melting curve: body-centered cubic (bcc) $\delta$-Fe at low $p$, face-centered cubic (fcc) $\gamma$-Fe at intermediate $p$, and hexagonal (hcp) $\epsilon$-Fe at high $p$. However, the fcc-hcp-liquid triple point is somewhat uncertain, and the melting curve especially for hcp-Fe is still controversial, as there are two quite different trends to be found in published experimental data at $p\gtrsim 30$ GPa. For this reason, we consider two alternative melting curves, one following a fairly flat $p$ gradient in the fcc stability field and another following a much steeper gradient. At this point, the even more controversial high-$p$ body-centered phase hypothesized by some authors to replace hcp-Fe above a certain high pressure is not included, and hcp-Fe is assumed to be stable to the highest $p$ considered. Experimental melting curve data are from Liu & Bassett (1975), Strong et al. (1973), and Swartzendruber (1982) for $\delta$-Fe, from Ahrens et al. (2002), Anzellini et al. (2013), Aquilanti et al. (2015), Boehler (1993), Boehler & Ross (2007), Boehler et al. (2008), Buono & Walker (2011), Jackson et al. (2013), Liu & Bassett (1975), Ringwood & Hibberson (1990), Shen et al. (1998, 2004), Strong et al. (1973), Williams et al. (1987), and Zhang et al. (2016) for fcc-Fe, and from Anzellini et al. (2013), Aquilanti et al. (2015), Boehler (1993), Boehler & Ross (2007), Boehler et al. (2008), Ma et al. (2004), Murphy et al. (2011), Nguyen & Holmes (2004), Shen et al. (1998), and Tateno et al. (2010) for hcp-Fe.
#
# In order to construct the "flat" and "steep" fits, the data have been divided into two subsets that have some overlap at $p\lesssim 30$ GPa as well as at $p>300$ GPa, and some points were discarded. In the data subset used for describing the "flat" gradient, it seems useful to subdivide the $p$ into three segments according to the stable solid phase and construct peacewise fits, i.e.,
# \begin{equation}
# T_\mathrm{m,Fe}=
# \begin{cases}
# 1811+30.9724p&p\leq 8.044\,\text{GPa (bcc/$\delta$)}\\
# 1907.47+19.8675p-0.1103p^2&8.044\,\mathrm{GPa}<p\leq 87.06\,\text{GPa (fcc/$\gamma$)}\\
# 2130.55+7.2054p+0.00571047p^2&87.06\,\mathrm{GPa}<p\text{ (hcp/$\epsilon$).}
# \end{cases}
# \end{equation}
# As an exception to the rule, the hcp interval of this fit is slightly bent *upward* and should therefore not be used beyond 400 GPa; excluding the points at $p>300$ GPa from Tateno et al. (2010) from the fit would result in a markedly stronger upward bent, however, and is therefore further from the expected trend. The data subset used for the "steep" gradient has an overall appearance that seems to allow for a uniform fit, ignoring possible kinks at phase boundaries, which cannot be identified with any certainty:
# \begin{equation}
# T_\mathrm{m,Fe}=1811+24.7307p-0.0627041p^2+6.14455\cdot 10^{-5}p^3.
# \end{equation}
# The zero-pressure melting point of $\delta$-Fe has been fixed at 1811 K (Swartzendruber, 1982) in both cases.
# In[ ]:
def Tliq_Fe(p,mode):
if mode == 'f':
# flat gradient
if p <= 8.044:
# delta phase (bcc)
Tm=1811.+30.9724*p
elif p <= 87.06:
# gamma phase (fcc)
Tm=1907.47+(19.8675-0.1103*p)*p
else:
# epsilon phase (hcp)
Tm=2130.55+(7.2054+0.00571047*p)*p
else:
# steep gradient
Tm=1811+(24.7307-(0.0627041-6.14455e-5*p)*p)*p
return Tm
# #### <a id='FeS'></a>FeS
# Melting point data for (stoichometric) FeS are much scarcer. The experimental data from Anderson & Ahrens (1996) and Boehler (1992) are supplemented by some points read from a phase diagram by Ono et al. (2008); additional data from a study by Williams & Jeanloz (1990) were not used, because they show strong scatter and are inconsistent with the other data. The available data suggest that in spite of the occurrence of several phase transitions, a single Simon-Glatzel-type function can be fitted over the entire pressure range up to 360 GPa:
# \begin{equation}
# T_\mathrm{m,FeS}=638.5613(p+13.2844)^{0.3216},
# \end{equation}
# whereby the zero-pressure melting point has been fixed at 1467 K.
# In[ ]:
def Tliq_FeS(p):
return 638.5613*(p+13.2844)**0.3216
# #### <a id='eutectic'></a>Eutectic
# The melting curve of an Fe-S alloy with a composition between Fe and FeS is not a monotonic function of composition but has a minimum at the intermediate eutectic composition, at which both solid phases and the melt are in equilibrium and which is itself a function of pressure. Experimental work indicates that the molar fraction of S at the eutectic is a monotonically falling function of $p$ and can be represented as
# \begin{equation}
# \bar{X}_\mathrm{eut}(p)=0.925(p+8.444)^{-0.373},
# \end{equation}
# based on data by Brett & Bell (1969), Buono & Walker (2011), Campbell et al. (2007), Chen et al. (2008), Chudinovskikh & Boehler (2007), Fei et al. (1997, 2000), Kamada et al. (2010), Li et al. (2001), Morard et al. (2007, 2008), Mori et al. (2017), Ryzhenko & Kennedy (1973), Stewart et al. (2007), and Usselman (1975). Campbell et al. (2007) had suggested that $\bar{X}_\mathrm{eut}(p)$ approaches a value of 0.25 at $p\gtrsim 25$ GPa, but as this is in conflict with other data, a value of 0.25 was assumed to apply only to their data up to 36.6 GPa. This fit replaces the exponential fit by Ruedas et al. (2013, App.B), which works only for $p\lesssim 100$ GPa.
# In[ ]:
def Xeut(p):
return 0.925/(p+8.444)**0.373
# The corresponding temperature $T_\mathrm{eut}$ at the eutectic point is also a function of $p$. We follow Mori et al. (2017) in basing the model on a Simon-Glatzel-type function, but experimental work by Fei et al. (1997, 2000) has shown that the intricacies of the Fe-S phase diagram introduce a minimum somewhere between 10 and 15 GPa, which is here represented by a Gaussian-shaped additional term that is significant only in the pressure range up to $\sim 20$ GPa:
# \begin{gather}
# T_\mathrm{eut}(p)=255(p+8)^{\tau(p)}+600\exp(-0.012(p+2)^2)\\
# \tau(p)=0.49\left\lbrace 1-0.085\left[\frac{1}{2}-\frac{\arctan(0.007(250-p))}{\pi}\right]\right\rbrace;
# \end{gather}
# the monstrosity $\tau(p)$ is needed to reduce the exponent of the first term of $T_\mathrm{eut}$ at high $p$ and introduce a slight, smooth bend between approximately 100 and 150 GPa that ensures that the eutectic does not become greater than the pure-phase liquidi at high $p$. The experimental data for the (in this case visually performed) fit are from the same sources and from Boehler (1996), but in the studies predating Fei et al. (1997), data above 6 GPa had to be dismissed, because they did not account for the aforementioned complexities of the phase diagram under those conditions. As with $\bar{X}_\mathrm{eut}(p)$, this fit replaces the corresponding exponential fit by Ruedas et al. (2013, App.B), which overestimates $T_\mathrm{eut}$ for $p\gtrsim 150$ GPa.
# In[ ]:
def Teut(p):
from math import exp,atan,pi
fxp=0.49*(1-0.085*(0.5-atan(0.007*(250-p))/pi))
return 255*(p+8)**fxp+600*exp(-0.012*(p+2)**2)
# #### <a id='intpol'></a>Interpolated melting curves
# For alloys with an arbitrary binary composition $\bar{X}_\mathrm{S}$, the melting curve of the alloy can be determined from the melting curve of a pure phase. Starting from thermodynamical arguments by Stevenson (1981), Anderson (1998) proposed to determine the solidus depression at a given pressure caused by $N$ solutes in impure Fe by
# \begin{equation}
# \Delta T_\mathrm{m,Fe}(p)=T_\mathrm{m,Fe}(p) \sum\limits_{i=1}^N \ln(1-\bar{X}_i);
# \end{equation}
# in the case considered here, $N=1$ and the solute is S. However, the presence of a eutectic and the existence of FeS as a second endmember limit the applicability to at most the compositions on the iron-rich side of the eutectic. On the sulfur-rich side, one would instead consider the depression of the FeS melting curve by an amount $\bar{X}_\mathrm{Fe+}=\bar{X}_\mathrm{Fe}-\frac{1}{2}$ of Fe in addition to the amount of the stoichometric composition of FeS as a "solute":
# \begin{equation}
# \Delta T_\mathrm{m,FeS}(p)=T_\mathrm{m,FeS}(p) \ln(1-\bar{X}_\mathrm{Fe+}).
# \end{equation}
# In terms of the total amount of sulfur, $\bar{X}_\mathrm{Fe+}=1-\bar{X}_\mathrm{S}-\frac{1}{2}$, and substituting we have
# \begin{equation}
# \Delta T_\mathrm{m,FeS}(p)=T_\mathrm{m,FeS}(p) \ln\left[1-\left(\frac{1}{2}-\bar{X}_\mathrm{S}\right)\right]=
# T_\mathrm{m,FeS}(p) \ln\left(\frac{1}{2}+\bar{X}_\mathrm{S}\right).
# \end{equation}
# Ideally, these expressions would hold on their respective sides of the eutectic up to the eutectic itself, which would then be defined by their intersection, but it turns out that the experimentally determined eutectic runs substantially closer to the pure Fe side of the Fe-FeS join at most $p$.
#
# Only imposing the experimentally determined $\bar{X}_\mathrm{eut}(p)$ as the boundary between the two domains would violate the condition that the depressed solidi on both sides of the eutectic do not have discontinuities but must be identical at the eutectic itself. Hence we rescale the $\Delta T_\mathrm{m,Fe}$ and $\Delta T_\mathrm{m,FeS}$ using the experimentally constrained $\bar{X}_\mathrm{eut}(p)$ and $T_\mathrm{eut}(p)$ in order to enforce continuity and a correct position of the eutectic (Ruedas et al., 2013):
# \begin{equation}
# T_\mathrm{m}(p,\bar{X}_\mathrm{S})=
# \begin{cases}
# T_\mathrm{m,Fe}(p)-\frac{\ln(1-\bar{X}_\mathrm{S})}{\ln(1-\bar{X}_\mathrm{eut})}(T_\mathrm{m,Fe}(p)-T_\mathrm{eut}(p)),&\bar{X}_\mathrm{S}\leq\bar{X}_\mathrm{eut}\\
# T_\mathrm{m,FeS}(p)-\frac{\ln\left(\frac{1}{2}+\bar{X}_\mathrm{S}\right)}{\ln\left(\frac{1}{2}+\bar{X}_\mathrm{eut}\right)}(T_\mathrm{m,FeS}(p)-T_\mathrm{eut}(p)),&\text{else.}
# \end{cases}
# \end{equation}
# In[ ]:
def Tmalloy(p,X,mode):
from math import log
if X <= Xeut(p):
# Fe-rich side of eutectic
Tm=Tliq_Fe(p,mode)
Tm=Tm-(Tm-Teut(p))*log(1-X)/log(1-Xeut(p))
else:
# Fe-poor side of eutectic
Tm=Tliq_FeS(p)
Tm=Tm-(Tm-Teut(p))*log(0.5+X)/log(0.5+Xeut(p))
return Tm
# In practice, the fraction of the second component is usually given as a mass fraction $X_\mathrm{S}$ rather than as a mole fraction $\bar{X}_\mathrm{S}$. We convert from the former to the latter by
# \begin{equation}
# \bar{X}_\mathrm{S}=\left[1+\frac{u_\mathrm{S}}{u_\mathrm{Fe}}\left(\frac{1}{X_\mathrm{S}}-1\right)\right]^{-1},
# \end{equation}
# where $u_\mathrm{Fe}$ and $u_\mathrm{S}$ are the atomic masses of Fe and S, respectively. The atomic masses of elements here and elsewhere are taken from the IUPAC 2013 evaluation (Meija et al., 2016).
# In[ ]:
# mass/mole conversions for binary composites
def mass2mol(el,xel):
if el in u.keys():
return 1/(1.+u[el]/u['Fe']*(1./xel-1))
else:
raise NameError('Atomic mass for element '+el+' unavailable.')
def mol2mass(el,xel):
if el in u.keys():
return xel*u[el]/(xel*u[el]+(1-xel)*u['Fe'])
else:
raise NameError('Atomic mass for element '+el+' unavailable.')
# atomic masses (IUPAC 2013)
u={'O': 15.999e-3,'Na': 22.98976928e-3,'Mg': 24.305e-3,'S': 32.06e-3, 'K': 39.0983e-3,'Fe': 55.845e-3}
# ### <a id='main'></a>Main program (user interface)
# One can select the material (peridotite, basalt/eclogite, iron/iron alloy) and for peridotite also the system (terrestrial or martian peridotite, CMAS, or chondritic material, or else specify "other"). If the peridotite system "other" is selected, the mass fractions of MgO, FeO, Na<sub>2</sub>O, and K<sub>2</sub>O must be entered to perform the interpolation. If the alloys are selected, the mass fraction of S must be entered; if $X_\mathrm{S}=0$, only the data for pure Fe and no eutectic parameterization are shown.
# In[ ]:
import sys
from collections import defaultdict
import matplotlib.pyplot as plt
ox={}
p_crv=[]; Ts_crv=[]; Tl_crv=[]; Ti_crv=[]
ps_dat=defaultdict(list); pl_dat=defaultdict(list); Ts_dat=defaultdict(list); Tl_dat=defaultdict(list)
mat=int(input("Materials: peridotite (1), basalt/eclogite (2), iron/Fe-S alloy (3)\nSelect material: "))
plt.figure(figsize=(12,7))
if mat == 1:
# peridotite
# common liquidus phases: experimental data
# forsterite
pl_dat={'Davis & England (1964)': [0.55,0.55,1.25,1.25,1.8,2.5,3.05,3.95,3.95,4.65,4.65],
'Navrotsky et al. (1989), Richet et al. (1993)': [0.,0.],
'Ohtani & Kumazawa, 1981': [4.7,6.2,7.7,8.4,10.],
'Presnall & Walter, 1993': [9.7,10.6,11.8,12.8,13.9,13.9,14.9,15.6,16.5,16.7]}
Tl_dat={'Davis & England (1964)': [2178.,2203.,2223.,2248.,2253.,2303.,2303.,2353.,2378.,2378.,2403.],
'Navrotsky et al. (1989), Richet et al. (1993)': [2163.,2174.],
'Ohtani & Kumazawa, 1981': [2373.,2458.,2458.,2503.,2518.],
'Presnall & Walter, 1993': [2483.,2493.,2513.,2533.,2583.,2543.,2553.,2583.,2553.,2588.]}
mark_dat={'Davis & England (1964)': 'o','Navrotsky et al. (1989), Richet et al. (1993)': 'v', 'Ohtani & Kumazawa, 1981': '^','Presnall & Walter, 1993': '<'}
for nm in pl_dat.keys():
plt.plot(pl_dat[nm],Tl_dat[nm],linestyle='None',marker=mark_dat[nm],color='''0.2''')
pl_dat={}; Tl_dat={}; mark_dat={}
# pyrope
pl_dat={'Boyd & England (1962)': [3.6,3.95,3.95,4.3,4.3,4.65,4.65],
'Irifune & Ohtani (1986)': [3.5,3.5,5.,5.,5.5,5.5,6.25,6.25,7.,7.,7.5,7.5,8.,8.5,9.,9.5,10.],
'Kudo & Ito (1996)': [15.,15.],
'Ohtani et al. (1981)': [4.72,4.72,6.27,6.27,7.72,7.80,10.02],
'Shen & Lazor (1995)': [8.4,11.,11.,14.,22.,20.,25.],
'Zhang & Herzberg (1994)': [7.,8.,8.,9.,10.,10.,13.,16.,16.]}
Tl_dat={'Boyd & England (1962)': [2048.,2048.,2073.,2073.,2098.,2098.,2123.],
'Irifune & Ohtani (1986)': [2003.,2023.,2073.,2123.,2123.,2173.,2173.,2248.,2223.,2248.,2223.,2273.,2273.,2273.,2303.,2313.,2348.],
'Kudo & Ito (1996)': [2523.,2573.],
'Ohtani et al. (1981)': [2069.,2126.,2256.,2273.,2343.,2348.,2407.],
'Shen & Lazor (1995)': [2300.,2450.,2480.,2650.,2830.,2800.,3100.],
'Zhang & Herzberg (1994)': [2253.,2253.,2283.,2343.,2323.,2383.,2573.,2653.,2773.]}
mark_dat={'Boyd & England (1962)': 'x','Irifune & Ohtani (1986)': 's', 'Kudo & Ito (1996)': '*','Ohtani et al. (1981)': '+', 'Shen & Lazor (1995)': 'D','Zhang & Herzberg (1994)': 'h'}
for nm in pl_dat.keys():
plt.plot(pl_dat[nm],Tl_dat[nm],linestyle='None',marker=mark_dat[nm],color='''0.5''')
pl_dat={}; Tl_dat={}; mark_dat={}
# periclase
pl_dat={'Alfè (2005)': [0.4,17.,47.,135.6]}
Tl_dat={'Alfè (2005)': [3070.,4590.,6047.,8144.]}
mark_dat={'Alfè (2005)': '>'}
for nm in pl_dat.keys():
plt.plot(pl_dat[nm],Tl_dat[nm],linestyle='None',marker=mark_dat[nm],color='''0.8''')
pl_dat={}; Tl_dat={}; mark_dat={}
# common liquidus phases: parameterizations
p_crv=np.linspace(0.,14.5,30)
Tl_crv=[Tliqfo(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,color='''0.2''',linestyle="-.",label="Forsterite liquidus")
p_crv=np.linspace(2.,26.,49)
Tl_crv=[Tliqpy(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,color='''0.5''',linestyle="-.",label="Pyrope liquidus")
p_crv=np.linspace(0.,140.,271)
Tl_crv=[Tliqpc(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,color='''0.8''',linestyle="-.",label="Periclase liquidus")
# lherzolite composition selection
compsys=int(input("Compositions: terrestrial (1), martian (2), CMAS (3), chondritic (4), other (5):\n" "Select composition: "))
if compsys == 1:
title="Terrestrial peridotite"
# experimental data and parameterization: batch liquidus
pl_dat={'Fiquet et al. (2010)': [36.1,52.4,61.3,80.1,95.1,137.9],
'Herzberg et al. (1990)': [14.,14.],
'Ito & Takahashi (1987), upper mantle': [16.,18.,18.,20.,20.],
'Ito & Takahashi (1987), lower mantle': [25.,25.],
'Ito & Kennedy (1967)': [2.04],
'Scarfe & Takahashi (1986)': [2.,2.,3.5,3.5,5.,5.,6.,7.,9.,9.,11.,13.],
'Takahashi (1986)': [0.0001,1.5,1.5,3.,3.,5.,7.5,8.,8.,10.5,12.5,12.5,14.,14.],
'Takahashi & Scarfe (1985)': [1.5,1.5,3.,5.,5.,7.5,11.,14.],
'Takahashi et al. (1993)': [4.6,6.5],
'Trønnes & Frost (2002)': [24.]}
Tl_dat={'Fiquet et al. (2010)': [3714.,4101.,4300.,4453.,4801.,5401.],
'Herzberg et al. (1990)': [2388.,2428.],
'Ito & Takahashi (1987), upper mantle': [2273.,2381.,2473.,2482.,2682.],
'Ito & Takahashi (1987), lower mantle': [2773.,3273.],
'Ito & Kennedy (1967)': [2123.],
'Scarfe & Takahashi (1986)': [1980.,2076.,2075.,2175.,2239.,2271.,2228.,2281.,2221.,2236.,2227.,2194.],
'Takahashi (1986)': [1873.,2073.,2123.,2173.,2223.,2273.,2273.,2283.,2393.,2273.,2283.,2423.,2273.,2323.],
'Takahashi & Scarfe (1985)': [2073.,2181.,2200.,2184.,2278.,2278.,2280.,2325.],
'Takahashi et al. (1993)': [2173.,2273.],
'Trønnes & Frost (2002)': [2613.]}
mark_dat={'Fiquet et al. (2010)': '<','Herzberg et al. (1990)': 'o', 'Ito & Takahashi (1987), upper mantle': 'v','Ito & Takahashi (1987), lower mantle': 'v', 'Ito & Kennedy (1967)': '*','Scarfe & Takahashi (1986)': '^', 'Takahashi (1986)': '+','Takahashi & Scarfe (1985)': 'x', 'Takahashi et al. (1993)': 's','Trønnes & Frost (2002)': 'H'}
for nm in pl_dat.keys():
plt.plot(pl_dat[nm],Tl_dat[nm],linestyle='None',marker=mark_dat[nm],color='red')
p_crv=np.linspace(0.,170.,341)
Tl_crv=[Tliqb_per(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,"r--",label="Batch liquidus")
# experimental data (lower mantle) and parameterization: solidus and fractional liquidus
ps_dat={'Fiquet et al. (2010)': [36.2,51.8,61.3,80.2,95.1,106.6,122.0],
'Hirose & Fei (2002)': [25.,25.,27.],
'Ito & Takahashi (1987)': [25.],
'Trønnes & Frost (2002)': [23.,24.5,23.,24.5,24.,24.5],
'Zerr et al. (1997,1998)': [27.81,31.05,32.31,25.04,29.30,41.21,57.70,26.80,39.72,47.25,58.93],
'Zhang & Herzberg (1994)': [22.5]}
Ts_dat={'Fiquet et al. (2010)': [2817.,3124.,3196.,3690.,3930.,3908.,4100.],
'Hirose & Fei (2002)': [2773.,2473.,2723.],
'Ito & Takahashi (1987)': [2773.],
'Trønnes & Frost (2002)': [2413.,2403.,2433.,2503.,2453.,2503.],
'Zerr et al. (1997,1998)': [2924.,2994.,3128.,2272.,2495.,2832.,3255.,2880.,3139.,3404.,3456.],
'Zhang & Herzberg (1994)': [2408.]}
mark_dat={'Fiquet et al. (2010)': 'o','Hirose & Fei (2002)': 's', 'Ito & Takahashi (1987)': 'v','Trønnes & Frost (2002)': '^', 'Zerr et al. (1997,1998)': 'x','Zhang & Herzberg (1994)': '+'}
Ts_crv=[Tsol_Earth(p) for p in p_crv]
Tl_crv=[Tliq_per(p) for p in p_crv]
pTrange=[0.,175.,1350.,8000.]
elif compsys == 2:
title="Martian peridotite"
# experimental data and parameterization: solidus and fractional liquidus
ps_dat={'Bertka & Holloway (1994)': [0.5,1.,1.,1.5,1.5,2.,2.,2.2,3.],
'Collinet et al. (2015)': [1.,1.5,2.,2.2,0.5,0.5,1.,1.,1.5,1.5,2.,2.,2.2,2.2],
'Matsukage et al. (2013)': [3.,4.5],
'Schmerr et al. (2001)': [6.,6.,15.,15.,20.,23.,23.,25.,25.]}
Ts_dat={'Bertka & Holloway (1994)': [1368.,1443.,1463.,1523.,1543.,1583.,1603.,1593.,1693.],
'Collinet et al. (2015)': [1473.,1543.,1623.,1623.,1387.5,1387.5,1453.4,1453.4,1528.4,1528.4,1604.,1604.,1584.9,1584.9],
'Matsukage et al. (2013)': [1639.,1743.],
'Schmerr et al. (2001)': [1893.,1943.,2243.,2293.,2343.,2373.,2423.,2473.,2573.]}
mark_dat={'Bertka & Holloway (1994)': 'o','Collinet et al. (2015)': 'x', 'Matsukage et al. (2013)': '^','Schmerr et al. (2001)': 'v'}
p_crv=np.linspace(0.,25.,101)
Tl_crv=[Tliq_per(p) for p in p_crv]
Ts_crv=[Tsol_Mars(p) for p in p_crv]
pTrange=[0.,25.,1350.,3000.]
elif compsys == 3:
title="CMAS composition"
# experimental data and parameterization: solidus and fractional liquidus
ps_dat={'Asahara et al. (1998)': [6.5],
'Gudfinnsson & Presnall (1996)': [2.4,2.6,2.8,3.0,3.2,3.3,3.3,3.4],
'Herzberg et al. (1990)': [9.2,14.,14.,15.3,15.3]}
Ts_dat={'Asahara et al. (1998)': [2073.],
'Gudfinnsson & Presnall (1996)': [1768.,1798.,1823.,1848.,1868.,1878.,1883.,1888.],
'Herzberg et al. (1990)': [2243.,2328.,2393.,2333.,2458.]}
mark_dat={'Asahara et al. (1998)': 'o','Gudfinnsson & Presnall (1996)': 'v', 'Herzberg et al. (1990)': '^'}
p_crv=np.linspace(0.,25.,101)
Tl_crv=[Tliq_per(p) for p in p_crv]
Ts_crv=[Tsol_CMAS(p) for p in p_crv]
pTrange=[0.,25.,1400.,3000.]
elif compsys == 4:
title="Chondritic composition"
# experimental data and parameterization: batch liquidus
pl_dat={'Agee & Draper (2004)': [4.7,5.],
'Andrault et al. (2011)': [20.7,27.9,35.5,42.5,45.1,50.2,59.5,65.0,68.1,69.9,84.9,110.6,141.1],
'Berthet et al. (2009)': [1.,1.],
'Ohtani (1987)': [12.,15.,17.5,20.],
'Takahashi (1983)': [1.5]}
Tl_dat={'Agee & Draper (2004)': [2123.,2123.],
'Andrault et al. (2011)': [2525.,2724.,2947.,3117.,3021.,3153.,3308.,3499.,3499.,3643.,3946.,4397.,4777.],
'Berthet et al. (2009)': [1823.,1873.],
'Ohtani (1987)': [2283.,2303.,2363.,2423.],
'Takahashi (1983)': [1923.]}
mark_dat={'Agee & Draper (2004)': '+','Andrault et al. (2011)': 'o',
'Berthet et al. (2009)': 's','Ohtani (1987)': '^','Takahashi (1983)': 'v'}
for nm in pl_dat.keys():
plt.plot(pl_dat[nm],Tl_dat[nm],linestyle='None',marker=mark_dat[nm],color='red')
p_crv=np.linspace(0.,150.,301)
Tl_crv=[Tliqb_chon(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,"r--",label="Batch liquidus")
# experimental data and parameterization: solidus and fractional liquidus
ps_dat={'Andrault et al. (2011)': [33.5,38.9,47.4,56.9,62.1,62.8,65.0,74.0,82.1,98.1,101.9,108.3,108.6,125.1,139.3,140.2],
'Andrault et al. (2018)': [1.,11.6,4.7,9.8,17.4,1.4,7.,9.,16.,24.,13.5,2.2,28.,5.,5.,5.,5.,5.,20.],
'Berthet et al. (2009)': [1.],
'Ohtani (1987)': [12.,15.,17.5],
'Takahashi (1983)': [0.8,1.5,2.25,0.6,0.6,1.5,1.5,2.5,2.5,3.]}
Ts_dat={'Andrault et al. (2011)': [2608.,2656.,2697.,2850.,3060.,3102.,3002.,3248.,3342.,3555.,3705.,3800.,3647.,3997.,4147.,4202.],
'Andrault et al. (2018)': [1519.,2030.,1823.,2006.,2155.,1600.,1839.,1917.,2107.,2268.,2062.,1616.,2351.,1865.,1840.,1840.,1815.,1815.,2185.],
'Berthet et al. (2009)': [1473.],
'Ohtani (1987)': [2073.,2173.,2198.],
'Takahashi (1983)': [1473.,1473.,1573.,1423.,1523.,1473.,1573.,1573.,1673.,1673.]}
mark_dat={'Andrault et al. (2011)': 'o','Andrault et al. (2018)': 'x',
'Berthet et al. (2009)': 's','Ohtani (1987)': '^','Takahashi (1983)': 'v'}
Tl_crv=[Tliq_per(p) for p in p_crv]
Ts_crv=[Tsol_chon(p) for p in p_crv]
pTrange=[0.,150.,1400.,7000.]
else:
print("Oxide mass fractions (in %): ")
sys.stdout.flush()
for nm in ["MgO","FeO","Na2O","K2O"]:
# convert % to fraction
print("%s (Earth: %.2f%%; Mars: %.2f%%)" % (nm,100*ox_E[nm],100*ox_M[nm]))
sys.stdout.flush()
ox.update({nm: 1e-2*float(input(" Enter target "+nm+": "))})
# parameterizations: interpolated solidus and fractional liquidus
p_crv=np.linspace(0.,25.,101)
Tl_crv=[Tliq_per(p) for p in p_crv]
Ts_crv=[Tsol_intp(p,ox) for p in p_crv]
pTrange=[0.,25.,1350.,3000.]
title="Peridotitic composition with Mg#=%.3f" % Mg0
# finish plotting of points and curves
plt.plot(p_crv,Tl_crv,"r-",label="Fractional liquidus")
if compsys <= 4:
for nm in ps_dat.keys():
plt.plot(ps_dat[nm],Ts_dat[nm],linestyle='None',marker=mark_dat[nm],color='blue')
else:
# parameterizations: terrestrial and martian or CMAS solidi as reference
Ti_crv=[Tsol_Earth(p) for p in p_crv]
plt.plot(p_crv,Ti_crv,"b--",label="terr. Solidus")
if Mg0 < Mg0_E:
Ti_crv=[Tsol_Mars(p) for p in p_crv]
plt.plot(p_crv,Ti_crv,"b-.",label="mart. Solidus")
else:
Ti_crv=[Tsol_CMAS(p) for p in p_crv]
plt.plot(p_crv,Ti_crv,"b-.",label="CMAS Solidus")
plt.plot(p_crv,Ts_crv,"b-",label="Solidus")
elif mat == 2:
# basalt/eclogite
title="Basalt/eclogite"
# common liquidus phases: experimental data
# forsterite
pl_dat={'Davis & England (1964)': [0.55,0.55,1.25,1.25,1.8,2.5,3.05,3.95,3.95,4.65,4.65],
'Navrotsky et al. (1989), Richet et al. (1993)': [0.,0.],
'Ohtani & Kumazawa, 1981': [4.7,6.2,7.7,8.4,10.],
'Presnall & Walter, 1993': [9.7,10.6,11.8,12.8,13.9,13.9,14.9,15.6,16.5,16.7]}
Tl_dat={'Davis & England (1964)': [2178.,2203.,2223.,2248.,2253.,2303.,2303.,2353.,2378.,2378.,2403.],
'Navrotsky et al. (1989), Richet et al. (1993)': [2163.,2174.],
'Ohtani & Kumazawa, 1981': [2373.,2458.,2458.,2503.,2518.],
'Presnall & Walter, 1993': [2483.,2493.,2513.,2533.,2583.,2543.,2553.,2583.,2553.,2588.]}
mark_dat={'Davis & England (1964)': 'o','Navrotsky et al. (1989), Richet et al. (1993)': 'v', 'Ohtani & Kumazawa, 1981': '^','Presnall & Walter, 1993': '<'}
for nm in pl_dat.keys():
plt.plot(pl_dat[nm],Tl_dat[nm],linestyle='None',marker=mark_dat[nm],color='''0.2''')
pl_dat={}; Tl_dat={}; mark_dat={}
# diopside
pl_dat={'Boyd & England (1963)': [0.54,0.52,1.08,1.08,1.79,1.79,2.51,2.51,2.51,2.51,2.51,3.23,3.23,3.23,3.23,3.95,3.95,4.67,4.67],
'Gasparik (1996)': [7.,8.,9.,10.,11.,12.,13.,14.,15.],
'Richet & Bottinga (1984)': [0.],
'Shen & Lazor (1995, di)': [10.,15.6],
'Yoder (1952)': [0.,0.0285,0.0428,0.0429,0.0498,0.0498,0.0579,0.0581,0.0781,0.0781,0.0812,0.1003,0.1004,0.1014,0.1014,0.1014,0.1051,0.1107,0.1109,0.111,0.1116,0.1196,0.1371,0.1472,0.1475,0.1512,0.1514,0.1514,0.1797,0.1798,0.1978,0.1986,0.2039,0.2039,0.2041,0.2585,0.2587,0.259,0.3248,0.325,0.325,0.3722,0.421,0.4212,0.4212,0.4622,0.463,0.4987,0.4992,0.501]}
Tl_dat={'Boyd & England (1963)': [1723.,1748.,1798.,1823.,1873.,1898.,1923.,1948.,1973.,1948.,1973.,1998.,2023.,1998.,2023.,2048.,2073.,2098.,2123.],
'Gasparik (1996)': [2213.,2263.,2303.,2343.,2363.,2393.,2423.,2423.,2423.],
'Richet & Bottinga (1984)': [1665.],
'Shen & Lazor (1995, di)': [2390.,2500.],
'Yoder (1952)': [1664.5,1669.,1670.2,1670.2,1671.,1671.,1672.4,1671.5,1674.3,1674.7,1674.9,1676.5,1675.9,1677.6,1679.4,1679.,1679.7,1676.9,1678.2,1676.9,1677.4,1680.7,1681.3,1682.3,1685.2,1683.9,1683.9,1684.3,1687.6,1687.6,1692.4,1691.9,1692.5,1691.3,1691.3,1697.5,1697.9,1699.9,1706.9,1708.5,1707.3,1712.2,1718.,1718.4,1718.8,1722.5,1724.9,1729.4,1728.6,1730.3]}
mark_dat={'Boyd & England (1963)': 'x','Gasparik (1996)': '+', 'Richet & Bottinga (1984)': 's','Shen & Lazor (1995, di)': '*','Yoder (1952)': 'D'}
for nm in pl_dat.keys():
plt.plot(pl_dat[nm],Tl_dat[nm],linestyle='None',marker=mark_dat[nm],color='''0.4''')
pl_dat={}; Tl_dat={}; mark_dat={}
# pyrope
pl_dat={'Boyd & England (1962)': [3.6,3.95,3.95,4.3,4.3,4.65,4.65],
'Irifune & Ohtani (1986)': [3.5,3.5,5.,5.,5.5,5.5,6.25,6.25,7.,7.,7.5,7.5,8.,8.5,9.,9.5,10.],
'Kudo & Ito (1996)': [15.,15.],
'Ohtani et al. (1981)': [4.72,4.72,6.27,6.27,7.72,7.80,10.02],
'Shen & Lazor (1995, py)': [8.4,11.,11.,14.,22.,20.,25.],
'Zhang & Herzberg (1994)': [7.,8.,8.,9.,10.,10.,13.,16.,16.]}
Tl_dat={'Boyd & England (1962)': [2048.,2048.,2073.,2073.,2098.,2098.,2123.],
'Irifune & Ohtani (1986)': [2003.,2023.,2073.,2123.,2123.,2173.,2173.,2248.,2223.,2248.,2223.,2273.,2273.,2273.,2303.,2313.,2348.],
'Kudo & Ito (1996)': [2523.,2573.],
'Ohtani et al. (1981)': [2069.,2126.,2256.,2273.,2343.,2348.,2407.],
'Shen & Lazor (1995, py)': [2300.,2450.,2480.,2650.,2830.,2800.,3100.],
'Zhang & Herzberg (1994)': [2253.,2253.,2283.,2343.,2323.,2383.,2573.,2653.,2773.]}
mark_dat={'Boyd & England (1962)': 'h','Irifune & Ohtani (1986)': '1', 'Kudo & Ito (1996)': '2','Ohtani et al. (1981)': '>', 'Shen & Lazor (1995, py)': '*','Zhang & Herzberg (1994)': '3'}
for nm in pl_dat.keys():
plt.plot(pl_dat[nm],Tl_dat[nm],linestyle='None',marker=mark_dat[nm],color='''0.6''')
pl_dat={}; Tl_dat={}; mark_dat={}
# Ca-perovskite
pl_dat={'Gasparik et al. (1994)': [15.2],
'Zerr et al. (1997)': [16.,17.7,20.,20.8,23.1,23.9,27.6,30.5,34.4,36.1,38.3,43.]}
Tl_dat={'Gasparik et al. (1994)': [2783.],
'Zerr et al. (1997)': [2860.,2920.,3005.,3035.,3190.,3185.,3450.,3590.,3905.,4030.,4005.,4120.]}
mark_dat={'Gasparik et al. (1994)': '4','Zerr et al. (1997)': 'p'}
for nm in pl_dat.keys():
plt.plot(pl_dat[nm],Tl_dat[nm],linestyle='None',marker=mark_dat[nm],color='''0.8''')
pl_dat={}; Tl_dat={}; mark_dat={}
# common liquidus phases: parameterizations
p_crv=np.linspace(0.,14.5,30)
Tl_crv=[Tliqfo(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,color='''0.2''',linestyle="-.",label="Forsterite liquidus")
p_crv=np.linspace(0.,16.,33)
Tl_crv=[Tliqdi(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,color='''0.4''',linestyle="-.",label="Diopside liquidus")
p_crv=np.linspace(2.,26.,49)
Tl_crv=[Tliqpy(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,color='''0.6''',linestyle="-.",label="Pyrope liquidus")
p_crv=np.linspace(14.5,45.,62)
Tl_crv=[TliqCapv(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,color='''0.8''',linestyle="-.",label="Ca-perovskite liquidus")
# experimental data and parameterizations: batch and fractional liquidi
pl_dat={'Hirose & Fei (2002)': [22.],
'Yasuda et al. (1994)': [3.,3.,5.,5.,6.,7.5,7.5,10.,10.,12.,14.,18.,20.,20.]}
Tl_dat={'Hirose & Fei (2002)': [2523.],
'Yasuda et al. (1994)': [1673.,1723.,1873.,1898.,1933.,2073.,2098.,2173.,2273.,2273.,2373.,2373.,2373.,2473.]}
mark_dat={'Hirose & Fei (2002)': 'o','Yasuda et al. (1994)': 'x'}
for nm in pl_dat.keys():
plt.plot(pl_dat[nm],Tl_dat[nm],linestyle='None',marker=mark_dat[nm],color='red')
p_crv=np.linspace(0.,150.,71)
Tl_crv=[Tliqb_bas(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,"r--",label="Batch liquidus")
Tl_crv=[Tliq_bas(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,"r-",label="Fractional liquidus")
# experimental data and parameterization: solidus
ps_dat={'Andrault et al. (2014)': [59.5,59.8,69.2,72.5,80.0,80.4,83.0,86.4,87.0,102.9,107.0,111.1,113.4,118.0,119.6,128.0,138.0,56.0,56.1,59.5,60.2,60.5,80.0,80.4,87.0,102.9,107.0,113.4,118.0,119.6,120.2,123.3,123.9,128.0,138.0],
'Hirose et al. (1999)': [16.6,20.7,26.0,31.9,31.9,32.0,36.0,36.2,39.8,44.8,46.3,46.3,46.3,54.7,55.4,61.2,61.2,64.2,22.4,27.2],
'Hirose & Fei (2002)': [22.,26.,27.,27.5],
'Pradhan et al. (2015)': [48.2,58.8,79.4,89.5,102.9,114.2,128.5],
'Yasuda et al. (1994)': [3.,3.,3.5,4.,5.,5.,7.5,7.5,10.,10.,12.,14.,14.,16.,18.,20.]}
Ts_dat={'Andrault et al. (2014)': [3000.,2950.,3260.,3150.,3400.,3350.,3330.,3350.,3370.,3550.,3520.,3700.,3540.,3670.,3650.,3750.,3730.,3070.,2930.,3050.,2980.,3100.,3550.,3450.,3550.,3700.,3650.,3750.,3750.,3900.,3850.,3800.,3850.,3800.,3900.],
'Hirose et al. (1999)': [2356.,2547.,2553.,2646.,2768.,2721.,2768.,2817.,2818.,2914.,2901.,3011.,3137.,3160.,3174.,3390.,3413.,3204.,2431.,2626.],
'Hirose & Fei (2002)': [2423.,2673.,2703.,2703.],
'Pradhan et al. (2015)': [2940.,3170.,3420.,3550.,3670.,3770.,3927.],
'Yasuda et al. (1994)': [1523.,1573.,1698.,1573.,1773.,1823.,1923.,1973.,2023.,2073.,2123.,2223.,2273.,2273.,2373.,2373.]}
mark_dat={'Andrault et al. (2014)': '3','Hirose et al. (1999)': '+', 'Hirose & Fei (2002)': 'o','Pradhan et al. (2015)': 's', 'Yasuda et al. (1994)': 'x'}
for nm in ps_dat.keys():
plt.plot(ps_dat[nm],Ts_dat[nm],linestyle='None',marker=mark_dat[nm],color='blue')
Ts_crv=[Tsol_bas(p) for p in p_crv]
plt.plot(p_crv,Ts_crv,"b-",label="Solidus")
pTrange=[0.,150.,1200.,4000.]
elif mat == 3:
# iron/iron alloys
title="Fe"
# experimental data: pure Fe
pl_dat_Fe={'Strong et al. (1973), delta': [1.8,2.3,2.5,2.6,3.1,3.7,3.8,4.4,5.2,5.3],
'Liu & Bassett (1975), delta': [2.86],
'Swartzendruber (1982)': [0.],
'Strong et al. (1973), fcc': [5.2,5.3,5.6,5.7,5.8],
'Liu & Bassett (1975), fcc': [7.20,10.86,11.87,15.38,15.89,18.83,19.74],
'Williams et al. (1987)': [12.,23.,38.,45.,55.,60.,64.,66.,68.,71.,82.,92.,10.,22.,68.],
'Ringwood & Hibberson (1990)': [16.],
'Boehler (1993), fcc': [6.,13.,15.,18.,19.,27.,34.,36.,39.,47.,48.,52.,64.,73.,81.],
'Jephcoat & Besedin (1996)': [47.],
'Shen et al. (1998), fcc': [12.0,23.2,23.2,29.2,30.3,39.2,42.2,46.3,46.3],
'Ahrens et al. (2002)': [67.,70.,76.,79.],
'Shen et al. (2004)': [27.,27.,42.,42.,50.,50.,57.,57.],
'Boehler et al. (2008), fcc': [35.,50.,67.,42.,50.,54.,68.],
'Buono & Walker (2011)': [6.],
'Jackson et al. (2013)': [20.,28.,37.,40.,49.,50.,61.,82.],
'Anzellini et al. (2013), fcc': [51.5,52.7,52.8,52.8,52.8,53.2,54.1,64.0,64.6,64.7,64.7,64.8,64.9,65.1,65.5,74.8,75.3,75.5,75.7,75.7,75.7,76.2,76.2,76.3,76.4,76.4,76.8,76.8,76.9,77.1,77.2,77.6,77.8,78.7,84.4,84.8,84.8,85.1,85.1,86.0,86.1,86.2,86.4,86.8,87.9],
'Aquilanti et al. (2015), fcc': [74.9,75.5,79.9,80.4],
'Zhang et al. (2016)': [19.0,19.1,55.9,57.1,59.8],
'Boehler (1993), hcp': [93.,95.,106.,111.,111.,115.,117.,121.,121.,125.,126.,132.,134.,137.,139.,139.,143.,146.,146.,149.,149.,154.,154.,158.,158.,160.,167.,178.,185.,195.,195.,196.,197.],
'Shen et al. (1998), hcp': [75.3,80.2],
'Ma et al. (2004)': [60.,105.],
'Nguyen & Holmes (2004)': [225.],
'Boehler et al. (2008), hcp': [90.,116.,68.,85.,116.,130.],
'Tateno et al. (2010)': [134.,135.5,210.,210.,330.,348.,356.,371.,376.],
'Murphy et al. (2011)': [106.,121.,133.,151.],
'Anzellini et al. (2013), hcp': [99.3,99.9,99.9,115.0,116.0,116.0,117.0,117.0,118.0,128.0,128.0,129.0,133.0,134.0,134.0,135.0,136.0,136.0,136.0,158.0,158.0,159.0,182.0,204.0],
'Aquilanti et al. (2015), hcp': [93.0,93.5,116.4,116.8],
'Boehler & Ross (2015)': [88.,93.,99.,101.]}
Tl_dat_Fe={'Strong et al. (1973), delta': [1859.,1878.,1878.,1880.,1904.,1922.,1925.,1949.,1978.,1983.],
'Liu & Bassett (1975), delta': [1912.3],
'Swartzendruber (1982)': [1811.],
'Strong et al. (1973), fcc': [1978.,1983.,1993.,2004.,2004.],
'Liu & Bassett (1975), fcc': [2073.9,2161.1,2188.0,2247.5,2263.2,2326.1,2350.3],
'Williams et al. (1987)': [2271.,2341.,2576.,2600.,2882.,2694.,3412.,3106.,3059.,3294.,3576.,3600.,2388.,2600.,3459.],
'Ringwood & Hibberson (1990)': [2218.],
'Boehler (1993), fcc': [2024.,2165.,2259.,2243.,2306.,2400.,2447.,2525.,2525.,2573.,2573.,2541.,2682.,2729.,2745.],
'Jephcoat & Besedin (1996)': [2750.],
'Shen et al. (1998), fcc': [2080.,2230.,2380.,2300.,2440.,2490.,2520.,2530.,2720.],
'Ahrens et al. (2002)': [2698.,2776.,2776.,2729.],
'Shen et al. (2004)': [2243.,2604.,2416.,2698.,2557.,2667.,2588.,2887.],
'Boehler et al. (2008), fcc': [2381.,2543.,2736.,2597.,2713.,2674.,2844.],
'Buono & Walker (2011)': [2078.],
'Jackson et al. (2013)': [2210.,2280.,2500.,2400.,2560.,2650.,2840.,3025.],
'Anzellini et al. (2013), fcc': [3031.,3137.,2994.,3026.,3205.,2938.,3139.,3225.,3277.,3084.,3202.,3314.,3044.,3342.,3312.,3534.,3337.,3323.,3352.,3400.,3529.,3532.,3566.,3460.,3507.,3504.,3560.,3400.,3368.,3319.,3473.,3425.,3425.,3376.,3635.,3771.,3664.,3549.,3763.,3458.,3598.,3535.,3618.,3597.,3650.],
'Aquilanti et al. (2015), fcc': [2735.,2840.,2705.,2810.],
'Zhang et al. (2016)': [2116.,2123.,2835.,2858.,2792.],
'Boehler (1993), hcp': [2824.,2839.,2918.,2933.,3027.,3043.,3043.,2949.,3075.,3012.,3217.,3184.,3106.,3059.,3153.,3216.,3216.,3294.,3357.,3294.,3357.,3310.,3420.,3310.,3357.,3482.,3373.,3655.,3796.,3733.,3796.,3796.,3859.],
'Shen et al. (1998), hcp': [2790.,2790.],
'Ma et al. (2004)': [2750.,3510.],
'Nguyen & Holmes (2004)': [5100.],
'Boehler et al. (2008), hcp': [2952.,3045.,2844.,2975.,3176.,3168.],
'Tateno et al. (2010)': [3244.,3083.,3853.,4109.,5135.,5327.,5519.,5423.,5679.],
'Murphy et al. (2011)': [3520.,3880.,3930.,4160.],
'Anzellini et al. (2013), hcp': [3521.,3615.,3842.,3846.,3892.,4069.,3791.,3841.,3911.,4064.,4119.,3961.,4097.,4024.,4269.,4095.,4114.,4175.,4087.,4269.,4463.,4612.,4669.,4923.],
'Aquilanti et al. (2015), hcp': [2755.,2835.,3010.,3090.],
'Boehler & Ross (2015)': [2792.,2871.,2933.,3043.]}
mark_dat_Fe={'Strong et al. (1973), delta': 'v','Liu & Bassett (1975), delta': 'd','Swartzendruber (1982)': 'p', 'Strong et al. (1973), fcc': 'v','Liu & Bassett (1975), fcc': 'd', 'Williams et al. (1987)': 'x','Ringwood & Hibberson (1990)': 'p', 'Boehler (1993), fcc': '2','Jephcoat & Besedin (1996)': 'H','Shen et al. (1998), fcc': '<', 'Ahrens et al. (2002)': '^','Shen et al. (2004)': '<','Boehler et al. (2008), fcc': '>', 'Buono & Walker (2011)': 'p','Jackson et al. (2013)': 's', 'Anzellini et al. (2013), fcc': '.','Aquilanti et al. (2015), fcc': '+', 'Zhang et al. (2016)': 'H', 'Shen et al. (1998), hcp': '<','Boehler (1993), hcp': '2', 'Ma et al. (2004)': '*','Nguyen & Holmes (2004)': 'p', 'Boehler et al. (2008), hcp': '>','Tateno et al. (2010)': 'D', 'Murphy et al. (2011)': 'o','Anzellini et al. (2013), hcp': '.', 'Boehler & Ross (2015)': 'h','Aquilanti et al. (2015), hcp': '+'}
for nm in pl_dat_Fe.keys():
plt.plot(pl_dat_Fe[nm],Tl_dat_Fe[nm],linestyle='None',marker=mark_dat_Fe[nm],color='red')
# parameterization
p_crv=np.linspace(0.,400.,801)
Tl_crv=[Tliq_Fe(p,'f') for p in p_crv]
plt.plot(p_crv,Tl_crv,'r-',label="Fe, flat gradient")
Tl_crv=[Tliq_Fe(p,'s') for p in p_crv]
plt.plot(p_crv,Tl_crv,'r--',label="Fe, steep gradient")
# decide whether to consider FeS and other alloys as well
xS=0.
wS=0.01*float(input("S content (wt.%): "))
if wS > 0.:
xS=mass2mol('S',wS)
if xS > 0.5:
print("WARNING: Alloy contains more S than Fe: xS=%.2f mol%%" % (xS*100))
xS=0.5
wS=mol2mass('S',xS)
print(" Calculating curves for pure FeS: wS=%.2f wt.%%" % (wS*100))
title+=", FeS, and eutectic"
else:
title+=", FeS, eutectic, and Fe + %.2f wt.%% (%.2f mol%%) S" % (wS*100,xS*100)
# melting curve experimental data: pure FeS
pl_dat_FeS={'STP': [0.],
'Anderson & Ahrens (1996)': [136.,330.,360.],
'Boehler (1992)': [6.4,10.8,11.7,22.1,31.1,43.0],
'Ono et al. (2008)': [43.,53.,64.,77.,95.,121.,154.,175.,223.,253.]}
Tl_dat_FeS={'STP': [1467.],
'Anderson & Ahrens (1996)': [3240.,4210.,4310.],
'Boehler (1992)': [1702.,1787.,1802.,2067.,2182.,2348.],
'Ono et al. (2008)': [2406.,2483.,2594.,2710.,2839.,3004.,3191.,3325.,3697.,3953.]}
mark_dat_FeS={'STP': 's','Anderson & Ahrens (1996)': 'x','Boehler (1992)': 'o',
'Ono et al. (2008)': '+'}
for nm in pl_dat_FeS.keys():
plt.plot(pl_dat_FeS[nm],Tl_dat_FeS[nm],linestyle='None',marker=mark_dat_FeS[nm],color='blue')
# parameterization
Tl_crv=[Tliq_FeS(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,'b-',label="FeS")
# eutectic temperature experimental data
pl_dat_eut={'Boehler (1996)': [5.],
'Brett & Bell (1969)': [3.],
'Buono & Walker (2011)': [6.],
'Campbell et al. (2007)': [29.6,31.4,36.6,52.1,54.,59.5,60.6,79.4,39.,59.7],
'Chen et al. (2008)': [10.,14.],
'Chudinovskikh & Boehler (2007)': [32.2,41.8,36.4],
'Fei et al. (1997, 2000)': [0.,7.,14.,18.,21.],
'Kamada et al. (2010)': [80.],
'Li et al. (2001)': [25.],
'Morard et al. (2007)': [6.1,8.85,10,10.6,13.21,13.35,13.6,15.3,16.5,16.97],
'Morard et al. (2008)': [29.33,31.87,32.17,43.,43.09,46.29,47.38,53.03,52.63,52.61,62.14,62.93,64.85],
'Mori et al. (2017)': [34.,46.,46.,60.,80.,130.,134.,200.,254.],
'Ryzhenko & Kennedy (1973)': [0.0001,0.0001,1.44,1.6,1.6,1.6,2.,2.,2.,2.,2.,2.,2.4,4.,4.,6.,6.,1.8,1.8,1.92,1.92,2.,2.,2.32,2.32,2.4,2.75,2.8,3.2,4.,4.,4.,4.,4.,6.,6.,6.,6.,3.,3.08,3.2,3.68,3.68,3.68,3.68,3.72,3.8,3.96,4.04,4.12,4.3,4.56,4.6,5.,5.2,5.36,5.48,5.68,5.76,5.8,5.96,6.06,2.,3.12,3.12,3.2,3.2,4.,4.,4.,4.,4.,4.,4.,4.,4.,4.,4.04,5.,5.12,5.2,5.8,5.96,6.,6.,6.,6.,6.,6.,6.],
'Stewart et al. (2007)': [23.,40.],
'Usselman (1975)': [3.,4.5,5.5]}
Tl_dat_eut={'Boehler (1996)': [1318.],
'Brett & Bell (1969)': [1263.],
'Buono & Walker (2011)': [1263.],
'Campbell et al. (2007)': [1520.,1510.,1520.,1710.,1780.,1930.,1870.,2040.,1820.,2000.],
'Chen et al. (2008)': [1191.,1132.],
'Chudinovskikh & Boehler (2007)': [1500.,1596.,1547.],
'Fei et al. (1997, 2000)': [1261.,1197.,1133.,1133.,1348.],
'Kamada et al. (2010)': [2197.],
'Li et al. (2001)': [1424.],
'Morard et al. (2007)': [1140.,1160.,1200.,1100.,1180.,1250.,1200.,1100.,1270.,1120.],
'Morard et al. (2008)': [1560.,1422.,1479.,1660.,1732.,1730.,1795.,1720.,1813.,1864.,1911.,1966.,1997.],
'Mori et al. (2017)': [1630.,1770.,1900.,1910.,2050.,2910.,2960.,3320.,3550.],
'Ryzhenko & Kennedy (1973)': [1260.,1262.,1260.,1283.,1264.,1273.,1262.,1266.,1279.,1285.,1284.,1282.,1276.,1292.,1292.,1283.,1283.,1287.,1291.,1288.,1292.,1275.,1292.,1293.,1293.,1293.,1291.,1295.,1288.,1287.,1287.,1284.,1284.,1278.,1265.,1269.,1269.,1274.,1288.,1289.,1292.,1290.,1290.,1290.,1290.,1285.,1285.,1291.,1291.,1293.,1289.,1287.,1279.,1276.,1275.,1272.,1281.,1286.,1277.,1273.,1274.,1274.,1284.,1292.,1292.,1301.,1292.,1296.,1302.,1302.,1299.,1297.,1304.,1291.,1292.,1289.,1286.,1299.,1294.,1281.,1281.,1283.,1274.,1283.,1285.,1283.,1284.,1285.,1286.,1286.],
'Stewart et al. (2007)': [1312.,1508.],
'Usselman (1975)': [1265.,1267.,1278.]}
mark_dat_eut={'Boehler (1996)': 'p','Brett & Bell (1969)': 'p','Buono & Walker (2011)': 'p', 'Campbell et al. (2007)': 'd','Chen et al. (2008)': 'h', 'Chudinovskikh & Boehler (2007)': 'o','Fei et al. (1997, 2000)': '+',
'Kamada et al. (2010)': 'p','Li et al. (2001)': 'p',\
'Morard et al. (2007)': 'x','Morard et al. (2008)': 'x',\
'Mori et al. (2017)': 's','Ryzhenko & Kennedy (1973)': '2',
'Stewart et al. (2007)': 'D','Usselman (1975)': 'v'}
for nm in pl_dat_eut.keys():
plt.plot(pl_dat_eut[nm],Tl_dat_eut[nm],linestyle='None',marker=mark_dat_eut[nm],color='cyan')
# eutectic temperature parameterization
Tl_crv=[Teut(p) for p in p_crv]
plt.plot(p_crv,Tl_crv,'c-',label="Eutectic")
if xS < 0.5:
Tl_crv=[Tmalloy(p,xS,'f') for p in p_crv]
plt.plot(p_crv,Tl_crv,'k-',label="Alloy, flat Fe liq.")
Tl_crv=[Tmalloy(p,xS,'s') for p in p_crv]
plt.plot(p_crv,Tl_crv,'k--',label="Alloy, steep Fe liq.")
pTrange=[0.,400.,1000.,6000.]
else:
pTrange=[0.,400.,1500.,6000.]
else:
raise IOError("Unknown material.")
# finish p-T plot
plt.axis(pTrange)
plt.xlabel('$p$ (GPa)')
plt.ylabel('$T_\mathrm{s,l}$ (K)')
plt.legend()
plt.title(title)
plt.show()
# also show information about the eutectic composition for alloys