Model and graph heatflow resulting from an emplacement for varying time before present and depth below msl of emplacment

model_emplacement_heat.py

# -*- 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()