In [1]:
###%pylab inline
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm
from matplotlib import colors as ml_colors
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from mpl_toolkits.mplot3d import Axes3D
import os,sys
import scipy as sp
sys.path.append('/home/ecoon/research/python')
#import colors
In [2]:
sys.path.append('/home/ecoon/research/arctic/ats-python/src/models')
import ewc
import wrm_vangenuchten
import permafrost_model_explicit_fpd
In [3]:
class TC(object):
def thermal_conductivity(self, poro, sl, si, T):
raise RuntimeError("class TC is not implemented")
def plot_end_members(self, poro, s_r, axs=None, color='k', label=None, npoints=1000):
if axs is None:
fig = plt.figure(figsize=(4,8))
axs = []
axs.append(fig.add_subplot(311))
axs.append(fig.add_subplot(312))
axs.append(fig.add_subplot(313))
# unfrozen
sle = np.linspace(0,1,npoints)
sl = sle*(1-s_r) + s_r
si = 0
tc_u = np.array([self.thermal_conductivity(poro, s, si, 275.) for s in sl])
axs[0].plot(sle, tc_u, color=color, label=label)
axs[0].set_xlabel("1 - reduced gas saturation")
axs[0].set_title("T > 0")
axs[0].set_ylabel("thermal conductivity")
# frozen
sie = np.linspace(0,1,1000)
si = sie*(1-s_r)
tc_f = np.array([self.thermal_conductivity(poro, s_r, s, 260.0) for s in si])
axs[1].plot(sie, tc_f, color=color, label=label)
axs[1].set_xlabel("1 - reduced gas saturation")
axs[1].set_title("T << 0")
axs[1].set_ylabel("thermal conductivity")
# transition
sie = np.linspace(0,1,1000)
si = sie*(1-s_r)
tc_t = np.array([self.thermal_conductivity(poro, 1-s, s, 260.0) for s in si])
axs[2].plot(sie, tc_t, color=color, label=label)
axs[2].set_xlabel("reduced ice saturation")
axs[2].set_ylabel("thermal conductivity")
axs[2].set_title("sg = 0")
def plot_pcolor(self, poro, s_r, ax=None):
sle = np.linspace(0,1,100)
sie = np.linspace(0,1,100)
sl = sle*(1-s_r) + s_r
si = sie*(1-s_r)
TC = np.array([[self.thermal_conductivity(poro, k,j, 260.0) for k in sl] for j in si])
if ax is None:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(TC, origin='lower')
ax.set_xlabel("reduced liquid saturation")
ax.set_ylabel("reduced ice saturation")
return TC
def plot_surf(self, poro, s_r, fig, ax=None):
sle = np.linspace(0,1,100)
sie = np.linspace(0,1,100)
sl = sle*(1-s_r) + s_r
si = sie*(1-s_r)
TC = np.array([[self.thermal_conductivity(poro, k,j, 260.0) for k in sl] for j in si])
if ax is None:
ax = fig.add_subplot(111,projection="3d")
cmap = cm.jet
vmin = np.min(np.where(np.isnan(TC),1.e10,TC))
vmax = np.max(np.where(np.isnan(TC),-1.e10,TC))
norml = ml_colors.Normalize(vmin=vmin, vmax=vmax, clip = True)
SLE,SIE = np.meshgrid(sle,sie)
surf = ax.plot_surface(SLE, SIE, TC, rstride=1, cstride=1, cmap=cmap, norm=norml,
linewidth=0, antialiased=False)
ax.set_xlabel("reduced liquid saturation")
ax.set_ylabel("reduced ice saturation")
ax.set_zlabel("thermal conductivity")
fig.colorbar(surf,shrink=0.5,aspect=5)
return TC, ax
In [4]:
class PetersLidard(TC):
d_ = 0.053
eps_ = 1.e-10
def __init__(self, a_u, a_f, k_soil, k_gas=0.024, k_liq=0.5611, k_ice=2.27742):
self.a_u = a_u
self.a_f = a_f
self.k_soil = k_soil
self.k_gas = k_gas
self.k_liq = k_liq
self.k_ice = k_ice
def thermal_conductivity(self, poro, sl, si, T):
if sl + si > 1: return np.nan
k_dry = (self.d_*(1-poro)*self.k_soil + self.k_gas*poro) / \
(self.d_*(1-poro) + poro)
k_sat_u = pow(self.k_soil,(1-poro)) * pow(self.k_liq,poro)
k_sat_f = pow(self.k_soil,(1-poro)) * pow(self.k_ice,poro)
kersten_u = pow(sl + self.eps_, self.a_u)
kersten_f = pow(si + self.eps_, self.a_f)
return kersten_f * k_sat_f + kersten_u * k_sat_u + (1.0 - kersten_f - kersten_u) * k_dry
In [ ]:
class ThreePhaseWetDry(TC):
def __init__(self, tc_sat_uf, tc_dry, a_f, a_u):
self.a_u = a_u
self.a_f = a_f
self.tc_sat_uf = tc_sat_uf
self.tc_dry = tc_dry
def thermal_conductivity(self, poro, sl, si, T):
Ki = 831.51 * std::pow(T, -1.0552)
Kl = 0.5611;
double k_sat_f = k_sat_u_ * std::pow(Ki/Kl, poro);
double kersten_u = std::pow(sat_liq + eps_, alpha_u_);
double kersten_f = std::pow(sat_ice + eps_, alpha_f_);
return kersten_f * k_sat_f + kersten_u * k_sat_u_ + (1.0 - kersten_f - kersten_u) * k_dry_;
In [5]:
class PetersLidard2(TC):
d_ = 0.053
eps_ = 1.e-10
def __init__(self, a_u, a_f, k_soil, k_gas=0.024, k_liq=0.5611, k_ice=2.27742):
self.a_u = a_u
self.a_f = a_f
self.k_soil = k_soil
self.k_gas = k_gas
self.k_liq = k_liq
self.k_ice = k_ice
def thermal_conductivity(self, poro, sl, si):
if sl + si > 1: return np.nan
k_dry = (self.d_*(1-poro)*self.k_soil + self.k_gas*poro) / \
(self.d_*(1-poro) + poro)
k_sat_u = pow(self.k_soil,(1-poro)) * pow(self.k_liq,poro)
k_sat_f = pow(self.k_soil,(1-poro)) * pow(self.k_ice,poro)
kersten_u = pow(sl + self.eps_, self.a_u)
kersten_f = pow(si + self.eps_, self.a_f)
k_f = kersten_f * k_sat_f + (1-kersten_f)*k_dry
k_u = kersten_u * k_sat_u + (1-kersten_u)*k_dry
uf = sl / (sl + si)
return uf*k_u + (1-uf)*k_f
In [6]:
class PetersLidard3(TC):
d_ = 0.053
eps_ = 1.e-10
def __init__(self, a_u, a_f, k_soil, k_gas=0.024, k_liq=0.5611, k_ice=2.27742):
self.a_u = a_u
self.a_f = a_f
self.k_soil = k_soil
self.k_gas = k_gas
self.k_liq = k_liq
self.k_ice = k_ice
def thermal_conductivity(self, poro, sl, si):
if sl + si > 1: return np.nan
sg = 1-sl-si
k_dry = (self.d_*(1-poro)*self.k_soil + self.k_gas*poro) / \
(self.d_*(1-poro) + poro)
k_sat_u = pow(self.k_soil,(1-poro)) * pow(self.k_liq,poro)
k_sat_f = pow(self.k_soil,(1-poro)) * pow(self.k_ice,poro)
kersten_u = pow(1-sg + self.eps_, self.a_u)
kersten_f = pow(1-sg + self.eps_, self.a_f)
k_f = kersten_f * k_sat_f + (1-kersten_f)*k_dry
k_u = kersten_u * k_sat_u + (1-kersten_u)*k_dry
uf = sl / (sl + si)
return uf*k_u + (1-uf)*k_f
In [7]:
fig = plt.figure(figsize=(6,12))
axs = []
axs.append(fig.add_subplot(311))
axs.append(fig.add_subplot(312))
axs.append(fig.add_subplot(313))
tc_peat = PetersLidard(0.397947, 2.00553, 0.107283)
#tc_peat.plot_end_members(0.876418, 0.345733, axs, color='b')
#tc_peat.plot_end_members(0.876418, 0., axs, color='k')
#tc_peat2 = PetersLidard2(0.5, 1.0, 0.107283)
#tc_peat2.plot_end_members(0.876418, 0.345733, axs, color='b')
#tc_peat2.plot_end_members(0.876418, 0., axs, color='b')
#tc_peat3 = PetersLidard3(0.5, 1.0, 0.107283)
#tc_peat3.plot_end_members(0.876418, 0.345733, axs, color='m')
#
#tc_peat3.plot_end_members(0.876418, 0., axs, color='r')
#plt.show()
In [8]:
fig = plt.figure(figsize=(6,12))
axs = []
axs.append(fig.add_subplot(211))
axs.append(fig.add_subplot(212))
tc_peat = PetersLidard(0.5, 1.0, 0.107283)
#tc_peat.plot_end_members(0.876418, 0.345733, axs, color='k')
tc_peat.plot_pcolor(0.876418, 0., ax=axs[0])
axs[0].set_title("PL - 1")
tc_peat3 = PetersLidard3(0.5, 1.0, 0.107283)
#tc_peat3.plot_end_members(0.876418, 0.345733, axs, color='r')
tc_peat3.plot_pcolor(0.876418, 0., ax=axs[1])
axs[1].set_title("PL - 3")
Out[8]:
In [10]:
fig = plt.figure()
#tc_peat = PetersLidard(0.397947, 2.00553, 0.107283)
tc_peat = PetersLidard(0.397947, 2.00553, 2.30283)
vals,ax = tc_peat.plot_surf(0.876418, 0.345733, fig)
tc_soil = PetersLidard(0.797355, 0.732066, 2.32081)
tc_soil.plot_surf(0.59611, 0.199097, fig, ax)
#fig = plt.figure()
#tc_peat3 = PetersLidard3(0.5, 1.0, 0.107283)
#tc_peat3.plot_end_members(0.876418, 0.345733, axs, color='r')
#tc_peat3.plot_surf(0.876418, 0., fig)
plt.show()
In [10]:
plt.show()
In [47]:
dat = np.loadtxt("sampleset.matk", skiprows=3)
#dat = dat[0:100]
print len(dat)
In [8]:
peat = dict(poro=1, sr=7, ksoil=9, af=11, au=12)
soil = dict(poro=2, sr=8, ksoil=10, af=13, au=14)
In [45]:
fig = plt.figure(figsize=(12,12))
axs = []
axs.append(fig.add_subplot(321))
axs.append(fig.add_subplot(323))
axs.append(fig.add_subplot(325))
axs.append(fig.add_subplot(322))
axs.append(fig.add_subplot(324))
axs.append(fig.add_subplot(326))
cfield = peat['ksoil']
cm = colors.cm_mapper(dat[:,cfield].min(), dat[:,cfield].max())
for i in range(len(dat)):
color = cm(dat[i,cfield])
tc_peat = PetersLidard(dat[i,peat['au']], dat[i,peat['af']], dat[i,peat['ksoil']])
tc_peat.plot_end_members(dat[i,peat['poro']], dat[i,peat['sr']], axs, color=cm(dat[i,cfield]), npoints=100)
#tc_peat.plot_end_members(dat[i,peat['poro']], 0., axs, color=color, npoints=100)
#tc_peat.plot_end_members(0., 0., axs, color=cm(dat[i,cfield]), npoints=100)
tc_soil = PetersLidard(dat[i,soil['au']], dat[i,soil['af']], dat[i,soil['ksoil']])
tc_soil.plot_end_members(dat[i,soil['poro']], dat[i,soil['sr']], axs[3:], color=cm(dat[i,cfield]), npoints=100)
#tc_soil.plot_end_members(dat[i,soil['poro']], 0., axs[3:], color=color, npoints=100)
plt.show()
In [30]:
for i in range(len(dat)):
tc_peat = PetersLidard(dat[i,peat['au']], dat[i,peat['af']], dat[i,peat['ksoil']])
tc_peat.thermal_conductivity(dat[i,peat['poro']], 0., 0.)
In [46]:
len(dat)
Out[46]:
In [ ]: