A single confined layer. GHB on the LHS representing an ocean bcn. WEL on the RHS representing the freshwater coming from land. One stress period is defined of perlen days and is devided into nstp time steps and only the solution at the end of the stress period is written as output. bdestombe@gmail.com
This version shows an confined aquifer, freshwater heads are around 0m. The initial zeta surface is placed a little bit above the bottom of the single aquifer. The interface is at steady state (tscale 100 changed to 1000 gave no difference). The interface is not able to move upwards, but (with playing with zini) the interface does move down if placed above the BGH solution.
In [2]:
%matplotlib inline
import os
import sys
import numpy as np
import flopy.modflow as mf
import flopy.utils as fu
import matplotlib.pyplot as plt
In [3]:
os.chdir('C:\\Users\\Bas\\Google Drive\\USGS\\FloPy\\slope1D')
sys.path.append('C:\\Users\\Bas\\Google Drive\\USGS\\FloPy\\basScript') # location of gridObj
modelname = 'run1swi2'
exe_name = 'mf2005'
workspace = 'data'
In [4]:
ml = mf.Modflow(modelname, exe_name=exe_name, model_ws=workspace)
In [5]:
tscale = 1000
nstp = tscale * 100 #[]
perlen = tscale * 365. * 200. #[d]
ssz = 0.2 #[]
Q = 0.001 #[m3/d]
In [6]:
nlay = 1
nrow = 1
ncol = 200#4
delr = 1.#[20.,25.,30.,35.]
delc = 1.
topL = -20.
topR = -15.
botL = -21.
iGHB = np.arange(30)
TOP = np.linspace(topL,topL,ncol).reshape((nrow,ncol))
BOTM = np.array(botL+TOP-topL).reshape((nlay,nrow,ncol))
TOP[0,np.arange(iGHB.max(),ncol)] = np.linspace(topL,topR,ncol-iGHB.max())
BOTM = np.array(botL+TOP-topL).reshape((nlay,nrow,ncol))
save_head_every=1
In [7]:
lrcQ1 = np.array([(0, 0, 199, Q)]) #LRCQ, Q[m**3/d]
In [8]:
lrchc = np.zeros((30, 5))
lrchc[:, [0, 1, 3, 4]] = [0, 0, 0., 0.8 / 2.0 * delc]
lrchc[:, 2] = np.arange(0, 30)
In [9]:
#zini = np.hstack(( topL * np.ones(24), np.arange(topL, botL, -0.5), botL * np.ones(94)))[np.newaxis, :] # (nrow,ncol)
zini = botL*np.ones((1, ncol), dtype=np.float)+0.2
isource = np.zeros((1, ncol), dtype=np.int)
isource[:, :30] = -2
In [10]:
ml = mf.Modflow(modelname, version='mf2005', exe_name=exe_name)
discret = mf.ModflowDis(ml, nrow=nrow, ncol=ncol, nlay=nlay, delr=delr, delc=delc,
laycbd=[0], top=TOP, botm=BOTM,
nper=1, perlen=perlen, nstp=nstp)
bas = mf.ModflowBas(ml, ibound=1, strt=1.0)
bcf = mf.ModflowBcf(ml, laycon=[0], tran=[4.0])
wel = mf.ModflowWel(ml, stress_period_data={0:lrcQ1})
ghb = mf.ModflowGhb(ml, stress_period_data={0:lrchc})
swi = mf.ModflowSwi2(ml, nsrf=1, istrat=1, toeslope=0.02, tipslope=0.04, nu=[0, 0.025],
zeta=[zini], ssz=ssz, isource=isource, nsolver=1)
oc = mf.ModflowOc(ml, save_head_every=nstp)
pcg = mf.ModflowPcg(ml)
In [11]:
ml.write_input() #--write the model files
In [12]:
m = ml.run_model(silent=True, report=True)
only one head entry and one zeta entry in binary files
In [13]:
headfile = modelname + '.hds'
hdobj = fu.HeadFile(headfile)
head = hdobj.get_data(idx=0)
In [14]:
zetafile = modelname + '.zta'
zobj = fu.CellBudgetFile(zetafile)
zeta = zobj.get_data(idx=0, text=' ZETASRF 1')[0]
In [15]:
import gridobj as grd
gr = grd.gridobj(discret)
In [16]:
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(gr.cGr[0,:-1],TOP.squeeze(),drawstyle='steps')
ax.plot(gr.cGr[0,:-1],BOTM.squeeze(),drawstyle='steps')
ax.plot(gr.cGr[0,:-1],zini[0,:],drawstyle='steps', label='Initial')
ax.plot(gr.cGr[0,:-1],zeta[0,0,:],drawstyle='steps', label='SWI2')
ax.plot(gr.cGr[0,:-1],head[0, 0, :], label='feshw head')
ax.plot(gr.cGr[0,:-1],TOP[:,0] - 40. * (head[0, 0, :]), label='Ghyben-Herzberg')
ax.axis(gr.limLC())
leg = ax.legend(loc='lower left', numpoints=1)
leg._drawFrame = False
In [18]:
fig = plt.figure(figsize=(12, 6), dpi=300, facecolor='w', edgecolor='k')
ax = fig.add_subplot(111)
ax.plot(gr.cGr[0,:-1],TOP.squeeze(),drawstyle='steps')
ax.plot(gr.cGr[0,:-1],BOTM.squeeze(),drawstyle='steps')
ax.plot(gr.cGr[0,:-1],zini[0,:],drawstyle='steps', label='Initial')
ax.plot(gr.cGr[0,:-1],zeta[0,0,:],drawstyle='steps', label='SWI2')
ax.plot(gr.cGr[0,:-1],head[0, 0, :], label='feshw head')
ax.plot(gr.cGr[0,:-1],TOP[:,0] - 40. * (head[0, 0, :]), label='Ghyben-Herzberg')
ax.axis(gr.limLC([20,-160.,-0.5,-4.5]))
leg = ax.legend(loc='lower left', numpoints=1)
leg._drawFrame = False
In [17]: