Skip to content
Snippets Groups Projects
Commit 0d952a64 authored by Miro Pütz's avatar Miro Pütz
Browse files

Replace strength_envelope.py

parent 012fe9bf
Branches master
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 81 additions and 107 deletions
strength_envelope.py
+ 81
− 107
View file @ 0d952a64
......@@ -3,102 +3,72 @@
# strength_envelope.py - A Python script for plotting strength envelopes in the
# continental lithosphere
#
# dwhipp 10.15
# dwhipp 10.15
# modified by mpuetz
#--- Import libraries ---
import numpy as np
import matplotlib.pyplot as plt
temp_array = np.load('temp_array.npy') #load the calculated temp_array
#--- User-defined variables ---------------------------------------------------#
# Thermal model - 1D conduction with heat production (flux basal BC)
qmax=40. # Basal heat flow [mW m-2]
S=0.15 # Heat production [microW m-3]
k=2.75 # Thermal conductivity [W m-1 K-1]
#--- plastic deformation - Byerlee's law ---
rhoc=2750. # Crust density [kg m-3]
rhom=3300. # Mantle density [kg m-3]
mu1=0.85 # Friction coefficient below 500 MPa
mu2=0.6 # Friction coefficient above 500 MPa
coh1=0. # Cohesion below 500 MPa
coh2=50. # Cohesion above 500 MPa
lam=0.4 # Pore fluid to lithostatic stress ratio
zc=35. # Crustal thickness [km]
zl=125. # Lithospheric thickness [km]
#--- Ductile deformation - Dorn's law ---
edot=1.e-14
# Wet quartzite
Aq=5.e-6 # Power law constant [MPa-n s-1]
Qq=1.9e5 # Activation energy [J mol-1]
nq=3. # Power law exponent
# Olivine
Ao=7.e4 # Power law constant [Pa-n s-1]
Qo=5.2e5 # Activation energy [J mol-1]
no=3. # Power law exponent
# Dorn's law
QD=5.454e5 # Activation energy for olivine creep [J mol-1]
edotD=5.7e11 # Strain rate [s-1]
sigmaD=8500. # Critical stress [MPa]
# Temperature
#Ts=0 # Surface temperature [deg. C]
#Tl=1300 # Temperature of LAB [deg. C]
#--- End user-defined variables -----------------------------------------------#
# Constants
R=8.314 # Universal gas constant [J mol-1 K-1]
sigmaB=500. # Pressure cutoff for Byerlee's law
# Conversions
zc=zc*1000. # km -> m
zl=zl*1000. # km -> m
coh1=coh1*1.e6 # MPa -> Pa
coh2=coh2*1.e6 # MPa -> Pa
Aq=(Aq/1.e6**nq) # MPa-n s-1 -> Pa-n s-1
Ao=(Ao/1.e6**no) # MPa-n s-1 -> Pa-n s-1
#print(Ao)
sigmaB=sigmaB*1.e6 # MPa -> Pa
sigmaD=sigmaD*1.e6 # MPa -> Pa
qmax=qmax/1000. # mW m-2 -> W m-2
S=S/1.e6 # microW m-3 -> W m-3
# Ranges
z = np.linspace(0, 150000, 201) # the program needs that information right now
#z=np.linspace(0.,zl,1000)
#T=pylab.linspace(Ts,Tl,1000)
#strength(temperature, depth vector, depth of the crust)
def strength(T, z, ZC):
#--- plastic deformation - Byerlee's law ---
rhoc=2750. # Crust density [kg m-3]
rhom=3300. # Mantle density [kg m-3]
mu1=0.85 # Friction coefficient below 500 MPa
mu2=0.6 # Friction coefficient above 500 MPa
coh1=0. # Cohesion below 500 MPa
coh2=50. # Cohesion above 500 MPa
lam=0.4 # Pore fluid to lithostatic stress ratio
#--- Ductile deformation - Dorn's law ---
edot=1.e-14
# Wet quartzite
Aq=5.e-6 # Power law constant [MPa-n s-1]
Qq=1.9e5 # Activation energy [J mol-1]
nq=3. # Power law exponent
# Olivine
Ao=7.e4 # Power law constant [Pa-n s-1]
Qo=5.2e5 # Activation energy [J mol-1]
no=3. # Power law exponent
# Dorn's law
QD=5.454e5 # Activation energy for olivine creep [J mol-1]
edotD=5.7e11 # Strain rate [s-1]
sigmaD=8500. # Critical stress [MPa]
#--- End user-defined variables -----------------------------------------------#
# Constants
R=8.314 # Universal gas constant [J mol-1 K-1]
sigmaB=500. # Pressure cutoff for Byerlee's law
# Conversions
coh1=coh1*1.e6 # MPa -> Pa
coh2=coh2*1.e6 # MPa -> Pa
Aq=(Aq/1.e6**nq) # MPa-n s-1 -> Pa-n s-1
Ao=(Ao/1.e6**no) # MPa-n s-1 -> Pa-n s-1
sigmaB=sigmaB*1.e6 # MPa -> Pa
sigmaD=sigmaD*1.e6 # MPa -> Pa
def strength(T, z):
plastic_c=np.zeros(np.size(z))
plastic_t=np.zeros(np.size(z))
visc=np.zeros(np.size(z))
viscD=np.zeros(np.size(z))
viscplot=np.zeros(np.size(z))
strength_c=np.zeros(np.size(z))
strength_t=np.zeros(np.size(z))
pressure=np.zeros(np.size(z))
# Calculate temperature NOT necessary anymore
#T = -(S * z**2.)/(2.*k) + (qmax + S * zl)/k * z
#print(-(S * zc**2.)/(2.*k) + (qmax + S * zl)/k * zc)
T = T + 273.15
for i in range(0,np.size(z)):
if z[i] < zc:
if z[i] < ZC:
pressure[i]=rhoc*9.81*z[i]
plastic_c[i]=(2.*(coh1 + mu1*pressure[i]*(1.-lam)))/(np.sqrt(mu1**2. + 1.) - mu1)
plastic_t[i]=(-2.*(coh1 - mu1*pressure[i]*(1.-lam)))/(np.sqrt(mu1**2. + 1.) + mu1)
if pressure[i] > 2.e8:
plastic_c[i]=(2.*(coh2 + mu2*pressure[i]*(1.-lam)))/(np.sqrt(mu2**2. + 1.) - mu2)
#plastic_t[i]=(-2.*(-coh2 - mu2*pressure[i]*(1.-lam)))/(np.sqrt(mu2**2. + 1.) + mu2)
plastic_t[i]=(-2.*(-coh2 - mu2*pressure[i]*(1.-lam)))/(np.sqrt(mu2**2. + 1.) + mu2)
visc[i]=(edot/Aq)**((1./nq))*np.exp(Qq/(nq*R*T[i]))
else:
pressure[i]=(rhoc*9.81*zc)+(rhom*9.81*(z[i]-zc))
pressure[i]=(rhoc*9.81*ZC)+(rhom*9.81*(z[i]-ZC))
plastic_c[i]=(2.*(coh2 + mu2*pressure[i]*(1.-lam)))/(np.sqrt(mu2**2. + 1.) - mu2)
plastic_t[i]=(-2.*(coh2 - mu2*pressure[i]*(1.-lam)))/(np.sqrt(mu2**2. + 1.) + mu2)
visc[i]=(edot/Ao)**((1./no))*np.exp(Qo/(no*R*T[i]))
......@@ -115,42 +85,46 @@ def strength(T, z):
else:
strength_t[i]=visc[i]
print((edot/Aq)**((1./nq))*np.exp(Qq/(nq*R*773.15))*1.e-6)
print((edot/Ao)**((1./no))*np.exp(Qo/(no*R*773.15))*1.e-6)
print(sigmaD*(1.-np.sqrt((R*(773.15)/QD*np.log(edotD/edot))))*1.e-6)
total_strength = np.sum(strength_t)
return total_strength
#print((edot/Aq)**((1./nq))*np.exp(Qq/(nq*R*773.15))*1.e-6)
#print((edot/Ao)**((1./no))*np.exp(Qo/(no*R*773.15))*1.e-6)
#print(sigmaD*(1.-np.sqrt((R*(773.15)/QD*np.log(edotD/edot))))*1.e-6)
# Scale for plotting
T = T - 273.15
strength_c = strength_c * 1.e-6
strength_t = strength_t * 1.e-6
plastic_c = plastic_c * 1.e-6
plastic_t = plastic_t * 1.e-6
visc = visc * 1.e-6
z=z/1000.
#T = T - 273.15
#strength_c = strength_c * 1.e-6
#strength_t = strength_t * 1.e-6
#plastic_c = plastic_c * 1.e-6
#plastic_t = plastic_t * 1.e-6
#visc = visc * 1.e-6
#z=z/1000.
# plt.figure()
# plt.subplot(121)
# plt.plot(T,z,'r')
# plt.xlabel('Temperature [${\circ}$C]')
# plt.ylabel('Depth [km]')
# plt.plot([min(T),max(T)],[zc/1000.,zc/1000.],'k--')
# plt.grid()
# plt.axis([0., 1.05*max(T), min(z), max(z)])
# plt.gca().invert_yaxis()
# plt.subplot(122)
# plt.plot(strength_c,z,'k',label='Compression')
# plt.plot(strength_t,z,'b',label='Tension')
# plt.gca().axes.yaxis.set_ticklabels([])
# plt.xlabel('Differential stress [MPa]')
# plt.grid()
# plt.legend(loc=4)
# plt.axis([0., 1.05*max(strength_c), min(z), max(z)])
plt.figure()
plt.subplot(121)
plt.plot(T,z,'r')
plt.xlabel('Temperature [${\circ}$C]')
plt.ylabel('Depth [km]')
plt.plot([min(T),max(T)],[zc/1000.,zc/1000.],'k--')
plt.grid()
plt.axis([0., 1.05*max(T), min(z), max(z)])
plt.gca().invert_yaxis()
# plt.gca().invert_yaxis()
# plt.show()
plt.subplot(122)
plt.plot(strength_c,z,'k',label='Compression')
plt.plot(strength_t,z,'b',label='Tension')
plt.gca().axes.yaxis.set_ticklabels([])
plt.xlabel('Differential stress [MPa]')
plt.grid()
plt.legend(loc=4)
plt.axis([0., 1.05*max(strength_c), min(z), max(z)])
plt.gca().invert_yaxis()
plt.show()
# --------------------------- strength calculation ----------------------------#
strength(temp_array[:,0], z)
strength(temp_array[:,100], z)
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment