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 flow is increased such, that the 50% front reaches the bottom of the confined aquifer. From this point, the SWI solution is only defined within the aquifer and therefore remains at the bottom of the aquifer.



In [59]:

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

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

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

Test variables



In [62]:

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

Discretization data



In [63]:

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

topL = -9.
botL = -50.

save_head_every = 1

BCN: WEL



In [64]:

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

BCN: GHB



In [65]:

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

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

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

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

Run the model



In [69]:

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

Read the output

only one head entry and one zeta entry in binary files



In [70]:

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



In [71]:

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

Plot



In [72]:

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



In [73]:

    
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.axis([gr.cGr[0,0],gr.cGr[0,-1],gr.lGr[-1]-10.,gr.lGr[0]+2.])
leg = ax.legend(loc='lower left', numpoints=1)
leg._drawFrame = False