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

The model is not approaching steady state just yet, therefore the SWI solution and the Badon Ghyben Herzberg cannot be compared. We will do this anyway.



In [89]:

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

    
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 [91]:

    
ml = mf.Modflow(modelname, exe_name=exe_name, model_ws=workspace)

Test variables



In [92]:

    
tscale	= 750
nstp 	= tscale * 100 			#[]
perlen 	= tscale * 365.			#[d]
ssz 	= 0.2 					#[]
Q 		= 0.005 	 			#[m3/d]

Discretization data



In [93]:

    
nlay = 1
nrow = 1
ncol = 200
delr = 20.
delc = 1.

topL = -9.
botL = -50.

save_head_every = 1

BCN: WEL



In [94]:

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

BCN: GHB



In [95]:

    
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 [96]:

    
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 [97]:

    
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=topL, botm=botL,
                        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 [98]:

    
ml.write_input() #--write the model files

Run the model



In [99]:

    
m = ml.run_model(silent=True, report=True)

Read the output

only one head entry and one zeta entry in binary files



In [100]:

    
headfile = modelname + '.hds'
hdobj = fu.HeadFile(headfile)
head = hdobj.get_data(idx=0)



In [101]:

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

Plot



In [102]:

    
import gridobj as grd
gr = grd.gridobj(discret)



In [103]:

    
fig = plt.figure()
ax = fig.add_subplot(111)
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],topL - 40. * (head[0, 0, :]), label='Ghyben-Herzberg')
ax.plot(gr.cGr[0,:-1],head[0, 0, :], label='Freshwater head')
ax.axis([gr.cGr[0,0],gr.cGr[0,-1],gr.lGr[-1]-2.,gr.lGr[0]+12.])
leg = ax.legend(loc='lower left', numpoints=1)
leg._drawFrame = False