notebook.community

Edit and run

FloPy

MODPATH example

This notebook demonstrates how to create and run forward and backward tracking with MODPATH. The notebooks also shows how to create subsets of pathline and endpoint information, plot MODPATH results on ModelMap objects, and export endpoints and pathlines as shapefiles.





In [1]:



    


%matplotlib inline
import sys
import shutil
import os
import glob
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import flopy

print(sys.version)
print('numpy version: {}'.format(np.__version__))
print('matplotlib version: {}'.format(mpl.__version__))
print('pandas version: {}'.format(pd.__version__))
print('flopy version: {}'.format(flopy.__version__))



    


















    






3.6.0 |Anaconda 4.3.0 (x86_64)| (default, Dec 23 2016, 13:19:00) 
[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]
numpy version: 1.11.3
matplotlib version: 2.0.0
pandas version: 0.19.2
flopy version: 3.2.6

Copy modflow datasets to scratch directory





In [2]:



    


mffiles = glob.glob('../data/mp6/EXAMPLE.*')
for f in mffiles:
    shutil.copy(f, 'temp/' + os.path.split(f)[1])

Load MODFLOW model





In [3]:



    


model_ws = 'temp/'
m = flopy.modflow.Modflow.load('EXAMPLE.nam', model_ws=model_ws)
m.get_package_list()



    


















    
Out[3]:








['DIS', 'BAS6', 'WEL', 'RIV', 'RCH', 'OC', 'PCG', 'LPF']

















In [4]:



    


nrow, ncol, nlay, nper = m.nrow_ncol_nlay_nper
nrow, ncol, nlay, nper



    


















    
Out[4]:








(25, 25, 5, 3)

















In [5]:



    


m.dis.steady.array



    


















    
Out[5]:








array([ True, False,  True], dtype=bool)

















In [6]:



    


m.write_input()



    











In [7]:



    


hdsfile = flopy.utils.HeadFile(model_ws + 'EXAMPLE.HED')
hdsfile.get_kstpkper()



    


















    
Out[7]:








[(0, 0),
 (0, 1),
 (1, 1),
 (2, 1),
 (3, 1),
 (4, 1),
 (5, 1),
 (6, 1),
 (7, 1),
 (8, 1),
 (9, 1),
 (0, 2)]

















In [8]:



    


hds = hdsfile.get_data(kstpkper=(0, 2))

Plot RIV bc and head results





In [9]:



    


plt.imshow(hds[4, :, :])
plt.colorbar();



    


















    
























In [10]:



    


fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=m, layer=4)
quadmesh = modelmap.plot_ibound()
linecollection = modelmap.plot_grid()
riv = modelmap.plot_bc('RIV', color='g', plotAll=True)
quadmesh = modelmap.plot_bc('WEL', kper=1, plotAll=True)
contour_set = modelmap.contour_array(hds, 
                                     levels=np.arange(np.min(hds),np.max(hds),0.5), colors='b')
plt.clabel(contour_set, inline=1, fontsize=14)



    


















    
Out[10]:








<a list of 5 text.Text objects>

Create forward particle tracking simulation where particles are released at the top of each cell in layer 1

  • specifying the recharge package in create_mpsim releases a single particle on iface=6 of each top cell
  • start the particles at begining of per 3, step 1, as in example 3 in MODPATH6 manual




In [11]:



    


mp = flopy.modpath.Modpath(modelname='ex6',
                           exe_name='mp6',
                           modflowmodel=m,
                           model_ws='temp',
                           dis_file=m.name+'.dis',
                           head_file=m.name+'.hed',
                           budget_file=m.name+'.bud')

mpb = flopy.modpath.ModpathBas(mp, hdry=m.lpf.hdry, laytyp=m.lpf.laytyp, ibound=1, prsity=0.1)

# start the particles at begining of per 3, step 1, as in example 3 in MODPATH6 manual
# (otherwise particles will all go to river)
sim = mp.create_mpsim(trackdir='forward', simtype='pathline', packages='RCH', start_time=(2, 0, 1.))
mp.write_input()

mp.run_model(silent=False)



    


















    






FloPy is using the following executable to run the model: /Users/JosephHughes/Development/bin/mp6
Processing basic data ...
Checking head file ...
Checking budget file and building index ...

Run particle tracking simulation ...
Processing Time Step     1 Period     3.  Time =  4.10000E+05                                                                       
Particle tracking complete. Writing endpoint file ...                                                                               
End of MODPATH simulation. Normal termination.










    
Out[11]:








(True, [])

Read in the endpoint file and plot particles that terminated in the well





In [12]:



    


epobj = flopy.utils.EndpointFile('temp/ex6.mpend')
well_epd = epobj.get_destination_endpoint_data(dest_cells=[(4, 12, 12)]) 
# returns record array of same form as epobj.get_all_data()



    











In [13]:



    


well_epd[0:2]



    


















    
Out[13]:








rec.array([ (50, 0, 2, 0.0, 29565.619140625, 1, 0, 2, 0, 6, 1, 0.5, 0.5, 1.0, 200.0, 9000.0, 339.12310791015625, 1, 4, 12, 12, 2, 1, 1.0, 0.917884886264801, 0.09755218774080276, 5200.0, 5167.15380859375, 9.755219459533691, b'rch'),
 (75, 0, 2, 0.0, 26106.58984375, 1, 0, 3, 0, 6, 1, 0.5, 0.5, 1.0, 200.0, 8600.0, 339.12030029296875, 1, 4, 12, 12, 4, 1, 0.7848777770996094, 1.0, 0.13873140513896942, 5113.951171875, 5200.0, 13.873140335083008, b'rch')], 
          dtype=[('particleid', '<i8'), ('particlegroup', '<i8'), ('status', '<i8'), ('initialtime', '<f4'), ('finaltime', '<f4'), ('initialgrid', '<i8'), ('k0', '<i8'), ('i0', '<i8'), ('j0', '<i8'), ('initialcellface', '<i8'), ('initialzone', '<i8'), ('xloc0', '<f4'), ('yloc0', '<f4'), ('zloc0', '<f4'), ('x0', '<f4'), ('y0', '<f4'), ('z0', '<f4'), ('finalgrid', '<i8'), ('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('finalcellface', '<i8'), ('finalzone', '<i8'), ('xloc', '<f4'), ('yloc', '<f4'), ('zloc', '<f4'), ('x', '<f4'), ('y', '<f4'), ('z', '<f4'), ('label', 'O')])

















In [14]:



    


fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=m, layer=2)
quadmesh = modelmap.plot_ibound()
linecollection = modelmap.plot_grid()
riv = modelmap.plot_bc('RIV', color='g', plotAll=True)
quadmesh = modelmap.plot_bc('WEL', kper=1, plotAll=True)
contour_set = modelmap.contour_array(hds, 
                                     levels=np.arange(np.min(hds),np.max(hds),0.5), colors='b')
plt.clabel(contour_set, inline=1, fontsize=14)
modelmap.plot_endpoint(well_epd, direction='starting', colorbar=True)



    


















    
Out[14]:








<matplotlib.collections.PathCollection at 0x1104591d0>

Write starting locations to a shapefile





In [15]:



    


epobj.write_shapefile(well_epd, direction='starting', shpname='temp/starting_locs.shp', sr=m.sr)



    


















    






wrote temp/starting_locs.shp

Read in the pathline file and subset to pathlines that terminated in the well





In [16]:



    


# make a selection of cells that terminate in the well cell = (4, 12, 12)
pthobj = flopy.utils.PathlineFile('temp/ex6.mppth')
well_pathlines = pthobj.get_destination_pathline_data(dest_cells=[(4, 12, 12)])
well_pathlines # this is a record array of same form as pthobj._data



    


















    
Out[16]:








rec.array([ (50, 0, 0, 12, 0.0, 200.0, 9000.0, 339.12310791015625, 0, 2, 0, 1, 0.5, 0.5, 1.0, 0),
 (50, 0, -1, 12, 1.0, 200.02200317382812, 8999.9970703125, 339.0539855957031, 0, 2, 0, 1, 0.5000550746917725, 0.49999260902404785, 0.9990000128746033, 1),
 (50, 0, 1, 12, 1.0, 200.02200317382812, 8999.9970703125, 339.0539855957031, 0, 2, 0, 1, 0.5000550746917725, 0.49999260902404785, 0.9990000128746033, 0),
 ...,
 (550, 0, -1, 12, 29544.23046875, 5332.0380859375, 4800.0, 10.603289604187012, 4, 12, 13, 1, 0.3300943076610565, 0.0, 0.10603290051221848, 1),
 (550, 0, -1, 12, 29566.94921875, 5200.0, 4832.8828125, 9.754863739013672, 4, 12, 12, 1, 1.0, 0.08220770210027695, 0.09754862636327744, 1),
 (550, 0, -1, 12, 29566.94921875, 5200.0, 4832.8828125, 9.754863739013672, 4, 12, 12, 1, 1.0, 0.08220770210027695, 0.09754862636327744, 1)], 
          dtype=[('particleid', '<i8'), ('particlegroup', '<i8'), ('timepointindex', '<i8'), ('cumulativetimestep', '<i8'), ('time', '<f4'), ('x', '<f4'), ('y', '<f4'), ('z', '<f4'), ('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('grid', '<i8'), ('xloc', '<f4'), ('yloc', '<f4'), ('zloc', '<f4'), ('linesegmentindex', '<i8')])

Plot the pathlines that terminate in the well and the starting locations of the particles





In [17]:



    


fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=m, layer=2)
quadmesh = modelmap.plot_ibound()
linecollection = modelmap.plot_grid()
riv = modelmap.plot_bc('RIV', color='g', plotAll=True)
quadmesh = modelmap.plot_bc('WEL', kper=1, plotAll=True)
contour_set = modelmap.contour_array(hds, 
                                     levels=np.arange(np.min(hds),np.max(hds),0.5), colors='b')
plt.clabel(contour_set, inline=1, fontsize=14)

modelmap.plot_endpoint(well_epd, direction='starting', colorbar=True)
# for now, each particle must be plotted individually 
#(plot_pathline() will plot a single line for recarray with multiple particles)
for pid in np.unique(well_pathlines.particleid):
    modelmap.plot_pathline(pthobj.get_data(pid), layer='all', colors='red');

Write pathlines to a shapefile





In [18]:



    


# one line feature per particle
pthobj.write_shapefile(well_pathlines,
                       direction='starting', shpname='temp/pathlines.shp')

# one line feature for each row in pathline file 
# (can be used to color lines by time or layer in a GIS)
pthobj.write_shapefile(well_pathlines, one_per_particle=False, 
                       shpname='temp/pathlines_1per.shp')



    


















    






Warning: no grid spacing or lower-left corner supplied. Setting the offset with xul, yul requires arguments for delr and delc. Origin will be set to zero.
wrote temp/pathlines.shp
Warning: no grid spacing or lower-left corner supplied. Setting the offset with xul, yul requires arguments for delr and delc. Origin will be set to zero.
wrote temp/pathlines_1per.shp

Replace WEL package with MNW2; create backward tracking simulation using particles released at MNW well





In [19]:



    


model_ws = 'temp/'
m2 = flopy.modflow.Modflow.load('EXAMPLE.nam', model_ws=model_ws, exe_name='mf2005')
m2.get_package_list()



    


















    
Out[19]:








['DIS', 'BAS6', 'WEL', 'RIV', 'RCH', 'OC', 'PCG', 'LPF']

















In [20]:



    


m2.nrow_ncol_nlay_nper



    


















    
Out[20]:








(25, 25, 5, 3)

















In [21]:



    


m2.wel.stress_period_data.data



    


















    
Out[21]:








{0: 0, 1: rec.array([(4, 12, 12, -150000.0, 0.0)], 
           dtype=[('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('flux', '<f4'), ('iface', '<f4')]), 2: rec.array([(4, 12, 12, -150000.0, 0.0)], 
           dtype=[('k', '<i8'), ('i', '<i8'), ('j', '<i8'), ('flux', '<f4'), ('iface', '<f4')])}

















In [22]:



    


node_data = np.array([(3, 12, 12, 'well1', 'skin', -1, 0, 0, 0, 1., 2., 5., 6.2),
                      (4, 12, 12, 'well1', 'skin', -1, 0, 0, 0, 0.5, 2., 5., 6.2)], 
                     dtype=[('k', np.int), ('i', np.int), ('j', np.int), 
                            ('wellid', np.object), ('losstype', np.object), 
                            ('pumploc', np.int), ('qlimit', np.int), 
                            ('ppflag', np.int), ('pumpcap', np.int), 
                            ('rw', np.float), ('rskin', np.float), 
                            ('kskin', np.float), ('zpump', np.float)]).view(np.recarray)

stress_period_data = {0: np.array([(0, 'well1', -150000.0)], dtype=[('per', np.int), ('wellid', np.object), 
                                                            ('qdes', np.float)])}



    











In [23]:



    


m2.name = 'Example_mnw'
m2.remove_package('WEL')
mnw2 = flopy.modflow.ModflowMnw2(model=m2, mnwmax=1,
                                 node_data=node_data, 
                                 stress_period_data=stress_period_data, 
                                 itmp=[1, -1, -1])
m2.get_package_list()



    


















    
Out[23]:








['DIS', 'BAS6', 'RIV', 'RCH', 'OC', 'PCG', 'LPF', 'MNW2']

Write and run MODFLOW





In [24]:



    


m2.write_input()

m2.run_model(silent=False)



    


















    






FloPy is using the following executable to run the model: /Users/JosephHughes/USGS/Development/bin/mf2005

                                  MODFLOW-2005     
    U.S. GEOLOGICAL SURVEY MODULAR FINITE-DIFFERENCE GROUND-WATER FLOW MODEL
                             Version 1.12.00 2/3/2017                        

 Using NAME file: Example_mnw.nam 
 Run start date and time (yyyy/mm/dd hh:mm:ss): 2017/03/18 19:48:04

 Solving:  Stress period:     1    Time step:     1    Ground-Water Flow Eqn.
 Solving:  Stress period:     2    Time step:     1    Ground-Water Flow Eqn.
 Solving:  Stress period:     2    Time step:     2    Ground-Water Flow Eqn.
 Solving:  Stress period:     2    Time step:     3    Ground-Water Flow Eqn.
 Solving:  Stress period:     2    Time step:     4    Ground-Water Flow Eqn.
 Solving:  Stress period:     2    Time step:     5    Ground-Water Flow Eqn.
 Solving:  Stress period:     2    Time step:     6    Ground-Water Flow Eqn.
 Solving:  Stress period:     2    Time step:     7    Ground-Water Flow Eqn.
 Solving:  Stress period:     2    Time step:     8    Ground-Water Flow Eqn.
 Solving:  Stress period:     2    Time step:     9    Ground-Water Flow Eqn.
 Solving:  Stress period:     2    Time step:    10    Ground-Water Flow Eqn.
 Solving:  Stress period:     3    Time step:     1    Ground-Water Flow Eqn.
 Run end date and time (yyyy/mm/dd hh:mm:ss): 2017/03/18 19:48:04
 Elapsed run time:  0.098 Seconds

  Normal termination of simulation










    
Out[24]:








(True, [])

Create new MODPATH object





In [25]:



    


mp = flopy.modpath.Modpath(modelname='ex6mnw',
                           exe_name='mp6',
                           modflowmodel=m2,
                           model_ws='temp',
                           dis_file=m.name+'.dis',
                           head_file=m.name+'.hds',
                           budget_file=m.name+'.cbc')

mpb = flopy.modpath.ModpathBas(mp, hdry=m2.lpf.hdry, laytyp=m2.lpf.laytyp, ibound=1, prsity=0.1)
sim = mp.create_mpsim(trackdir='backward', simtype='pathline', packages='MNW2')

mp.write_input()

mp.run_model(silent=False)



    


















    






FloPy is using the following executable to run the model: /Users/JosephHughes/Development/bin/mp6
Processing basic data ...
Checking head file ...
Checking budget file and building index ...

Run particle tracking simulation ...
Processing Time Step     1 Period     1.  Time =  2.00000E+05                                                                       
Particle tracking complete. Writing endpoint file ...                                                                               
End of MODPATH simulation. Normal termination.
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL










    
Out[25]:








(True, [])

Read in results and plot





In [26]:



    


pthobj = flopy.utils.PathlineFile('temp/ex6mnw.mppth')
epdobj = flopy.utils.EndpointFile('temp/ex6mnw.mpend')
well_epd = epdobj.get_alldata()
well_pathlines = pthobj.get_alldata() # returns a list of recarrays; one per pathline



    











In [27]:



    


fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(1, 1, 1, aspect='equal')
modelmap = flopy.plot.ModelMap(model=m2, layer=2)
quadmesh = modelmap.plot_ibound()
linecollection = modelmap.plot_grid()
riv = modelmap.plot_bc('RIV', color='g', plotAll=True)
quadmesh = modelmap.plot_bc('MNW2', kper=1, plotAll=True)
contour_set = modelmap.contour_array(hds, 
                                     levels=np.arange(np.min(hds),np.max(hds),0.5), colors='b')
plt.clabel(contour_set, inline=1, fontsize=14)

modelmap.plot_pathline(well_pathlines, travel_time='<10000',
                       layer='all', colors='red');



    


















    
























In [ ]:

Content source: brclark-usgs/flopy

Similar notebooks:

notebook.community | gallery | about