Reproduction of the Zaidel Model


In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import os
import sys
import shutil

In [2]:
#import flopy
flopypath = '/Users/langevin/langevin/dev/flopy_site.svn/flopy_site/trunk'
if flopypath not in sys.path:
    sys.path.append(flopypath)
import flopy

In [3]:
#setup model parameters
nrow = 1
ncol = 200
nlay = 1
delr = 5.
delc = 1.
h1 = 23.
h2 = 1.

#ibound
ibound = np.ones((nlay, nrow, ncol), dtype=np.int)
ibound[:, :, 0] = -1
ibound[:, :, -1] = -1

#botm
botm = 25 * np.ones((nlay + 1, nrow, ncol), dtype=np.float)
base = 20.
for j in xrange(ncol):
    botm[1, :, j] = base
    if j > 0 and j % 40 == 0:
        base -= 5

#strt
strt = h1 * np.ones((nlay, nrow, ncol), dtype=np.float)
strt[:, :, -1] = h2

In [4]:
#make the flopy model
mf = flopy.modflow.Modflow(modelname='zaidel', exe_name='mfnwt')
dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, top=botm[0, :, :], botm=botm[1:, :, :], perlen=1, nstp=1, steady=True)
bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)
lpf = flopy.modflow.ModflowUpw(mf, hk=0.0001, laytyp=1)
oc = flopy.modflow.ModflowOc(mf, words=['pbudget', 'phead', 'head', 'budget'], save_head_every=1)
nwt = flopy.modflow.ModflowNwt(mf, linmeth=2, iprnwt=1, options='COMPLEX')
mf.write_input()

In [5]:
#shutil.copyfile('simple.nwt', 'zaidel.nwt')
#help(flopy.modflow.ModflowNwt)

In [6]:
#help(flopy.modflow.ModflowNwt)
try:
    os.remove('zaidel.hds')
except:
    pass
mf.run_model2()


Out[6]:
[True, []]

In [7]:
import flopy.utils.binaryfile as bf

#Create the headfile object
headfile = os.path.join('zaidel.hds')
headobj = bf.HeadFile(headfile)
times = headobj.get_times()
print times
head = headobj.get_data(totim=times[-1])


[1.0]

In [8]:
plt.plot(head[0, 0, :])
plt.plot(botm[1, 0, :])


Out[8]:
[<matplotlib.lines.Line2D at 0x1052e9a10>]

In [8]:


In [8]: