In [1]:
import numpy as np
import pandas as pd
import vG_conv as vG
In [2]:
#get soil data for plot
soilsamples=pd.read_csv('soilsamplesVI.dat')
soilsamples
Out[2]:
In [3]:
#van Genuchten retention functions with fixed m
#def theta_f(psi, theta_res, theta_sat, alpha, n):
# m=1.-1./n
# theta=theta_res+(theta_sat-theta_res)*(1./(1.+(alpha*psi)**n))**m
# return theta
#def psi_f(theta, theta_res, theta_sat, alpha, n):
# m=1.-1./n
# psi=1./alpha * (1. / ((theta-theta_res)/(theta_sat-theta_res))**(1./m)-1.)**(1./n)
# return psi
In [4]:
def theta_ff(psi,sample,soilsamples):
return vG.theta_psi(psi, soilsamples.theta_sat[sample], soilsamples.theta_res[sample], soilsamples.alpha[sample], soilsamples.n[sample])
def psi_ff(theta,sample,soilsamples):
return vG.psi_theta(theta, soilsamples.theta_sat[sample], soilsamples.theta_res[sample], soilsamples.alpha[sample], soilsamples.n[sample])
In [5]:
psi = 10**np.arange(-3,8,0.1)
In [6]:
for i in soilsamples.index:
plot(theta_ff(psi,i,soilsamples),psi,label=''.join(['soil',str(i)]))
yscale('log')
legend(loc=1,ncol=2)
title('Retention curves')
xlabel('Theta [m3/m3]')
ylabel('psi [m]')
Out[6]:
In [7]:
#round depth to integrate close values to one layer
soilsamples['rdepth']=soilsamples.depth.round(decimals=1)
#group by depth
soil_zmean=soilsamples.groupby(by='rdepth', axis=0).mean()
for i in soil_zmean.index:
plot(theta_ff(psi,i,soil_zmean),psi,label=''.join([str(i),' m']))
yscale('log')
legend(loc=1,ncol=2)
title('Retention curves')
xlabel('Theta [m3/m3]')
ylabel('psi [m]')
Out[7]:
In [8]:
soil_zmean['rdepth']=soil_zmean.index
soil_zmean.index=np.arange(len(soil_zmean))
#rename columns according to EchoRD
dummy=soil_zmean[['Ks_meas','theta_sat','theta_res','alpha','n','BD','rdepth']].copy()
dummy.columns=['ks','ts','tr','alpha','n','roh','z']
dummy.index=dummy.index+1
dummy.z=-1.*dummy.z
#dummy.alpha=10000.*dummy.alpha
dummy.roh=dummy.roh*1000.
dummy.tr[dummy.tr<0.001]=0.001
#save to file
dummy.to_csv('matrix_VI.dat',float_format='%2.6f',sep=' ',index_label='no')
In [9]:
dummy
Out[9]:
In [10]:
#define soil type representation in matrixdef.dat manually - room for improvement
In [10]:
In [11]:
import scipy as sp
import matplotlib.pyplot as plt
import scipy.constants as const
import mcini_VI as mc
import dataread as dr
In [12]:
mc=dr.dataread_caos(mc)
In [13]:
import mcpickle as mcp
mcp.mcpick_in(mc,'mc_dump_VI.pickle')
In [13]:
In [ ]:
In [ ]:
In [14]:
#get soil data for plot
soilsamples=pd.read_csv('soilsamplesXI.dat',na_values='NaN')
soilsamples
Out[14]:
In [15]:
#round depth to integrate close values to one layer
soilsamples['rdepth']=soilsamples.depth.round(decimals=1)
In [16]:
#group by depth
soil_zmean=soilsamples[0:14].groupby(by='rdepth', axis=0).mean()
for i in soil_zmean.index:
plot(theta_ff(psi,i,soil_zmean),psi,label=''.join([str(i),' m']))
yscale('log')
legend(loc=3,ncol=2)
title('Retention curves')
xlabel('Theta [m3/m3]')
ylabel('psi [m]')
Out[16]:
In [17]:
soil_zmean['rdepth']=soil_zmean.index
soil_zmean.index=np.arange(len(soil_zmean))
#rename columns according to EchoRD
dummy=soil_zmean[['Ks_meas','theta_sat','theta_res','alpha','n','BD','rdepth']]
dummy.columns=['ks','ts','tr','alpha','n','roh','z']
dummy.index=dummy.index+1
dummy.z=-dummy.z.values
dummy.alpha=dummy.alpha.values
dummy.roh=dummy.roh.values*1000
dummy.tr[dummy.tr<0.001]=0.001
#save to file
dummy.to_csv('matrix_XI.dat',float_format='%2.6f',sep=' ',index_label='no')
In [18]:
import mcini_XI as mcXI
mcXI=dr.dataread_caos(mcXI)
In [19]:
mcp.mcpick_in(mcXI,'mc_dump_XI.pickle')
In [24]:
D=np.empty((101,mc.soilmatrix.no.size))
Dcalc=np.empty(101)
ku=np.empty(101)
ku1=np.empty(101)
psi=np.empty(101)
psi1=np.empty(101)
theta=np.empty(101)
theta1=np.empty(101)
cH2O=np.empty(101)
i=1
for j in np.arange(101):
thetaS=float(j)/100.
ku1[j]=mc.soilmatrix.ks[i]*thetaS**0.5*(1.-(1.-thetaS**(1./mc.soilmatrix.m[i]))**mc.soilmatrix.m[i])**2.
psi[j]=vG.psi_theta(thetaS,mc.soilmatrix.ts[i],mc.soilmatrix.tr[i],mc.soilmatrix.alpha[i],mc.soilmatrix.n[i],mc.soilmatrix.m[i])
ku[j]=vG.ku_psi(psi[j], mc.soilmatrix.ks[i], mc.soilmatrix.alpha[i], mc.soilmatrix.n[i], mc.soilmatrix.m[i])
psi1[j]= -1 * ( (1-thetaS**(1./mc.soilmatrix.m[i])) / (thetaS**(1./mc.soilmatrix.m[i]) ))**(1./mc.soilmatrix.n[i]) / mc.soilmatrix.alpha[i]
psi[0]=-1.0e+11
psi1=psi/100. #convert to [m]
theta1[j]=thetaS*(mc.soilmatrix.ts[i]-mc.soilmatrix.tr[i])+mc.soilmatrix.tr[i]
theta[j]=vG.theta_thst(thetaS,mc.soilmatrix.ts[i],mc.soilmatrix.tr[i])
cH2O[j]=vG.c_psi(psi[j],mc.soilmatrix.ts[i],mc.soilmatrix.tr[i],mc.soilmatrix.alpha[i], mc.soilmatrix.n[i], mc.soilmatrix.m[i])
dummy=-mc.soilmatrix.m[i]*(1./(1.+np.abs((psi[j])*mc.soilmatrix.alpha[i])**mc.soilmatrix.n[i]))**(mc.soilmatrix.m[i]+1) *mc.soilmatrix.n[i]*(abs(psi[j])*mc.soilmatrix.alpha[i])**(mc.soilmatrix.n[i]-1.)*mc.soilmatrix.alpha[i]
cH2O[j]=-1.*(mc.soilmatrix.ts[i]-mc.soilmatrix.tr[i])*dummy
Dcalc[j]=vG.D_psi(psi[j],mc.soilmatrix.ks[i],mc.soilmatrix.ts[i],mc.soilmatrix.tr[i],mc.soilmatrix.alpha[i], mc.soilmatrix.n[i], mc.soilmatrix.m[i])
theta[theta<0.01]=0.01
DI=(ku[:-1]+(np.diff(ku)/2.))*np.diff(psi)/np.diff(theta)
DI[0]=0 #define very low diffusion at zero
D[:,i]=Dcalc#/10000. # convert cm2/s -> m2/s
D[-1,:]=D[-2,:]
In [27]:
plot(ku)
plot(ku1)
Out[27]:
In [ ]: