This is the 1D testcase of the echoRD particle model foundation

As presented in section 4.1 in Jackisch and Zehe (WRR, 2014), this is the testcase to reproduce a plot scale sprinkler experiment with the 1D version of echoRD.

Make sure to have numpy, pandas, scipy and echoRD installed.

Preamble



In [1]:

    
import numpy as np
import pandas as pd



In [2]:

    
import scipy as sp
import matplotlib.pyplot as plt
import os, sys



In [3]:

    
#connect echoRD Tools
pathdir='../' #path to echoRD
lib_path = os.path.abspath(pathdir)
sys.path.append(lib_path)
import vG_conv as vG
from hydro_tools import hydroplot

Read Observations



In [4]:

    
obs=pd.read_csv('brspecht.dat',delimiter='\t')
obs.index=np.arange(-0.05,-1.,-0.1)

plt.figure(figsize=(18,6))
plt.subplot(131)
plt.imshow(obs.iloc[:,:10],cmap='Blues')
plt.colorbar()
plt.title('Br- Tracer Recovery\nProfile 1')
plt.subplot(132)
plt.imshow(obs.iloc[:,10:],cmap='Blues')
plt.title('Br- Tracer Recovery\nProfile 2')
plt.colorbar()
plt.subplot(133)
plt.plot(obs.iloc[:,:10].sum(axis=1),obs.index,label='Profile 1')
plt.plot(obs.iloc[:,10:].sum(axis=1),obs.index,label='Profile 2')
xlabel('c(Br-)')
ylabel('depth [m]')
plt.title('Br- Tracer Recovery\nLateral Sum over Depth')
plt.tight_layout()

Prepare echoRD



In [5]:

    
#connect to echoRD
import run_echoRD as rE
#connect and load project
[dr,mc,mcp,pdyn,cinf,vG]=rE.loadconnect(pathdir='../',mcinif='mcini_specht4y',oldvers=True)
mc.advectref='obs'
[mc,particles,npart,precTS]=rE.pickup_echoRD(mc, mcp, dr, 'specht4y.pickle')
mc.prects=True









    



updated mc from specht4y.pickle



In [6]:

    
#load 1D particle model
from echoRD_1D import gridupdate_thS1D,part_diffusion_1D,part_advect_1D,infilt_1Dobs,CAOSpy_run1D_adv
from hydro_tools import hydroprofile,ProgressBar



In [7]:

    
#initial moisture already set in prepared setup
[thS,npart]=pdyn.gridupdate_thS(particles.lat,particles.z,mc)



In [8]:

    
#as sprinkling we have a block rain with constant rate:
precTS









    Out[8]:






  
    
      
      tstart
      tend
      total
      intense
      conc
    
  
  
    
      0
       200
       4880
       0.021
       0.000004
       0.000165
    
  

1 rows × 5 columns



In [9]:

    
#prepare particles in 1D domain
thSx=thS[:,1]/100.
ths_part=mc.part_sizefac
npart=np.floor(mc.part_sizefac*thSx).astype(int)
#dgrid=np.append(-np.diff(mc.zgrid[:,1]),-np.diff(mc.zgrid[0:2,1]))
#npart=np.floor(dgrid*ths_part*vG.thst_theta(theta,mc.soilmatrix.ts[mc.soilgrid[:,1]-1], mc.soilmatrix.tr[mc.soilgrid[:,1]-1])).values.astype(int)
particles=pd.DataFrame(np.zeros(int(np.sum(npart))*9).reshape(int(np.sum(npart)),9),columns=['lat', 'z', 'conc', 'temp', 'age', 'flag', 'fastlane', 'advect','lastZ'])
particles['cell']=pd.Series(np.zeros(int(np.sum(npart)),dtype=int), index=particles.index)
k=0

cells=len(npart)

for i in np.arange(len(npart)):
    particles.cell[k:(k+npart[i])]=int(i)
    particles.lat[k:(k+npart[i])]=(np.random.rand(npart[i]))*mc.mgrid.latfac.values
    particles.z[k:(k+npart[i])]=(i+np.random.rand(npart[i]))*mc.mgrid.vertfac.values
    k+=npart[i]
    
thS=gridupdate_thS1D(particles.cell,mc,pdyn)



In [10]:

    
plt.figure(figsize=(6,6))
plt.plot(thS/100.,mc.zgrid[:,1])
plt.xlabel('Theta [m3/m3]')
plt.ylabel('depth [m]')
plt.title('Initial Soil Moistue')









    Out[10]:





<matplotlib.text.Text at 0x106b83b50>



In [26]:

    
#further preparations for 1D
mc.mgrid.latgrid=1
mc.mgrid.width=mc.mgrid.latfac*40.
drained=pd.DataFrame(np.array([]))
leftover=0
runname='echoRD_1D_sprinkler'



In [12]:

    
t_end=86400.
output=200.
dummy=np.floor(t_end/output)
dummi=dummy.astype(int)
leftover=0
drained=pd.DataFrame(np.array([]))
timenow=0.
p = ProgressBar(dummi)
mc.prects=False
ad_diff=False



In [ ]:

    
mc.soilmatrix



In [14]:

    
soil=mc.soilgrid[:,1]-1



In [15]:

    
zi=np.append(mc.zgrid[:,1],mc.soildepth)
TSstore_all=pd.DataFrame(np.zeros((dummi,len(zi))),columns=zi)
TSstore_new=pd.DataFrame(np.zeros((dummi,len(zi))),columns=zi)



In [16]:

    
particles.lastZ=particles.z

Run Lagrange Model



In [ ]:

    
for i in np.arange(dummi):
    particles.lastZ=particles.z
    [particles,thS,leftover,drained,timenow,infiltp] = CAOSpy_run1D_adv(particles,thS,leftover,drained,i*output,(i+1)*output,precTS,mc,pdyn,cinf,vG)
    
    oldp=np.bincount(np.round(np.append(-particles.z[particles.flag==0].values,-mc.zgrid[:,1])*100.).astype(int))-1
    allp=np.bincount(np.round(np.append(-particles.z[particles.flag<2].values,-mc.zgrid[:,1])*100.).astype(int))-1
    newp=allp[0:len(oldp)]-oldp
    TSstore_all.iloc[i]=allp
    TSstore_new.iloc[i]=newp
    p.animate(i+1)

Prepare and Compare Results and Observations



In [ ]:

    
sim_prof=TSstore_new.iloc[[9,36]].T
sim_prof.columns=['Model 30min','Model 120min']
obs_prof=pd.DataFrame([obs.iloc[:,:10].sum(axis=1)*300.,obs.iloc[:,10:].sum(axis=1)*300.]).T
obs_prof.columns=['Ref Profile for v','Obs Profile']
hydroprofile(sim_prof,obs_prof,fsize=(3, 6),xbound=[0.,70.,3],ybound=[-1.,0.,5],ptitle='Tracer Recovery [mg]\nModel vs. Observation\n',colors=['#87CEFA','#00BFFF','#FFA07A','#FF8C00'],saving='./results/Tracer_Recovery1.pdf')



In [151]:

    
sim_prof=TSstore_new.iloc[[90,180]].T
sim_prof.columns=['Model 5h','Model 10h']
obs_prof=pd.DataFrame([obs.iloc[:,:10].sum(axis=1)*300.,obs.iloc[:,10:].sum(axis=1)*300.]).T
obs_prof.columns=['Ref Profile for v','Obs Profile']
hydroprofile(sim_prof,obs_prof,fsize=(3, 6),xbound=[0.,70.,3],ybound=[-1.,0.,5],ptitle='\n\n\n',colors=['#87CEFA','#4169E1','#FFA07A','#FF8C00'],saving='./results/Tracer_Recovery2.pdf')
#hydroprofile(sim_prof,obs_prof,fsize=(3, 6),xbound=[0.,70.,3],ybound=[-1.,0.,5],ptitle='Tracer Recovery [mg]\nModel vs. Observation',colors=['#87CEFA','#4169E1','#FFA07A','#FF8C00'],saving='./results/Tracer_Recovery2.pdf')



In [154]: