SWI - single layer

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 describes an inclined confined aquifer in the opposite direction. There is no interface, everything is salt. isource is 0 for the well.


In [38]:
%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

Paths


In [39]:
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'

Model input


In [40]:
ml = mf.Modflow(modelname, exe_name=exe_name, model_ws=workspace)

Test variables


In [41]:
tscale	= 1000
nstp 	= tscale * 100 			#[]
perlen 	= tscale * 365. * 200. 	#[d]
ssz 	= 0.2 			#[]
Q 		= 0.01 	 	#[m3/d]

Discretization data


In [42]:
nlay = 1
nrow = 1
ncol = 200#4
delr = 20.#[20.,25.,30.,35.]
delc = 1.

topL = -10.
topR = -410.
botL = topL-41.

TOP  = np.linspace(topL,topR,ncol).reshape((nrow,ncol))
BOTM = np.array(botL+TOP-topL).reshape((nlay,nrow,ncol))

save_head_every=1

BCN: WEL


In [43]:
lrcQ1 = np.array([(0, 0, 199, Q)]) #LRCQ, Q[m**3/d]

BCN: GHB


In [44]:
lrchc 					= np.zeros((30, 5))
lrchc[:, [0, 1, 3, 4]] 	= [0, 0, 0., 0.8 / 2.0 * delc]
lrchc[:, 2] 			= np.arange(0, 30)

SWI2


In [45]:
zini 	= np.hstack(( topL * np.ones(24), np.arange(topL, botL, -0.5), botL * np.ones(94)))[np.newaxis, :]
isource = np.zeros((1, 200), dtype=np.int)
isource[:, :30] = -2

Model objects


In [46]:
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=[40.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 [47]:
ml.write_input() #--write the model files

Run the model


In [48]:
m = ml.run_model(silent=True, report=True)

Read the output

only one head entry and one zeta entry in binary files


In [49]:
headfile = modelname + '.hds'
hdobj = fu.HeadFile(headfile)
head = hdobj.get_data(idx=0)

In [50]:
zetafile = modelname + '.zta'
zobj = fu.CellBudgetFile(zetafile)
zeta = zobj.get_data(idx=0, text='      ZETASRF  1')[0]

Plot


In [51]:
import gridobj as grd
gr = grd.gridobj(discret)

In [52]:
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([0.,0.,0.,50.]))
leg = ax.legend(loc='lower left', numpoints=1)
leg._drawFrame = False