Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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