pythonFLAC/P1figS1.py at main · chingchenn/pythonFLAC · GitHub

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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 26 11:09:22 2023

@author: chingchen
"""

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Helvetica"
fontsize=25
labelsize=23
bwith = 5
fig6,(ax,ax2,ax3) = plt.subplots(1,3,figsize=(18,10))
axdep1 = ax.twinx()
depth_limit = 200


# ---------------------------------- chlorite ------------------------------
pres = np.linspace(0,7,100) # GPa
TTT=np.zeros(len(pres))
for q,dd in enumerate(pres):
    ss1 = 700+150*dd
    ss2 = 764.545+53.7*dd
    ss3 = 973.447-43.478*dd
    ss4 = -63.64*dd+1040.016
    ss5 = -200*dd+1640
    TTT[q] = min(ss1,ss2,ss3,ss4,ss5)
lab1=ax.plot(TTT,pres,c='#D14309',lw=5,label='chlorite-peridotite')


# #### KATZ 2003, dry solidus
# x=np.array([0,1,2.6,4,6,7])
# y=np.array([1180,1200,1460,1560,1660,1750])
# pres = np.linspace(0,2000,1000)

# a,b,c = np.polyfit(y, x, 2)
# yyyy = a*pres**2+b*pres+c
# lab3=ax.plot(pres[pres>1000],yyyy[pres>1000],c='green',lw=5,label='dry solidus')

##-------------------------------- serpentinite ----------------------------
phase_change = (140/255, 20/255, 70/255)
x = np.linspace(500,730)
tt1 = 2.1 +(7.5-2.1)* (x - 730)/ (500-730)
lab4=ax.plot(x,tt1,c=phase_change,label = 'serpentinite-perdotite',lw=5)
x = np.linspace(620,730)
ttold = 2.1 + (0.2-2.1) * (x-730)/(650-730)
ax.plot(x,ttold,c=phase_change,lw=5)

depth = np.linspace(0,300000,100)
presss = np.linspace(0,9.9,100)
sss=np.zeros(len(depth))
for q,dd in enumerate(presss):
    ss1 = 980+0.6e-3*(dd*3.3e4-140e3)
    ss2 = 1090-178*(1-np.exp(-dd*4.125))
    sss[q] = max(ss1,ss2)
lab5=ax.plot(sss,presss,c='#FF9900',lw=5,label='solidus')
# lns = lab3+lab4+lab5+lab1#+lab2
# labs = [l.get_label() for l in lns]
#ax.legend(lns, labs, fontsize = fontsize-7, loc='lower right',bbox_to_anchor=(0.95, 0.15))

axdep2 = ax2.twinx()
ax2.vlines(x=650, ymin=0, ymax=3, colors='#6B8E23', ls='-', lw=5,)
ax2.hlines(y=3, xmin=200, xmax=650, colors='#6B8E23', ls='-', lw=5,)

depth = np.linspace(0,300000,100)
sss=np.zeros(len(depth))
for q,dd in enumerate(depth):
    ss1 = 680+0.6e-3*(dd-140e3)
    ss2 = 930-313*(1-np.exp(-dd/7e3))
    sss[q] = max(ss1,ss2)
lab3=axdep2.plot(sss,depth/1e3,c='#FF9900',lw=5,label='solidus')

axdep3 = ax3.twinx()
tempr = np.linspace(0,513)
y = -0.0375 * tempr + 20.1
basalt_change = (50/255, 200/255, 180/255)
ax3.plot(tempr,y,c=basalt_change,lw=5)

x = np.linspace(515,1000)
y = 0.0022 * x - 0.3
lab3=ax3.plot(x,y,c=basalt_change,lw=5,label='basalt-eclogite')
## Dash line for eclogite
x2 = np.linspace(1000,1200)
y2 = 0.0022 * x2 - 0.3
ax3.plot(x2,y2,c=basalt_change,lw=5,linestyle=(5, (3, 1)))

pressure_limit = 2.45 # GPa
pressure=np.linspace(0,pressure_limit,100)
sss=np.zeros(len(pressure))
for q,dd in enumerate(pressure):
    if dd<1:
        ss=1050-420*(1-np.exp(-dd*3.3))
    elif dd>2.38:
        ss=(dd+14)*43
    else:
        ss=630+26*((-dd)**2)/2
    sss[q] = ss
lab4=ax3.plot(sss,pressure,c='#FF9900',lw=5,label='solidus')

lns = lab3+lab4
labs = [l.get_label() for l in lns]
#ax3.legend(lns, labs, fontsize = fontsize-7, loc='lower right',bbox_to_anchor=(0.9, 0.4))

#ax2.text(770,1.5,'schist',fontsize=36)
#ax3.text(570,3.5,'Eclogite',fontsize=26)
ax.set_ylabel('Pressure (GPa)',fontsize=fontsize)
axdep3.set_ylabel('Depth (km)',fontsize=fontsize)
xmajor_ticks=np.array([300,600,900,1200])
ymajor_ticks=np.array([40,80,120,160,200])
for aa in [ax,ax2,ax3]:
    for axis in ['top','bottom','left','right']:
        aa.spines[axis].set_linewidth(bwith)
    aa.tick_params(labelsize=labelsize,width=3,length=10,right=False, top=True,direction='in')
    aa.set_xlim(300,1200)
    aa.grid()
    aa.set_ylim(depth_limit*1000*10*3300/1e9,0)
    aa.set_xlabel('Temperature ($^\circ$C)',fontsize=fontsize)
    aa.set_xticks(xmajor_ticks)
for axdep in [axdep1,axdep2,axdep3]:
    axdep.set_ylim(depth_limit,0)
    axdep.tick_params(axis='y',labelsize=labelsize,width=3,length=10,right=True, top=True,direction='in')
    axdep.set_ylim(depth_limit,0)
    axdep.set_yticks(ymajor_ticks)
    #axdep.grid()
fig6.savefig('/Users/chingchen/Desktop/FLAC_Works/Eclogite_flat_slab/phase_diagram_v2.pdf')