Define density and viscosity functions
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
################## Load NIST tables ###########################
###############################################################
nist30 = pd.read_csv("../data/nist1-50_0-450.out",header=None,delim_whitespace=True)
nist30.columns = ['pressure','temperature','density','drhodp','drhodt','enth/1000','dhdp/1000','dhdt/1000','viscosity','dmudp','dmudt']
###############################################################
##################### Density #################################
# https://towardsdatascience.com/machine-learning-polynomial-regression-with-python-5328e4e8a386
###############################################################
d = np.array(nist30.loc[(nist30['pressure']==30) & (nist30['temperature']<=300),'density'])
x = np.array(nist30.loc[(nist30['pressure']==30) & (nist30['temperature']<=300),'temperature'])
X = x.reshape(-1, 1)
den_func_poly = PolynomialFeatures(degree=4)
X_poly = den_func_poly.fit_transform(X)
den_func_pol = LinearRegression()
den_func_pol.fit(X_poly, d)
###############################################################
##################### Viscosity ###############################
###############################################################
mu = np.array(nist30.loc[(nist30['pressure']==30) & (nist30['temperature']<=300),'viscosity'])
vis_func_poly = PolynomialFeatures(degree=6)
X_poly = vis_func_poly.fit_transform(X)
vis_func_pol = LinearRegression()
vis_func_pol.fit(X_poly, mu)
def hydroStat(radius='radius',
height='height',
aq_thickness='aq_thickness',
volume_flux='volume_flux',
fraction_flux='fraction_flux',
gravity='gravity',
heat_flux='heat_flux',
sedTC='sedTC',
aqTC='aqTC',
perm='perm',
T_sf='T_sf',
spec_heat='spec_heat'):
###################################################################################
######### Define geometry #########################################################
###################################################################################
a_dis = radius**2
hb_dis = height + aq_thickness
###################################################################################
######### Calculate specfic discharge #############################################
###################################################################################
qSpec_dis = (volume_flux/1000)/a_dis/fraction_flux # specific discharge of discharge site (m/s)
###################################################################################
######### Calculate temperatures based on heat flux and seafloor temperatues ######
###################################################################################
Tsbi_dis = T_sf + (heat_flux/1000)*(height/sedTC) # Temperature at Sediment-Basement Interface in discharge zone (C)
Taq_base_dis = Tsbi_dis + (heat_flux/1000)*(aq_thickness/aqTC) # Temperature at base of upflow zone (C)
delT_dis = np.average([Tsbi_dis,Taq_base_dis]) # Temperature difference (C)
rhoSW_dis = den_func_pol.predict(den_func_poly.fit_transform([[Taq_base_dis]])) # Mean density of seawater (kg/m3)
beta = -spec_heat * rhoSW_dis * qSpec_dis * hb_dis/aqTC # Dimensionaless parameter from B&P
###################################################################################
############# Calculate temperature profile #######################################
###################################################################################
depth = np.linspace(0,hb_dis,hb_dis+1)
Tcalc_prof = T_sf + delT_dis * (np.exp((beta * depth) / hb_dis)-1)/(np.exp(beta)-1)
###################################################################################
######### Density and viscocity profile ###########################################
###################################################################################
X = Tcalc_prof.reshape(-1, 1)
DensityCalc = den_func_pol.predict(den_func_poly.fit_transform(X))
VisCalc = vis_func_pol.predict(vis_func_poly.fit_transform(X))
###################################################################################
### Calculate change in density and change in pressure ############################
###################################################################################
d_depth = np.diff(depth)
avgDens = np.convolve(DensityCalc, np.ones(2), 'valid') / 2
avgVisc = np.convolve(VisCalc, np.ones(2), 'valid') / 2
delDens = np.empty(d_depth.size)
delDensLoss = np.empty(d_depth.size)
last = 0
lastLoss = 0
for i in np.arange(0,len(d_depth)):
delDens[i] = last = last + d_depth[i] * gravity * avgDens[i]
delDensLoss[i] = lastLoss = lastLoss + qSpec_dis* d_depth[i] * gravity * avgVisc[i]/perm
delDens = np.insert(delDens,0,0)
delDensLoss = np.insert(delDensLoss,0,0)
return {'Depth':depth,'TempProfile':Tcalc_prof,'Density':DensityCalc,'DelDens':delDens,'DelDensLoss':delDensLoss}
radius = 1000 # radius of discharge site (m)
height = 500 # Height to seafloor above regional basement (m)
aq_thickness = 1000 # Thickness of lateral flow zone in regional basement (m)
volume_flux = -10 # volume flux (L/s)
fraction_flux = 0.5 # fraction of ediface which flow occurs
gravity = 9.81 # gravity in m/s2
heat_flux = 180 # Heat flux at discharge site (mW/m2)
sedTC = 1.3 # Thermal conductivity of sediment (W/m C)
aqTC = 1.6 # Thermal conductivity of basalt (W/m C)
perm = 1e-12 # Permeability of the up flow site (m2)
T_sf = 2 # Temperature at top of upflow zone
spec_heat = 4200 # Specific heat of water (J/kg C)
r = hydroStat(radius,height,aq_thickness,volume_flux,fraction_flux,gravity,heat_flux,sedTC,aqTC,perm,T_sf,spec_heat)
radius = 1000 # radius of discharge site (m)
height = 500 # Height to seafloor above regional basement (m)
aq_thickness = 1000 # Thickness of lateral flow zone in regional basement (m)
volume_flux = 10 # volume flux (L/s)
fraction_flux = 1 # fraction of ediface which flow occurs
gravity = 9.81 # gravity in m/s2
heat_flux = 180 # Heat flux at discharge site (mW/m2)
sedTC = 1.3 # Thermal conductivity of sediment (W/m C)
aqTC = 1.6 # Thermal conductivity of basalt (W/m C)
perm = 1e-10 # Permeability of the up flow site (m2)
T_sf = 2 # Temperature at top of upflow zone
spec_heat = 4200 # Specific heat of water (J/kg C)
d = hydroStat(radius,height,aq_thickness,volume_flux,fraction_flux,gravity,heat_flux,sedTC,aqTC,perm,T_sf,spec_heat)
plt.plot(d['TempProfile'],d['Depth'],color='red')
plt.plot(r['TempProfile'],r['Depth'],color='blue')
plt.xlabel('Temperature (ºC)')
plt.ylabel('Depth (msb)')
plt.gca().invert_yaxis()
plt.plot(d['Density'],d['Depth'],color='red')
plt.plot(r['Density'],r['Depth'],color='blue')
plt.xlabel('Density (kg/m$^3$)')
plt.ylabel('Depth (msb)')
plt.gca().invert_yaxis()
plt.plot(((r['DelDens']-d['DelDens'])/1000),d['Depth'],color='green')
plt.plot(((r['DelDens']-d['DelDens']-d['DelDensLoss']-r['DelDensLoss'])/1000),d['Depth'],color='purple')
plt.xlabel('$\Delta$ Pressure (kPa)')
plt.ylabel('Depth (msb)')
plt.legend(['No energy loss','With energy loss'])
plt.gca().invert_yaxis()
radius = 1000 # radius of discharge site (m)
height = 500 # Height to seafloor above regional basement (m)
aq_thickness = 1000 # Thickness of lateral flow zone in regional basement (m)
volume_flux = -0.1 # volume flux (L/s)
fraction_flux = 0.5 # fraction of ediface which flow occurs
gravity = 9.81 # gravity in m/s2
heat_flux = 180 # Heat flux at discharge site (mW/m2)
sedTC = 1.3 # Thermal conductivity of sediment (W/m C)
aqTC = 1.6 # Thermal conductivity of basalt (W/m C)
perm = 1e-12 # Permeability of the up flow site (m2)
T_sf = 2 # Temperature at top of upflow zone
spec_heat = 4200 # Specific heat of water (J/kg C)
r = {}
vel = [1e-2,1e-1,1,5,10,20,30,40]
for v in vel:
d = hydroStat(radius=radius,
height=height,
aq_thickness=aq_thickness,
volume_flux=v,
fraction_flux=fraction_flux,
gravity=gravity,
heat_flux=heat_flux,
sedTC=sedTC,
aqTC=aqTC,
perm=perm,
T_sf=T_sf,
spec_heat=spec_heat)
r.update({str(v):d})
# print(v)
vel = list(map(str, r.keys()))
# loop through inner_keys
for v in vel:
# create a list of values for inner key
y_axis_values = [v['TempProfile'] for v in r.values()]
x_axis_values =[v['Depth'] for v in r.values()]
x_axis_values[0]
Temp_df = pd.concat([pd.DataFrame(x_axis_values[0]),pd.DataFrame.from_dict(y_axis_values).transpose()],axis=1, ignore_index=True)
Temp_df
pd.DataFrame.from_dict(y_axis_values).transpose().plot()
