Skip to content
Snippets Groups Projects
Commit dd20916f authored by Aleksi Rantanen's avatar Aleksi Rantanen
Browse files

New MELTS compositions

Materials are not yet defined to have these new compositions
parent 0cc5a963
No related branches found
No related tags found
No related merge requests found
......@@ -49,56 +49,154 @@ l3 = Layer(Z=150, k=2.75, H=0, RHO=3300, CP=1300)
# ----------------------------- Melts Data ------------------------------------#
# Melts data: First column contain temperatures, and first and second rows contain pressures and depths respectively
# Temperatures decrease from up to down, and pressures and depths increase from left to right
# Rock data sheets
dryMafic = gft("dryMafic.csv", delimiter=",")
dryFelsic = gft("dryFelsic.csv", delimiter=",")
wet2Mafic = gft("wetMafic2.csv", delimiter=",")
wetImed = gft("wetImed3.csv", delimiter=",")
wetFelsic = gft("wetFelsic4.csv", delimiter=",")
# interpolating the MELTS data: pressure vs temperature (gives melt fraction)
dryMF = interp.RectBivariateSpline(dryMafic[2:,0], dryMafic[0,1:], dryMafic[2:,1:], kx=3, ky=3)
dryFF = interp.RectBivariateSpline(dryFelsic[2:,0], dryFelsic[0,1:], dryFelsic[2:,1:], kx=3, ky=3)
wet2MF = interp.RectBivariateSpline(wet2Mafic[2:,0], wet2Mafic[0,1:], wet2Mafic[2:,1:], kx=3, ky=3)
wetIF = interp.RectBivariateSpline(wetImed[2:,0], wetImed[0,1:], wetImed[2:,1:], kx=3, ky=3)
wetFF = interp.RectBivariateSpline(wetFelsic[2:,0], wetFelsic[0,1:], wetFelsic[2:,1:], kx=3, ky=3)
# Creating Pressure vs Melts fraction grid which gives the values for themperatures
Tn = np.linspace(0, 1, 1921) # Last number defines the number of melt fraction steps
Tn[-1] = 0.99999999 # Highest temperature value
interpDryM = [] # lists of melt fractions at interpolated pressure ranges
interpDryF = []
interpWet2M = []
interpWetI = []
interpWetF = []
for i in range(len(dryMafic[0,1:])):
interpDryM.append(np.interp(Tn, dryMafic[2:,i+1], dryMafic[2:,0])) # x(melt), y(T)
for i in range(len(dryFelsic[0,1:])):
interpDryF.append(np.interp(Tn, dryFelsic[2:,i+1], dryFelsic[2:,0])) # x(melt), y(T)
for i in range(len(wet2Mafic[0,1:])):
interpWet2M.append(np.interp(Tn, wet2Mafic[2:,i+1], wet2Mafic[2:,0])) # x(melt), y(T)
for i in range(len(wetImed[0,1:])):
interpWetI.append(np.interp(Tn, wetImed[2:,i+1], wetImed[2:,0])) # x(melt), y(T)
for i in range(len(wetFelsic[0,1:])):
interpWetF.append(np.interp(Tn, wetFelsic[2:,i+1], wetFelsic[2:,0])) # x(melt), y(T)
# Creating 2D array of Pressure (x axis) vs melt fraction (y axis)
interpDryM = np.transpose(np.array(interpDryM))
interpDryF = np.transpose(np.array(interpDryF))
interpWet2M = np.transpose(np.array(interpWet2M))
interpWetI = np.transpose(np.array(interpWetI))
interpWetF = np.transpose(np.array(interpWetF))
# Interpolating the melt fraction vs pressure data (gives temperature)
dryMT = interp.RectBivariateSpline(Tn, dryMafic[0,1:], interpDryM, kx=3, ky=3)
dryFT = interp.RectBivariateSpline(Tn, dryFelsic[0,1:], interpDryF, kx=3, ky=3)
wet2MT = interp.RectBivariateSpline(Tn, wet2Mafic[0,1:], interpWet2M, kx=3, ky=3)
wetIT = interp.RectBivariateSpline(Tn, wetImed[0,1:], interpWetI, kx=3, ky=3)
wetFT = interp.RectBivariateSpline(Tn, wetFelsic[0,1:], interpWetF, kx=3, ky=3)
def MELTS():
'''
- Reading in the melts data sheets
- Creating interpolated P vs T melt fraction arrays
- Creating P vs F temperature arrays
- Creating interpolated P vs F temperature arrays
'''
# Melts data: First column contain temperatures, and first and second rows contain pressures and depths respectively
# Temperatures decrease from up to down, and pressures and depths increase from left to right
# Rock data sheets (mantle, dry and wet ultramafic, dry and wet mafic, dry and wet intermadiate, and dry and wet felsic)
dryMantle = gft("dryMantle.csv", delimiter=",")
dryUltra = gft("dryUltra.csv", delimiter=",")
wet05Ultra = gft("wetUltra05.csv", delimiter=",")
dryMafic = gft("dryMafic.csv", delimiter=",")
wet05Mafic = gft("wetMafic05.csv", delimiter=",")
wet1Mafic = gft("wetMafic1.csv", delimiter=",")
wet2Mafic = gft("wetMafic2.csv", delimiter=",")
dryImed = gft("dryImed.csv", delimiter=",")
wet075Imed = gft("wetImed075.csv", delimiter=",")
wet15Imed = gft("wetImed15.csv", delimiter=",")
wet3Imed = gft("wetImed3.csv", delimiter=",")
dryFelsic = gft("dryFelsic.csv", delimiter=",")
wet1Felsic = gft("wetFelsic1.csv", delimiter=",")
wet2Felsic = gft("wetFelsic2.csv", delimiter=",")
wet4Felsic = gft("wetFelsic4.csv", delimiter=",")
# interpolating the MELTS data: pressure vs temperature (gives melt fraction)
# Mantle and ultramafic copmositions
dryMantleF = interp.RectBivariateSpline(dryMantle[2:,0], dryMantle[0,1:], dryMantle[2:,1:], kx=3, ky=3)
dryUF = interp.RectBivariateSpline(dryUltra[2:,0], dryUltra[0,1:], dryUltra[2:,1:], kx=3, ky=3)
wet05UF = interp.RectBivariateSpline(wet05Ultra[2:,0], wet05Ultra[0,1:], wet05Ultra[2:,1:], kx=3, ky=3)
# Mafic compositions
dryMF = interp.RectBivariateSpline(dryMafic[2:,0], dryMafic[0,1:], dryMafic[2:,1:], kx=3, ky=3)
wet05MF = interp.RectBivariateSpline(wet05Mafic[2:,0], wet05Mafic[0,1:], wet05Mafic[2:,1:], kx=3, ky=3)
wet1MF = interp.RectBivariateSpline(wet1Mafic[2:,0], wet1Mafic[0,1:], wet1Mafic[2:,1:], kx=3, ky=3)
wet2MF = interp.RectBivariateSpline(wet2Mafic[2:,0], wet2Mafic[0,1:], wet2Mafic[2:,1:], kx=3, ky=3)
# Intermediate compositions
dryIF = interp.RectBivariateSpline(dryImed[2:,0], dryImed[0,1:], dryImed[2:,1:], kx=3, ky=3)
wet075IF = interp.RectBivariateSpline(wet075Imed[2:,0], wet075Imed[0,1:], wet075Imed[2:,1:], kx=3, ky=3)
wet15IF = interp.RectBivariateSpline(wet15Imed[2:,0], wet15Imed[0,1:], wet15Imed[2:,1:], kx=3, ky=3)
wet3IF = interp.RectBivariateSpline(wet3Imed[2:,0], wet3Imed[0,1:], wet3Imed[2:,1:], kx=3, ky=3)
# Felsic compositions
dryFF = interp.RectBivariateSpline(dryFelsic[2:,0], dryFelsic[0,1:], dryFelsic[2:,1:], kx=3, ky=3)
wet1FF = interp.RectBivariateSpline(wet1Felsic[2:,0], wet1Felsic[0,1:], wet1Felsic[2:,1:], kx=3, ky=3)
wet2FF = interp.RectBivariateSpline(wet2Felsic[2:,0], wet2Felsic[0,1:], wet2Felsic[2:,1:], kx=3, ky=3)
wet4FF = interp.RectBivariateSpline(wet4Felsic[2:,0], wet4Felsic[0,1:], wet4Felsic[2:,1:], kx=3, ky=3)
# Creating Pressure vs Melts fraction grid which gives the values for themperatures
Tn = np.linspace(0, 1, 1921) # Last number defines the number of melt fraction steps
Tn[-1] = 0.99999999 # Highest temperature value
interpDryMantle = [] # lists of melt fractions at interpolated pressure ranges
interpDryU = []
interpWet05U = []
interpDryM = []
interpWet05M = []
interpWet1M = []
interpWet2M = []
interpDryI = []
interpWet075I = []
interpWet15I = []
interpWet3I = []
interpDryF = []
interpWet1F = []
interpWet2F = []
interpWet4F = []
for i in range(len(dryMantle[0,1:])): # Mantle
interpDryMantle.append(np.interp(Tn, dryMantle[2:,i+1], dryMantle[2:,0])) # x(melt), y(T)
for i in range(len(dryUltra[0,1:])): # dry ultramafic
interpDryU.append(np.interp(Tn, dryUltra[2:,i+1], dryUltra[2:,0]))
for i in range(len(wet05Ultra[0,1:])): # wet ultramafic
interpWet05U.append(np.interp(Tn, wet05Ultra[2:,i+1], wet05Ultra[2:,0]))
for i in range(len(dryMafic[0,1:])): # mafic compositions
interpDryM.append(np.interp(Tn, dryMafic[2:,i+1], dryMafic[2:,0]))
for i in range(len(wet05Mafic[0,1:])):
interpWet05M.append(np.interp(Tn, wet05Mafic[2:,i+1], wet05Mafic[2:,0]))
for i in range(len(wet1Mafic[0,1:])):
interpWet1M.append(np.interp(Tn, wet1Mafic[2:,i+1], wet1Mafic[2:,0]))
for i in range(len(wet2Mafic[0,1:])):
interpWet2M.append(np.interp(Tn, wet2Mafic[2:,i+1], wet2Mafic[2:,0]))
for i in range(len(dryImed[0,1:])): # intermediate compositions
interpDryI.append(np.interp(Tn, dryImed[2:,i+1], dryImed[2:,0]))
for i in range(len(wet075Imed[0,1:])):
interpWet075I.append(np.interp(Tn, wet075Imed[2:,i+1], wet075Imed[2:,0]))
for i in range(len(wet15Imed[0,1:])):
interpWet15I.append(np.interp(Tn, wet15Imed[2:,i+1], wet15Imed[2:,0]))
for i in range(len(wet3Imed[0,1:])):
interpWet3I.append(np.interp(Tn, wet3Imed[2:,i+1], wet3Imed[2:,0]))
for i in range(len(dryFelsic[0,1:])): # felsic compositions
interpDryF.append(np.interp(Tn, dryFelsic[2:,i+1], dryFelsic[2:,0]))
for i in range(len(wet1Felsic[0,1:])):
interpWet1F.append(np.interp(Tn, wet1Felsic[2:,i+1], wet1Felsic[2:,0]))
for i in range(len(wet2Felsic[0,1:])):
interpWet2F.append(np.interp(Tn, wet2Felsic[2:,i+1], wet2Felsic[2:,0]))
for i in range(len(wet4Felsic[0,1:])):
interpWet4F.append(np.interp(Tn, wet4Felsic[2:,i+1], wet4Felsic[2:,0]))
# Creating 2D array of Pressure (x axis) vs melt fraction (y axis)
interpDryMantle = np.transpose(np.array(interpDryMantle))
interpDryU = np.transpose(np.array(interpDryU))
interpWet05U = np.transpose(np.array(interpWet05U))
interpDryM = np.transpose(np.array(interpDryM))
interpWet05M = np.transpose(np.array(interpWet05M))
interpWet1M = np.transpose(np.array(interpWet1M))
interpWet2M = np.transpose(np.array(interpWet2M))
interpDryI = np.transpose(np.array(interpDryI))
interpWet075I = np.transpose(np.array(interpWet075I))
interpWet15I = np.transpose(np.array(interpWet15I))
interpWet3I = np.transpose(np.array(interpWet3I))
interpDryF = np.transpose(np.array(interpDryF))
interpWet1F = np.transpose(np.array(interpWet1F))
interpWet2F = np.transpose(np.array(interpWet2F))
interpWet4F = np.transpose(np.array(interpWet4F))
# Interpolating the melt fraction vs pressure data (gives temperature)
dryMantleT = interp.RectBivariateSpline(Tn, dryMantle[0,1:], interpDryMantle, kx=3, ky=3)
dryUT = interp.RectBivariateSpline(Tn, dryUltra[0,1:], interpDryU, kx=3, ky=3)
wet05UT = interp.RectBivariateSpline(Tn, wet05Ultra[0,1:], interpWet05U, kx=3, ky=3)
dryMT = interp.RectBivariateSpline(Tn, dryMafic[0,1:], interpDryM, kx=3, ky=3)
wet05MT = interp.RectBivariateSpline(Tn, wet05Mafic[0,1:], interpWet05M, kx=3, ky=3)
wet1MT = interp.RectBivariateSpline(Tn, wet1Mafic[0,1:], interpWet1M, kx=3, ky=3)
wet2MT = interp.RectBivariateSpline(Tn, wet2Mafic[0,1:], interpWet2M, kx=3, ky=3)
dryIT = interp.RectBivariateSpline(Tn, dryImed[0,1:], interpDryI, kx=3, ky=3)
wet075IT = interp.RectBivariateSpline(Tn, wet075Imed[0,1:], interpWet075I, kx=3, ky=3)
wet15IT = interp.RectBivariateSpline(Tn, wet15Imed[0,1:], interpWet15I, kx=3, ky=3)
wet3IT = interp.RectBivariateSpline(Tn, wet3Imed[0,1:], interpWet3I, kx=3, ky=3)
dryFT = interp.RectBivariateSpline(Tn, dryFelsic[0,1:], interpDryF, kx=3, ky=3)
wet1FT = interp.RectBivariateSpline(Tn, wet1Felsic[0,1:], interpWet1F, kx=3, ky=3)
wet2FT = interp.RectBivariateSpline(Tn, wet2Felsic[0,1:], interpWet2F, kx=3, ky=3)
wet4FT = interp.RectBivariateSpline(Tn, wet4Felsic[0,1:], interpWet4F, kx=3, ky=3)
return dryMantleF,dryMantleT,dryUF,dryUT,wet05UF,wet05UT,dryMF,dryMT,wet05MF,wet05MT,wet1MF,wet1MT,wet2MF,wet2MT,dryIF,dryIT,wet075IF,wet075IT,wet15IF,wet15IT,wet3IF,wet3IT,dryFF,dryFT,wet1FF,wet1FT,wet2FF,wet2FT,wet4FF,wet4FT
# interpolated MELTS as F(P,T) and T(P,F) arrays
dryMantleF,dryMantleT,dryUF,dryUT,wet05UF,wet05UT,dryMF,dryMT,wet05MF,wet05MT,wet1MF,wet1MT,wet2MF,wet2MT,dryIF,dryIT,wet075IF,wet075IT,wet15IF,wet15IT,wet3IF,wet3IT,dryFF,dryFT,wet1FF,wet1FT,wet2FF,wet2FT,wet4FF,wet4FT = MELTS()
def material_plots():
......@@ -108,8 +206,8 @@ def material_plots():
dryMaficCurves = dryMT(melt_fractions, pressures)
wet2MaficCurves = wet2MT(melt_fractions, pressures)
wetIntermediateCurves = wetIT(melt_fractions, pressures)
wetFelsicCurves = wetFT(melt_fractions, pressures)
wet3IntermediateCurves = wet3IT(melt_fractions, pressures)
wet4FelsicCurves = wet4FT(melt_fractions, pressures)
# Scale for plotting
pressures /= 1.e9
......@@ -144,10 +242,10 @@ def material_plots():
plt.subplot(223).set_title("Wet Intermediate", fontweight="bold", size=20)
plt.subplot(223).tick_params(labelsize=15)
plt.plot(wetIntermediateCurves[0,:], pressures, 'blue', label="Solidus")
plt.plot(wetIntermediateCurves[-1,:], pressures,'red', label="Liquidus")
plt.plot(wet3IntermediateCurves[0,:], pressures, 'blue', label="Solidus")
plt.plot(wet3IntermediateCurves[-1,:], pressures,'red', label="Liquidus")
for i in range(len(melt_fractions[1:-1])):
plt.plot(wetIntermediateCurves[i+1], pressures, color=(melt_fractions[i+1], 0.0, 1-melt_fractions[i], 0.3))
plt.plot(wet3IntermediateCurves[i+1], pressures, color=(melt_fractions[i+1], 0.0, 1-melt_fractions[i], 0.3))
plt.grid()
plt.gca().set_ylim([0, 1.5])
plt.gca().set_xlim([500, 1400])
......@@ -159,10 +257,10 @@ def material_plots():
plt.subplot(224).set_title("Wet Felsic", fontweight="bold", size=20)
plt.subplot(224).tick_params(labelsize=15)
plt.plot(wetFelsicCurves[0,:], pressures, 'blue', label="Solidus")
plt.plot(wetFelsicCurves[-1,:], pressures,'red', label="Liquidus")
plt.plot(wet4FelsicCurves[0,:], pressures, 'blue', label="Solidus")
plt.plot(wet4FelsicCurves[-1,:], pressures,'red', label="Liquidus")
for i in range(len(melt_fractions[1:-1])):
plt.plot(wetFelsicCurves[i+1], pressures, color=(melt_fractions[i+1], 0.0, 1-melt_fractions[i], 0.4))
plt.plot(wet4FelsicCurves[i+1], pressures, color=(melt_fractions[i+1], 0.0, 1-melt_fractions[i], 0.4))
plt.grid()
plt.gca().set_ylim([0, 1.5])
plt.gca().set_xlim([500, 1400])
......@@ -184,8 +282,8 @@ class Material(object):
Common class for all Layers.
Arguments:
k : Thermal conductivity [W/(m*K)]
H : Heat Production [microW/m**3] Values from: "D. Hasterock, J. Webb (2017). On the radiogenic heat production of igneous rocks"
RHO : Density [kg/m**3] Values from: "D. Hasterock, J. Webb (2017). On the radiogenic heat production of igneous rocks"
H : Heat Production [microW/m**3] D. Hasterock, M. Gard, and J. Webb (2018)
RHO : Density [kg/m**3] D. Hasterock, M. Gard, and J. Webb (2018)
CP : Heat capacity [J/(kg*K)]
K : Thermal diffusivity [m**2/s]
F : a melt fraction array of T vs P
......@@ -194,14 +292,13 @@ class Material(object):
Q : Activation energy [kJ/mol] Values from Burov (2015) and references within
n : Power law exponent Values from Burov (2015) and references within """
def __init__(self, k, H, RHO, CP, L, A, n, Q, scaling, melts, F, T, name, index):
def __init__(self, k, H, RHO, CP, L, A, n, Q, scaling, F, T, name, index):
self.k = k
self.H = H/1.e6
self.RHO = RHO
self.CP = CP
self.L = L
self.K = k/(RHO*CP)
self.melts = melts
self.F = F
self.T = T
self.n = n
......@@ -213,63 +310,81 @@ class Material(object):
"""
Materials For the Country Rock (upper crust, lower crust, and mantle):
Upper crust always uses deformation parameters for wet quartzite (even when the material is intermediary)
Lower Crust always uses deformation parameters for dry Maryland Diabase scaled by 0.1
Moho uses Predetermined deformation parameters found in strength_envelope script (A, Q, and n are dummy variables) """
Materials For the Lithosphere:
Upper crust always uses deformation parameters for wet quartzite
Lower Crust always uses deformation parameters for dry Maryland Diabase scaled by 0.1 or Pikwitonei granulite
Mantle uses Predetermined deformation parameters found in strength_envelope script (A, Q, and n here are dummy variables)
Deformation parameters:
A=4.5e5, n=4, Q=223 : Wet Quartzite (Gleason and Tullis (1995))
A=8.0, n=4.7, Q=485 : Maryland Diabase (Strong) (Mackwell et al. (1998))
A=1.4e4, n=4.7, Q=485 : Undried Pikwitonei granulite (Wilkin and Carter (1990))
A=2.0e-4, n=1.9, Q=140 : Wet Granite (Mackwell et al. (1998)) (wet felsic composition)
A=3.2e-2, n=2.4, Q=212 : Wet Diorite (Ranalli (1995)) (wet intermediate composition)
A=5.0e5, n=3.4, Q=497: Microgabbro (Wilks and Carter (1990))
A=5.0, n=3.5, Q=530: Dry Olivine (Hirth and Kohlstedt (2003))
Heat Production (D. Hasterock, M. Gard, and J. Webb (2018)):
RHO=2657, H=3.461 : granite (felsic)
RHO=2718, H=1.838 : granodiorite (felsic)
RHO=2792, H=1.221 : diorite (intermediate)
RHO=2905, H=0.725 : gabbroic diorite (mafic)
RHO=2990, H=0.498 : subalkalic gabbro (mafic)
RHO=3194, H=0.466 : peridotgabbro (ultramafic)
RHO=3299, H=0.692 : mantle peridotite (ultramafic)
RHO=2658, H=3.513 : syenite (felsic)
RHO=2683, H=2.872 : quartz monzonite (felsic)
RHO=2761, H=2.196 : monzonite (intermediate)
RHO=2838, H=1.521 : monzodiorite (intermediate)
RHO=2910, H=1.482 : monzogabbro (intemediate)
RHO=2958, H=0.743 : alkalic gabbro (mafic)
# Upper Crust:
# Dry felsic rock (material with dry felsic melts data and wet quartzite disclocation creep parameters)
mat1 = Material(k=2.75, H=1.5, RHO=2700, CP=1000, L=4.5e5, A=1.1e-4, n=4, Q=223, scaling=1, name="Dry Felsic Upper Crust", melts=dryFelsic, F=dryFF, T=dryFT, index=1) # A, Q, n: Wet Quartzite (Gleason and Tullis (1995))
# Wet felsic Rock (material with wet felsic melts data and wet quartzite disclocation creep parameters)
mat2 = Material(k=2.75, H=1.5, RHO=2700, CP=1000, L=4.5e5, A=1.1e-4, n=4, Q=223, scaling=1, name="Wet Felsic Upper Crust", melts=wetFelsic, F=wetFF, T=wetFT, index=2) # A, Q, n: Wet Quartzite (Gleason and Tullis (1995))
# Wet intermediary rocks (materials with wet intermadiary melts data and wet quartzite)
mat3 = Material(k=2.75, H=1.5, RHO=2700, CP=1000, L=4.0e5, A=1.1e-4, n=4, Q=223, scaling=1, name="Wet Intermediate Upper Crust", melts=wetImed, F=wetIF, T=wetIT, index=3) # A, Q, n: Wet Quartzite (Gleason and Tullis (1995))
"""
# Lower Crust:
# Dry mafic rock (materials with dry mafic melts data and Maryland Diabase)
mat4 = Material(k=3.4, H=0.2, RHO=2900, CP=1000, L=4.0e5, A=8.0, n=4.7, Q=485, scaling=0.1, name="Dry Mafic Lower Crust", melts=dryMafic, F=dryMF, T=dryMT, index=4) # A, Q, n: Maryland Diabase (Strong) (Mackwell et al. (1998))
# Dry mafic rock (materials with dry mafic melts data and Undried Adirondack granulite)
mat5 = Material(k=3.4, H=0.2, RHO=2900, CP=1000, L=4.0e5, A=1.4e4, n=4.2, Q=445, scaling=1, name="Dry Mafic Granulite Lower Crust", melts=dryMafic, F=dryMF, T=dryMT, index=5) # A, Q, n: Undried Pikwitonei granulite (Wilkin and Carter (1990))
# Felsic compositions
# Granite:
mat1 = Material(k=2.75, H=1.82, RHO=2650, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Dry Felsic Granodiorite", F=dryFF, T=dryFT, index=1)
mat2 = Material(k=3.1, H=3.54, RHO=2700, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Dry Felsic Granite", F=wet1FF, T=wet1FT, index=2)
mat3 = Material(k=2.75, H=1.82, RHO=2650, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Wet Felsic Granodiorite", F=wet4FF, T=wet4FT, index=3)
mat4 = Material(k=3.1, H=3.54, RHO=2700, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Wet Felsic Granite", F=wet4FF, T=wet4FT, index=4)
# Granodiorite:
mat1 = Material(k=2.75, H=1.82, RHO=2650, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Dry Felsic Granodiorite", F=dryFF, T=dryFT, index=1)
mat2 = Material(k=3.1, H=3.54, RHO=2700, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Dry Felsic Granite", F=wet1FF, T=wet1FT, index=2)
mat3 = Material(k=2.75, H=1.82, RHO=2650, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Wet Felsic Granodiorite", F=wet4FF, T=wet4FT, index=3)
mat4 = Material(k=3.1, H=3.54, RHO=2700, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Wet Felsic Granite", F=wet4FF, T=wet4FT, index=4)
# Intermediate compositions
mat5 = Material(k=3.2, H=1.20, RHO=2800, CP=1000, L=4.0e5, A=3.2e-2, n=2.4, Q=212, scaling=1, name="Wet Intermediate Diorite", F=wet3IF, T=wet3IT, index=5)
mat6 = Material(k=3.2, H=1.20, RHO=2800, CP=1000, L=4.0e5, A=3.2e-2, n=2.4, Q=212, scaling=1, name="Wet Intermediate Diorite", F=wet3IF, T=wet3IT, index=6)
mat7 = Material(k=3.2, H=1.20, RHO=2800, CP=1000, L=4.0e5, A=3.2e-2, n=2.4, Q=212, scaling=1, name="Wet Intermediate Diorite", F=wet3IF, T=wet3IT, index=7)
mat8 = Material(k=3.2, H=1.20, RHO=2800, CP=1000, L=4.0e5, A=3.2e-2, n=2.4, Q=212, scaling=1, name="Wet Intermediate Diorite", F=wet3IF, T=wet3IT, index=8)
# Mantle:
# Mafic compositions
mat9 = Material(k=3.0, H=0.67, RHO=2900, CP=1000, L=4.0e5, A=5.0e9, n=3.4, Q=497, scaling=1, name="Wet Mafic Gabbroic Diorite", F=wet2MF, T=wet2MT, index=9)
mat10 = Material(k=3.0, H=0.46, RHO=3000, CP=1000, L=4.0e5, A=5.0e9, n=3.4, Q=497, scaling=1, name="Wet Mafic Subalkalic Gabbro", F=wet2MF, T=wet2MT, index=10)
mat11 = Material(k=3.0, H=0.46, RHO=3000, CP=1000, L=4.0e5, A=5.0e9, n=3.4, Q=497, scaling=1, name="Dry Mafic Subalkalic Gabbro", F=dryMF, T=dryMT, index=11)
mat12 = Material(k=3.75, H=0.0, RHO=3200, CP=1300, L=3.5e5, A=5.0, n=3.5, Q=530, scaling=1, name="Dry Mafic Peridotite Gabbro", F=dryMF, T=dryMT, index=12)
# Dry mafic rock (materials with dry mafic melts data)
mat6 = Material(k=3.75, H=0.0, RHO=3300, CP=1300, L=3.5e5, A=5.0, n=3.5, Q=530, scaling=1, name="Mantle", melts=dryMafic, F=dryMF, T=dryMT, index=6) # A, Q, n:
# Ultramafic compositions
mat13 = Material(k=3.75, H=0.0, RHO=3200, CP=1300, L=3.5e5, A=5.0, n=3.5, Q=530, scaling=1, name="Dry Mafic Peridotite Gabbro", F=dryMF, T=dryMT, index=13)
mat14 = Material(k=3.75, H=0.0, RHO=3200, CP=1300, L=3.5e5, A=5.0, n=3.5, Q=530, scaling=1, name="Dry Mafic Peridotite Gabbro", F=dryMF, T=dryMT, index=14)
"""
Materials for the intrusions
Intrusions at MOHO must always be Dry Mafic in composition
Intrusions in the Crust can have any Rhyolite-MELTS data """
# Upper Crust:
# Dry felsic intrusion
mat7 = Material(k=2.75, H=1.82, RHO=2650, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Dry Felsic Granodiorite", melts=dryFelsic, F=dryFF, T=dryFT, index=7) # A, Q, n: Wet Granite (Mackwell et al. (1998))
mat8 = Material(k=3.1, H=3.54, RHO=2700, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Dry Felsic Granite", melts=dryFelsic, F=dryFF, T=dryFT, index=8) # A, Q, n: Wet Granite (Mackwell et al. (1998))
# Wet felsic intrusion
mat9 = Material(k=2.75, H=1.82, RHO=2650, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Wet Felsic Granodiorite", melts=wetFelsic, F=wetFF, T=wetFT, index=9) # A, Q, n: Wet Granite (Mackwell et al. (1998))
mat10 = Material(k=3.1, H=3.54, RHO=2700, CP=1100, L=4.0e5, A=2.0e-4, n=1.9, Q=140, scaling=1, name="Wet Felsic Granite", melts=wetFelsic, F=wetFF, T=wetFT, index=10) # A, Q, n: Wet Granite (Mackwell et al. (1998))
# Wet intermediate intrusion
mat11 = Material(k=3.2, H=1.20, RHO=2800, CP=1000, L=4.0e5, A=3.2e-2, n=2.4, Q=212, scaling=1, name="Wet Intermediate Diorite", melts=wetImed, F=wetIF, T=wetIT, index=11) # A, Q, n: Wet Diorite (Ranalli (1995))
# Wet mafic intrusion
mat12 = Material(k=3.0, H=0.67, RHO=2900, CP=1000, L=4.0e5, A=5.0e9, n=3.4, Q=497, scaling=1, name="Wet Mafic Gabbroic Diorite", melts=wet2Mafic, F=wet2MF, T=wet2MT, index=12) # A, Q, n: Microgabbro (Wilks and Carter (1990))
mat13 = Material(k=3.0, H=0.46, RHO=3000, CP=1000, L=4.0e5, A=5.0e9, n=3.4, Q=497, scaling=1, name="Wet Mafic Subalkalic Gabbro", melts=wet2Mafic, F=wet2MF, T=wet2MT, index=13) # A, Q, n: Microgabbro (Wilks and Carter (1990))
# Dry mafic intrusion
mat14 = Material(k=3.0, H=0.46, RHO=3000, CP=1000, L=4.0e5, A=5.0e9, n=3.4, Q=497, scaling=1, name="Dry Mafic Subalkalic Gabbro", melts=dryMafic, F=dryMF, T=dryMT, index=14) # A, Q, n: Microgabbro (Wilks and Carter (1990))
mat15 = Material(k=3.75, H=0.0, RHO=3200, CP=1300, L=3.5e5, A=5.0, n=3.5, Q=530, scaling=1, name="Dry Mafic Peridotite Gabbro", melts=dryMafic, F=dryMF, T=dryMT, index=15) # A, Q, n: Dry Olivine (Hirth and Kohlstedt (2003))
# Dry felsic rock (material with dry felsic melts data and wet quartzite disclocation creep parameters)
mat15 = Material(k=2.75, H=1.5, RHO=2700, CP=1000, L=4.5e5, A=1.1e-4, n=4, Q=223, scaling=1, name="Dry Felsic Upper Crust", F=dryFF, T=dryFT, index=15)
# Lower Crust:
mat16 = Material(k=3.4, H=0.2, RHO=2900, CP=1000, L=4.0e5, A=1.4e4, n=4.2, Q=445, scaling=1, name="Dry Mafic Granulite Lower Crust", F=dryMF, T=dryMT, index=16)
# Test material
mat16 = Material(k=2.75, H=0, RHO=2900, CP=1000, L=4.0e5, A=5.0e9, n=3.4, Q=497, scaling=1, name="Dry Mafic Granulite Lower Crust", melts=dryMafic, F=dryMF, T=dryMT, index=16)
# Mantle
mat17 = Material(k=3.75, H=0.0, RHO=3300, CP=1300, L=3.5e5, A=5.0, n=3.5, Q=530, scaling=1, name="Mantle", F=dryMantleF, T=dryMantleT, index=17)
# a material dictionary mat which takes in all different Material Class objects
# a material dictionary which takes in all different Material Class objects
mat = {1: mat1, 2: mat2, 3: mat3, 4: mat4, 5: mat5, 6: mat6, 7: mat7, 8:mat8, 9:mat9, 10:mat10, 11:mat11, 12:mat12, 13:mat13, 14:mat14, 15:mat15, 16:mat16}
......@@ -673,7 +788,7 @@ for i in range(TIMESTEPS):
mask[y,x,c.index] = True
# Choosing materials (7 - 15) for the intrusions (upper, middle, and lower)
if c.index % 3 == 0:
index_mat = 9 #11
index_mat = 9
matnum[y,x] = index_mat
H[y,x] = mat[index_mat].H
CP[y,x] = mat[index_mat].CP
......@@ -687,7 +802,7 @@ for i in range(TIMESTEPS):
# scaling[y,x] = mat[index_mat].scaling
elif c.index % 3 == 1:
index_mat = 13 #13
index_mat = 13
matnum[y,x] = index_mat
H[y,x] = mat[index_mat].H
CP[y,x] = mat[index_mat].CP
......@@ -701,7 +816,7 @@ for i in range(TIMESTEPS):
# scaling[y,x] = mat[index_mat].scaling
else:
index_mat = 15 #15
index_mat = 15
matnum[y,x] = index_mat
H[y,x] = mat[index_mat].H
CP[y,x] = mat[index_mat].CP
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment