Skip to content

Instantly share code, notes, and snippets.

@wassname
Last active August 29, 2015 14:02
Show Gist options
  • Save wassname/b5f55a26bf022048f730 to your computer and use it in GitHub Desktop.
Save wassname/b5f55a26bf022048f730 to your computer and use it in GitHub Desktop.
Model and graph heatflow resulting from an emplacement for varying time before present and depth below msl of emplacment
# -*- coding: utf-8 -*-
"""
Created on Sun Feb 16 20:03:59 2014
This will generate a schematic of heat vs depth vs time for difference emplacement times
@author: wassname
"""
import numpy as np
from numpy import pi, sin, cos, exp
from matplotlib import pyplot as plt
from pylab import *
from scipy.special import erf
import collections
import functools
class memoized(object):
'''Decorator. Caches a function's return value each time it is called.
If called later with the same arguments, the cached value is returned
(not reevaluated).
'''
def __init__(self, func):
self.func = func
self.cache = {}
def __call__(self, *args):
if not isinstance(args, collections.Hashable):
# uncacheable. a list, for instance.
# better to not cache than blow up.
return self.func(*args)
if args in self.cache:
return self.cache[args]
else:
value = self.func(*args)
self.cache[args] = value
return value
def __repr__(self):
'''Return the function's docstring.'''
return self.func.__doc__
def __get__(self, obj, objtype):
'''Support instance methods.'''
return functools.partial(self.__call__, obj)
#==============================================================================
# PARAMETERS
#==============================================================================
# T and Q data
"""
Funnell2006
standard basal heat flow south islan (60+-4 mW/m2), resulting in a predicted present day surface heat flow of 75 mW/m2
Duendin has 85+-10% mW/m2
@Godfrey2001
The change in mantle density (35 kg m -3) represented by the load in the flexural
model is related to a change in temperature of---354 ø, assuming a
thermal expansion coefficient of 3 x 10 -s K -1 [Turcotte and Schubert,1982, p. 182].
This is consistent with temperatures predicted in the upper mantle from the
high heat flow measured over Dunedin (---860 degC at 30 km depth) compared with
the region outside Dunedin (---490 degC at 30 km depth) [Funnell and Allis, 1996; Cook et al., 1999]
"""
#==============================================================================
# # PHYSICAL PARAMS
#==============================================================================
# generate each point on a 2d cross section x vs z, then for each point work out the distance and therefore heat at a fixed time
Tm=490+273.15 # ambient tempreture in K
T0=833+273.15 # increased tempreture K # a change in tempt of 354 degC is predicted at 30km depth (Funell and Allis 1996;, Cook et al 1999)
lambd=2.4 #thermal cond W/m/K .2.4 W/m/k value used in Godfrey et al 2001 for granites and gabbros with no ref but is constisten with:
#1.7-4 W/m/k from granite thermal cond http://en.wikipedia.org/wiki/List_of_thermal_conductivities
pc=4.0E6 #volume heat capacity J/m^3/K, used in Godfrey et al 2001 with no ref but is constisten with:
# volumetric heat capacity of granite 2.17E6 J/cm^3/K http://en.wikipedia.org/wiki/Heat_capacity
# average density of granite is between 2.65 and 2.75 g/cm3 http://en.wikipedia.org/wiki/List_of_thermal_conductivities
# specific heat capacity is 790 J/kg/K
# process the physical params...
k=lambd/pc # thermal diffusivity. W/m/K
# convert to km and Ma to match our input data
SperMa=(60.0*60.0*24.0*365.0*1.0E6) # second per million year, for conversion
k=k *SperMa/1.0E6 # m^2/s *km/1000m *km/1000m *(60*60*24*365*1E6)/Ma
lambd=lambd *1000.0 # W/m/K=W/m/K =*1000m/km W/km/K
# sanity check
assert(np.round(k,1)==18.9) # check k is about 18.9, as in Gofery et al 2001
#==============================================================================
# # CALCULATION PARAMS
#==============================================================================
# depth ranges for generating and contouring (km)
zmin=0 # km
zmax=50 # km
zres=200 # calculation resolution in z
# emplacement time and depth, leave both at zero for most cases
z0=0 # depth km
t0=0 # emplace 13 Ma
# time range for calculating (in million years) (not zero as scale and calc are logarithmic)
tmin=0.001
tmax=40
tres=60 # calculation resolution in t
contour_res=15 # contour resolution
#==============================================================================
# # PLOT PARAMS
#==============================================================================
zlabel="Depth (km)"
tlabel="Time (Ma)"
title="Heatflow Increase (mW/m2)"
title2="Temperature Increase (K)"
rc('figure', figsize=(8.27*1.5/2,8.27*1.5/2)) # figure size etc for plotting
# MARKUP PARAMS
# Dates of dunedin volcanics
vstart=10.1 # start of valcanics in dunedin
vstop=16
#Depths
moho=25 # depth
d0=13.4
d1=19
fill_alpha=0.1 # how transparent the markup is
lw2=2 # linewidth
vmax=250 # max heat to contour
vmin=1 # min heat to contour
# heat flow max and min to mark on plot
hmax=25+9
hmin=25-9
#==============================================================================
# Functions
#==============================================================================
# function for heat flux
#@memoized # this remembers past results to speed up function for repeated calls
def q(r,t):
"""
eq 7.36 in Fowler solid earth pg 282, differentiating the simple 7.34 via r
r is distance from it, t is time since emplacement
also 7.5 in Beardsmore 's Crustal Heat flow pg 241
"""
if t<=0:
return 0
else:
return lambd*(T0-Tm)/sqrt(pi*t*k)*exp(-r**2/(4*k*t)) * 1E3/1E6# comes out in W/km^2 so 1km/1000m *1km/1000m also 1000 to make it mW/m^2
#@memoized
def dT(r,t):
"""
# this is in fowler solid earth page 282 eq 7.34
this is for a point source I think
r is distance from it, t is time since emplacement
"""
if t<=0:
return Tm
else:
return T0+(Tm-T0)*erf(r/sqrt(4*k*t))
#==============================================================================
# calculate
#==============================================================================
za=np.linspace(zmin,zmax,zres)
ta=np.linspace(tmin,tmax,tres)
t0i = ta.searchsorted(t0) # time sice clsoe to t0 Ma
z0i = za.searchsorted(z0) # time sice clsoe to t0 Ma
labels=(zlabel,tlabel)
inputs=(za,ta)
emplacment=(z0,t0)
slicei=(z0i,t0i)
# now populate the q and T matrixes
qr=np.zeros((len(za),len(ta)))
Tr=np.zeros((len(za),len(ta)))
for zi in xrange(len(za)):
for ti in xrange(len(ta)):
d=za[zi]
t=ta[ti]
qr[zi,ti]=q(d,t)
Tr[zi,ti]=dT(d,t)
#==============================================================================
# # diagram
#==============================================================================
#fig=plt.figure(0)
ax=plt.gca()
plt.ylabel("Emplacement depth (km)")
plt.xlabel("Emplacement time (Ma)")
plt.title(r"""Excess heat flow
for varying initial times and depths""")
#==============================================================================
# PLOTTING
#==============================================================================
cmap=cm.RdBu_r
levels = np.logspace(np.log10(vmin), np.log10(vmax), contour_res) # make logarithmic levels
norm=matplotlib.colors.LogNorm()
# plot filled contours
colors=plt.contourf(ta,za,qr[:,:],contour_res,norm=norm,cmap=cmap,levels=levels,zorder=0)
#==============================================================================
# Plot the colorbar and title
#==============================================================================
cbar2=plt.colorbar(colors, format='%.1f')
cbar2.ax.text(0,1.04,r"$mW/m^2$",fontsize=rcParams['font.size']*1.4) # this put vertical text in middle of cbar http://stackoverflow.com/questions/17475619/how-do-i-adjust-offset-colorbar-title-in-matplotlib
#==============================================================================
# mark observed excess heat
#==============================================================================
# mark extra surface heatflow
cmap=cm.Greys#cm.coolwarm#cm.rainbow
contour_res=2
levels = np.linspace(hmin, hmax, contour_res) # make logarithmic levels
norm=matplotlib.colors.Normalize()
contour1=plt.contour(ta,za,qr[:,:],contour_res,colors='black',levels=levels,zorder=0,alpha=1,linestyles="solid",linewidths=lw2)
contourf1=plt.contourf(ta,za,qr[:,:],contour_res,norm=norm,cmap=cmap,levels=levels,zorder=0,alpha=fill_alpha*2)
locx=22
conmin,segmin,imin,locx,locy,dmin=contour1.find_nearest_contour(locx,vstop,pixel=False)
ax.annotate('Observed \nexcess heat flow', xy=(locx,locy),
xytext=(30, 30), textcoords='offset points', va='center',ha='center',
arrowprops=dict(arrowstyle='->'))
cbar2.add_lines(contour1) # add lines to colorbar
#==============================================================================
# # MARK PEAK HEATFLOW
#==============================================================================
plt.plot(ta,np.sqrt(2*ta*k),'--',lw=lw2)
# label post surface and presurface
loc=30 # Ma # loc of labes
ax.annotate('Peak heat flow', xy=(loc,np.sqrt(2*loc*k)),
xytext=(-60, 40), textcoords='offset points', va='top',ha='left',
arrowprops=dict(arrowstyle='->'))
#==============================================================================
# # MARK volcanism times
#==============================================================================
plt.fill_between([vstart,vstop],[zmin,zmin],[zmax,zmax],color='black',alpha=fill_alpha,linestyle="dashdot",lw=0)
plt.vlines(vstop,zmin,zmax,color='black',alpha=1,linestyle="dashdot",lw=lw2)
plt.vlines(vstart,zmin,zmax,color='black',alpha=1,linestyle="dashdot",lw=lw2)
#plt.text(np.mean([vstart,vstop])-4,3,"Dunedin\nVolcanics",verticalalignment='center')
ax.annotate('Age of Dunedin \nvolcanics', xy=(vstop,3),
xytext=(20, -0), textcoords='offset points', va='center',
arrowprops=dict(arrowstyle='->'))
#==============================================================================
# # MARK 13.4 19 25 km special depths
#==============================================================================
plt.fill_betweenx([d0,moho],[tmin,tmin],[tmax,tmax],color='black',alpha=fill_alpha,lw=0)
locy=d0
locx=30
ax.annotate('Depth of \nreflective crust', xy=(locx,locy),
xytext=(-0, -30), textcoords='offset points', va='center',ha='center',
arrowprops=dict(arrowstyle='->'))
locx=37
locy=d0
plt.hlines(locy,tmin,tmax,color='black',linestyle="dotted",lw=lw2)
locy=moho
plt.hlines(locy,tmin,tmax,color='black',linestyle="dotted",lw=lw2)
vmean=np.mean([vstart,vstop])
# stats
print "Threshold depth (km) for ", vmean, "Ma is ", np.sqrt(vmean*18.9)
print "Threshold depth (km) for ", vstart, "Ma is ", np.sqrt(vstart*18.9)
print "Threshold depth (km) for ", vstop, "Ma is ", np.sqrt(vstop*18.9)
# heat stats
print "At t={} Ma and d={} km, excess heatflow is {}".format(vstop,d0,q(d0,vstop))
print "At t={} Ma and d={} km, excess heatflow is {}".format(vstop,moho,q(moho,vstop))
print "At t={} Ma and d={} km, excess heatflow is {}".format(vstart,moho,q(moho,vstart))
# print
plot_name="dq_for_t0_vs_z0"
#plt.subplots_adjust()
plt.savefig(plot_name+'.pdf',bbox_inches='tight')
plt.savefig(plot_name+'.png',dpi=450,bbox_inches='tight')
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment