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 the interface does move down if placed above the BGH solution. Especially remarkable is the small dip in the heel.
Changing the initial zeta condition such that the first zini value (all the way on the left) reaches the top of the aquifer, lets the SWI solution follow the BGH solution. If you move this cell to the right the swi solution gets off more ( zeta is lower)
In [136]:
%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 [137]:
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 [138]:
ml = mf.Modflow(modelname, exe_name=exe_name, model_ws=workspace)
In [139]:
tscale = 1.
nstp = tscale * 1000. #[]
perlen = tscale * 365. * 200. #[d]
ssz = 0.2 #[]
Q = 0.001 #[m3/d]
In [140]:
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 [141]:
lrcQ1 = np.array([(0, 0, 199, Q)]) #LRCQ, Q[m**3/d]
In [142]:
lrchc = np.zeros((30, 5))
lrchc[:, [0, 1, 3, 4]] = [0, 0, -topL*0.025, 0.8 / 2.0 * delc]
lrchc[:, 2] = np.arange(0, 30)
#print lrchc, topL
In [143]:
#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
zini[0,0] = topL #<<<========================================================================================================<<
isource = np.zeros((1, ncol), dtype=np.int)
isource[:, :30] = -2
imax = (BOTM[0,0,-1]-BOTM[0,0,-2])/delr
In [144]:
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 [145]:
ml.write_input() #--write the model files
In [146]:
m = ml.run_model(silent=True, report=True)
only one head entry and one zeta entry in binary files
In [147]:
headfile = modelname + '.hds'
hdobj = fu.HeadFile(headfile)
head = hdobj.get_data(idx=0)
In [148]:
zetafile = modelname + '.zta'
zobj = fu.CellBudgetFile(zetafile)
zeta = zobj.get_data(idx=0, text=' ZETASRF 1')[0]
In [149]:
import gridobj as grd
gr = grd.gridobj(discret)
In [150]:
fig = plt.figure(figsize=(20, 8), dpi=300, facecolor='w', edgecolor='k')
ax = fig.add_subplot(111)
ax.plot(gr.cGr[0,:-1],TOP.squeeze(),drawstyle='steps-post')
ax.plot(gr.cGr[0,:-1],BOTM.squeeze(),drawstyle='steps-post')
ax.plot(gr.cGr[0,:-1],zini[0,:], drawstyle='steps-post',label='Initial')
ax.plot(gr.cGr[0,:-1],zeta[0,0,:],drawstyle='steps-post', 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,-120.,-0.5,-4.5]))
leg = ax.legend(loc='lower left', numpoints=1)
leg._drawFrame = False
In [152]:
plt.plot(gr.cGr[0,:-1],zini[0,:]-zeta[0,0,:])
Out[152]:
VOLUMETRIC SWI ZONE BUDGET FOR ENTIRE MODEL
AT END OF TIME STEP 1000 IN STRESS PERIOD 1
ZONE 1
CUMULATIVE VOLUMES L**3 RATES FOR THIS TIME STEP L**3/T
------------------ ------------------------
IN: IN:
--- ---
BOUNDARIES = 73.0006 BOUNDARIES = 1.0000E-03
CONSTANT HEAD = 0.0000 CONSTANT HEAD = 0.0000
ZONE CHANGE = 5.8592 ZONE CHANGE = 0.0000
ZONE CHG TIP/TOE = 4.7836E-03 ZONE CHG TIP/TOE = 0.0000
ZONE MIXING = 0.0000 ZONE MIXING = 0.0000
TOTAL IN = 78.8645 TOTAL IN = 1.0000E-03
OUT: OUT:
---- ----
BOUNDARIES = 78.8409 BOUNDARIES = 1.0000E-03
CONSTANT HEAD = 0.0000 CONSTANT HEAD = 0.0000
ZONE CHANGE = 2.6695E-02 ZONE CHANGE = 0.0000
ZONE CHG TIP/TOE = 4.7836E-03 ZONE CHG TIP/TOE = 0.0000
ZONE MIXING = 0.0000 ZONE MIXING = 0.0000
TOTAL OUT = 78.8724 TOTAL OUT = 1.0000E-03
IN - OUT = -7.8430E-03 IN - OUT = -3.3528E-08
PERCENT DISCREPANCY = -0.01 PERCENT DISCREPANCY = -0.00
VOLUMETRIC SWI ZONE BUDGET FOR ENTIRE MODEL
AT END OF TIME STEP 1000 IN STRESS PERIOD 1
ZONE 2
CUMULATIVE VOLUMES L**3 RATES FOR THIS TIME STEP L**3/T
------------------ ------------------------
IN: IN:
--- ---
BOUNDARIES = 5.8406 BOUNDARIES = 3.4947E-08
CONSTANT HEAD = 0.0000 CONSTANT HEAD = 0.0000
ZONE CHANGE = 2.6311E-02 ZONE CHANGE = 0.0000
ZONE CHG TIP/TOE = 1.0994E-02 ZONE CHG TIP/TOE = 0.0000
ZONE MIXING = 0.0000 ZONE MIXING = 0.0000
TOTAL IN = 5.8779 TOTAL IN = 3.4947E-08
OUT: OUT:
---- ----
BOUNDARIES = 0.0000 BOUNDARIES = 0.0000
CONSTANT HEAD = 0.0000 CONSTANT HEAD = 0.0000
ZONE CHANGE = 5.8654 ZONE CHANGE = 0.0000
ZONE CHG TIP/TOE = 4.3995E-03 ZONE CHG TIP/TOE = 0.0000
ZONE MIXING = 0.0000 ZONE MIXING = 0.0000
TOTAL OUT = 5.8698 TOTAL OUT = 0.0000
IN - OUT = 8.0605E-03 IN - OUT = 3.4947E-08
PERCENT DISCREPANCY = 0.14 PERCENT DISCREPANCY = 200.00
HEAD WILL BE SAVED ON UNIT 51 AT END OF TIME STEP 1000, STRESS PERIOD 1 1
CUMULATIVE VOLUMES L**3 RATES FOR THIS TIME STEP L**3/T
------------------ ------------------------
IN: IN:
--- ---
STORAGE = 0.0000 STORAGE = 0.0000
CONSTANT HEAD = 0.0000 CONSTANT HEAD = 0.0000
WELLS = 73.0006 WELLS = 1.0000E-03
HEAD DEP BOUNDS = 5.8406 HEAD DEP BOUNDS = 3.4947E-08
SWIADDTOCH = 0.0000 SWIADDTOCH = 0.0000
TOTAL IN = 78.8411 TOTAL IN = 1.0000E-03
OUT: OUT:
---- ----
STORAGE = 0.0000 STORAGE = 0.0000
CONSTANT HEAD = 0.0000 CONSTANT HEAD = 0.0000
WELLS = 0.0000 WELLS = 0.0000
HEAD DEP BOUNDS = 78.8409 HEAD DEP BOUNDS = 1.0000E-03
SWIADDTOCH = 0.0000 SWIADDTOCH = 0.0000
TOTAL OUT = 78.8409 TOTAL OUT = 1.0000E-03
IN - OUT = 2.1362E-04 IN - OUT = 1.3970E-09
PERCENT DISCREPANCY = 0.00 PERCENT DISCREPANCY = 0.00
In [ ]: